oa_module.f

离散点插值程序

		!  Compute the angle from the positve x-direction to the line connecting the
		!  center of the curvature to a couple of different points - the observation and
		!  the (i,j) locations in the window.  The variable theta_ob is the same as the
		!  output from eqn 4.21 MH94, and theta_ij is eqn 4.22 from MH94.

		theta_ob = ATAN2 ( REAL(j_center)-yob     , REAL(i_center)-xob     )
		theta_ij = ATAN2 ( REAL(j_center)-REAL(j) , REAL(i_center)-REAL(i) )

		!  Compute the distance from the center of the curvature to the grid point
		!  (i,j).  This corresponds to the definition given for rij on the top of page 63
		!  in MH94.

		rij = SQRT ( (REAL(i_center)-REAL(i))**2 + (REAL(j_center)-REAL(j))**2 )

		!  All of the required pieces are assembled for use in MH94 eqns 4.15a, 
		!  4.14b and 4.14.  First, make sure that the values of theta_diff are bounded
		!  by +- 2 pi.  Then, bound the values between +- pi.

		theta_diff = theta_ob - theta_ij
		IF ( theta_diff .GT.  twopi ) theta_diff = theta_diff - REAL ( INT ( theta_diff/twopi ) ) * twopi
		IF ( theta_diff .LT. -twopi ) theta_diff = theta_diff + REAL ( INT ( theta_diff/twopi ) ) * twopi
		IF ( theta_diff .GT.     pi ) theta_diff = theta_diff - twopi
		IF ( theta_diff .LT.    -pi ) theta_diff = theta_diff + twopi

	!     x = radius_of_curvature * theta_diff
		x = radius_of_curvature * theta_diff / pi  ! this dividing by pi is home grown
		y = ABS(radius_of_curvature) - rij

		elongation = 1. + 0.7778 / vcrit * speed_ob
		distance = SQRT ( x*x/elongation + y*y )

		choice = 'banana  '

	   !  We have the fast enough speed (the observation speed is larger than the critical
	   !  speed), but too large radius of curvature (radius_curvature > 3X radius of influence).
	   !  This is, therefore, the ellipse scheme.

	   ELSE

		x = ( u_ob * ( REAL(i) - xob ) + v_ob * ( REAL(j) - yob ) ) / speed_ob
		y = ( v_ob * ( REAL(i) - xob ) - u_ob * ( REAL(j) - yob ) ) / speed_ob

		elongation = 1. + 0.7778 / vcrit * speed_ob
		distance = SQRT ( x*x/elongation + y*y )

		choice = 'ellipse '

	   END IF

	END SUBROUTINE dist

	!
	!---------------------------------------------------------------------------

	SUBROUTINE get_background_for_oa ( t , u , v , rh , slp_x , &
	dum2d , kp , name , &
	iew_alloc , jns_alloc , kbu_alloc )

	   INCLUDE 'first_guess_size.inc'
	   INCLUDE 'first_guess.inc'

	   REAL , DIMENSION ( iew_alloc , jns_alloc ) :: dum2d
	   INTEGER                                    :: kp
	   CHARACTER ( LEN = 8 )                      :: name


	   which_var_to_get : SELECT CASE ( name ) 

		CASE ( 'U       ' ) which_var_to_get
		   CALL yx2xy (     u(1,1,kp) , jns_alloc , iew_alloc , dum2d )

		CASE ( 'V       ' ) which_var_to_get
		   CALL yx2xy (     v(1,1,kp) , jns_alloc , iew_alloc , dum2d )

		CASE ( 'T       ' ) which_var_to_get
		   CALL yx2xy (     t(1,1,kp) , jns_alloc , iew_alloc , dum2d )

		CASE ( 'RH      ' ) which_var_to_get
		   CALL yx2xy (    rh(1,1,kp) , jns_alloc , iew_alloc , dum2d )

		CASE ( 'PSEALVLC' ) which_var_to_get
		   CALL yx2xy ( slp_x         , jns_alloc , iew_alloc , dum2d )

	   END SELECT which_var_to_get

	END SUBROUTINE get_background_for_oa

	!
	!---------------------------------------------------------------------------

	SUBROUTINE mqd ( obs , xob , yob , numobs , &
	gridded , iew , jns , &
	crsdot , name , passes , smooth_type , use_first_guess ) 

	!  Multiquadric interpolation based on Nuss and Titley (1994 MWR).
	!  The free parameters are lambda, the smoothing factor and c, the
	!  multiquadric parameter. If the procedure bombs, the most likely
	!  problem is that the value for c is incorrect. For long, narrow
	!  domains, it might be a problem. When there are lots of obs, this
	!  routine is a memory hog.

	   ! arguments:
	   !    obs:       difference between 1st guess and obs
	   !               ( need any duplicated obs removed from this array )
	   !    numobs:    number of station obs
	   !    xob,yob:   x, y locations of station obs
	   !    gridded:   background/first field as input, final analysis as output
	   !    iew,ins:   1st and 2nd dimensions of array 'gridded'
	   !    name:      variable name, T, RH, p, U, or V
	   !    use_first_guess: T/F use the first-guess fields in the objective analysis
    
	   ! other variables:
	   !    qi:        distance function for obs
	   !    qgi:       distance function for gridded data
	   !    dummy:     used for matrix calculation
	   !    workarray: to hold the input first guess array
	   !    errm:      error account for obs

	   IMPLICIT NONE

	   INTEGER, INTENT ( IN )                     :: numobs
	   INTEGER, INTENT ( IN )                     :: iew, jns
	   REAL,    INTENT ( IN ), DIMENSION ( numobs )      :: obs, xob, yob
	   REAL,    INTENT ( INOUT ), DIMENSION( iew , jns ) :: gridded
	   CHARACTER (8), INTENT ( IN )               :: name
	   INTEGER                                    :: crsdot
	   INTEGER , INTENT ( IN )                    :: passes , smooth_type 
	   LOGICAL , INTENT ( IN )                    :: use_first_guess

	   REAL, DIMENSION ( numobs, numobs )         :: qi
	   REAL                                       :: errm,           &
												   c
	   INTEGER                                    :: numpts,         &
												   i,  j,          &
												   ig, jg
	   REAL, DIMENSION ( iew,jns )                :: workarray2
	   REAL, PARAMETER                            :: lambda = 0.0025
	   INTEGER , DIMENSION ( numobs ) :: pvt
	   INTEGER , PARAMETER            :: job = 01
	   REAL    , DIMENSION ( numobs ) :: z , & 
									   work
	   REAL                           :: condition , det(2)
	   INTEGER                        :: imult , jmult , kmult


	   REAL, DIMENSION ( iew, numobs )            :: qgi
	   REAL, DIMENSION (numobs)                   :: dummy

	   INCLUDE 'error.inc'
	   INTERFACE 
		INCLUDE 'error.int'
	   END INTERFACE

	   c      = 0.01 * ( max ( iew, jns) )
	   numpts = iew * jns

	   !  Set errm, the mean error value for the variable being analyzed. 
	   !  The results aren't too sensitive to this parameter, but the 
	   !  values must be sane.

	   SELECT CASE ( name )
		CASE ( 'U       ' ) ; errm = 2.
		CASE ( 'V       ' ) ; errm = 2.
		CASE ( 'PSEALVLC' ) ; errm = 1.
		CASE ( 'T       ' ) ; errm = 0.5
		CASE ( 'RH      ' ) ; errm = 10.
	   END SELECT

	   !  Compute the radial basis function (how far is each observation
	   !  away from every other observation).  The diagonal entries are
	   !  identically zero, and are treated below.

	   do j = 1, numobs
		do i = 1, numobs
		   qi(i,j) = -1.* SQRT(((xob(j)-xob(i))**2 +               &
								(yob(j)-yob(i))**2)/(c*c)+1.)
		end do
	   end do

	   !  Account for observational uncertainty (Nuss and Titley 1994),
	   !  which are the diagonal entries.

	   do j = 1, numobs
		qi(j,j) = qi(j,j) + numobs * lambda * errm
	   end do

	   !  Find the inverse of Qi (radial basis function).  These are single
	   !  precision LINPAK routines.  Double precision routines are available.

	   CALL sgeco ( qi , numobs , numobs , pvt , condition , z )
	   CALL sgedi ( qi , numobs , numobs , pvt , det , work , job )

	   !  Compute the radial basis function for gridded points (similar to
	   !  how far away is each observation away from every grid point).

	   dummy = matmul(Qi, obs)

	   DO jg = 1, jns
		DO ig = 1, iew
		   DO i = 1, numobs
			  qgi(ig,i) = -1.* SQRT(((REAL(ig)-xob(i))**2 +          &
									(REAL(jg)-yob(i))**2)/(c*c)+1.)
		   END DO
		END DO
		workarray2(:,jg) = matmul(Qgi,dummy)
	   END DO


	   !  Which smoother will we use? Smooth the final analysis with a 5-point stencil 
	   !  (1-2-1 on the lateral boundaries), or the traditional smoother-desmoother.
	   !  Both accept a 2D field with the same arguments.

	   IF      ( smooth_type .EQ. 1 ) THEN
		CALL smooth_5            ( workarray2 , iew, jns, passes , crsdot ) 
	   ELSE IF ( smooth_type .EQ. 2 ) THEN
		CALL smoother_desmoother ( workarray2 , iew, jns, passes , crsdot )
	   END IF

	   !  Add the first guess to the perturbations from the observations at
	   !  each (i,j).

	   IF ( use_first_guess ) THEN
		gridded = gridded + workarray2
	   ELSE
		gridded = workarray2
	   END IF

	END SUBROUTINE mqd

	!
	!---------------------------------------------------------------------------

	SUBROUTINE put_background_from_oa ( t , u , v , rh , slp_x , &
	dum2d , kp , name , &
	iew_alloc , jns_alloc , kbu_alloc )

	   INCLUDE 'first_guess_size.inc'
	   INCLUDE 'first_guess.inc'

	   REAL , DIMENSION ( iew_alloc , jns_alloc ) :: dum2d
	   INTEGER                                    :: kp
	   CHARACTER ( LEN = 8 )                      :: name


	   which_var_to_get : SELECT CASE ( name ) 

		CASE ( 'U       ' ) which_var_to_get
		   CALL yx2xy ( dum2d , iew_alloc , jns_alloc ,     u(1,1,kp) ) 

		CASE ( 'V       ' ) which_var_to_get
		   CALL yx2xy ( dum2d , iew_alloc , jns_alloc ,     v(1,1,kp) ) 

		CASE ( 'T       ' ) which_var_to_get
		   CALL yx2xy ( dum2d , iew_alloc , jns_alloc ,     t(1,1,kp) ) 

		CASE ( 'RH      ' ) which_var_to_get
		   CALL yx2xy ( dum2d , iew_alloc , jns_alloc ,    rh(1,1,kp) ) 

		CASE ( 'PSEALVLC' ) which_var_to_get
		   CALL yx2xy ( dum2d , iew_alloc , jns_alloc , slp_x         ) 

	   END SELECT which_var_to_get

	END SUBROUTINE put_background_from_oa

	!
	!---------------------------------------------------------------------------

	SUBROUTINE smoother_desmoother ( slab , imx , jmx , passes , crsdot )

	   IMPLICIT NONE

	   INTEGER                        :: imx , jmx , passes , crsdot
	   REAL , DIMENSION ( imx , jmx ) :: slab , & 
									   slabnew

	   REAL , DIMENSION ( 2 )         :: xnu
	   INTEGER                        :: i , j , loop , n 

	   xnu  =  (/ 0.50 , -0.52 /)

	   !  The odd number passes of this are the "smoother", the even
	   !  number passes are the "de-smoother" (note the differnt signs on xnu).

	   smoothing_passes : DO loop = 1 , passes * 2

		n  =  2 - MOD ( loop , 2 )
 
		DO i = 2 , imx - 1 - crsdot
		   DO j = 2 , jmx - 1 - crsdot
			  slabnew(i,j) = slab(i,j) + xnu(n) *  & 
			  ((slab(i,j+1) + slab(i,j-1)) * 0.5-slab(i,j))
		   END DO
		END DO
 
		DO i = 2 , imx - 1 - crsdot
		   DO j = 2 , jmx - 1 - crsdot
			  slab(i,j) = slabnew(i,j)
		   END DO
		END DO
 
		DO j = 2 , jmx - 1 - crsdot
		   DO i = 2 , imx - 1 - crsdot
			  slabnew(i,j) = slab(i,j) + xnu(n) *  &
			  ((slab(i+1,j) + slab(i-1,j)) * 0.5-slab(i,j))
		   END DO
		END DO
 
		DO i = 2 , imx - 1 - crsdot
		   DO j = 2 , jmx - 1 - crsdot
			  slab(i,j) = slabnew(i,j)
		   END DO
		END DO
 
	   END DO smoothing_passes

	END SUBROUTINE smoother_desmoother

	!
	!---------------------------------------------------------------------------

	SUBROUTINE smooth_5 ( field , iew , jns , passes , crsdot )

	   INTEGER                        :: iew , jns , &
									   passes    , &
									   crsdot
	   REAL , DIMENSION ( iew , jns ) :: field

	   REAL , DIMENSION ( iew , jns ) :: temp
	   INTEGER                        :: i , j , num_passes

	   !  How may passes of this smoother are we using.

	   smoothing_passes : DO num_passes = 1 , passes

		!  Apply 5-point stencil smoother on interior of the domain.
   
		DO j = 2 , jns - 1 - crsdot
		   DO i = 2 , iew - 1 - crsdot
			  temp(i,j) = ( field(i  ,j  ) * 4. +  & 
							field(i+1,j  )      +  & 
							field(i-1,j  )      +  & 
							field(i  ,j+1)      +  & 
							field(i  ,j-1)      )  * 0.125
		   END DO
		END DO

		!  Apply 3-point stencil smoother on the boundaries.
   
		i = 1
		DO j = 2 , jns - 1 - crsdot
		   temp(i,j) = ( field(i  ,j  ) * 2. +  & 
						 field(i  ,j+1)      +  & 
						 field(i  ,j-1)      )  * 0.25
		END DO

		i = iew - crsdot
		DO j = 2 , jns - 1 - crsdot
		   temp(i,j) = ( field(i  ,j  ) * 2. +  & 
						 field(i  ,j+1)      +  & 
						 field(i  ,j-1)      )  * 0.25
		END DO
   
		j = 1
		DO i = 2 , iew - 1 - crsdot
		   temp(i,j) = ( field(i  ,j  ) * 2. +  & 
						 field(i+1,j  )      +  & 
						 field(i-1,j  )      ) * 0.25
		END DO
   
		j = jns - crsdot
		DO i = 2 , iew - 1 - crsdot
		   temp(i,j) = ( field(i  ,j  ) * 2. +  & 
						 field(i+1,j  )      +  & 
						 field(i-1,j  )      ) * 0.25
		END DO
   
		!  Store smoothed field back into original array.
   
		DO j = 2 , jns - 1 - crsdot
		   DO i = 2 , iew - 1 - crsdot
			  field(i,j) = temp(i,j)
		   END DO
		END DO
   
		!  Store smoothed boundary field back into original array.
   
		DO j = 2 , jns - 1 - crsdot
		   field(1         ,j) = temp(1         ,j)
		   field(iew-crsdot,j) = temp(iew-crsdot,j)
		END DO
   
		DO i = 2 , iew - 1 - crsdot
		   field(i,1         ) = temp(i,1         )
		   field(i,jns-crsdot) = temp(i,jns-crsdot)
		END DO

	   END DO smoothing_passes

	END SUBROUTINE smooth_5

	!---------------------------------------------------------------------------

	END MODULE obj_analysis