oa_module.f

离散点插值程序

	!客观分析程序
	MODULE obj_analysis

	   REAL , PARAMETER , PRIVATE :: pi = 3.14159265358
	   REAL , PARAMETER , PRIVATE :: twopi = 2. * 3.14159265358

	CONTAINS

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

	SUBROUTINE clean_rh ( rh , iew , jns , rh_min , rh_max )

	   INTEGER                        :: iew , jns
	   REAL , DIMENSION ( iew , jns ) :: rh
	   REAL                           :: rh_min , rh_max

	   WHERE ( rh .GT. rh_max ) rh = rh_max
	   WHERE ( rh .LT. rh_min ) rh = rh_min

	END SUBROUTINE clean_rh

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

	SUBROUTINE cressman ( obs_value , obs , xob , yob , numobs , &
	gridded , iew , jns , &
	crsdot , name , radius_influence , &
	dxd , u_banana , v_banana , pressure , passes , smooth_type , lat_center , &
	use_first_guess ) 

	   ! arguments:
	   !    obs_value: difference between 1st guess and obs
	   !               ( need any duplicated obs removed from this array )
	   !    obs:       difference of obs_value and the interpolated first_guess value
	   !    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
	   !    dxd        grid distance of this domain (m)
	   !    u_banana   first guess field of U for banana scheme
	   !    v_banana   first guess field of V for banana scheme
	   !    pressure   pressure of this level (Pa, 1001 signifies surface)
	   !    use_first_guess T/F to include first-guess in the objective analysis
    
	   IMPLICIT NONE

	   INTEGER, INTENT ( IN )                     :: numobs
	   INTEGER, INTENT ( IN )                     :: iew, jns
	   REAL,    INTENT ( IN ), DIMENSION ( numobs )      :: obs_value , obs, xob, yob
	   REAL,    INTENT ( INOUT ), DIMENSION( iew , jns ) :: gridded
	   INTEGER, INTENT ( IN )                     :: radius_influence
	   CHARACTER (8), INTENT ( IN )               :: name
	   REAL ,   INTENT ( IN ) , DIMENSION( iew , jns ) :: u_banana , v_banana
	   REAL ,   INTENT ( IN )                     :: dxd , pressure
	   INTEGER                                    :: crsdot
	   INTEGER , INTENT( IN )                     :: passes , smooth_type
	   REAL ,   INTENT ( IN )                     :: lat_center
	   LOGICAL , INTENT ( IN )                    :: use_first_guess

	   REAL, DIMENSION ( iew,jns )                :: denominator,    &
												   numerator ,     &
												   perturbation
	   REAL                                       :: dxob, dyob, distance , &
												   weight
	   INTEGER                                    :: iob, job
	   INTEGER                                    :: numpts,         &
												   i,  j , n ,     &
												   iew_min  ,      &
												   jns_min  ,      &
												   iew_max  ,      &
												   jns_max
	   CHARACTER (LEN=8)                          :: choice

	   INCLUDE 'error.inc'

	   INTERFACE 
		INCLUDE 'error.int'
	   END INTERFACE

	   !  Initialize the Cressman weights (numerator and denominator) and
	   !  the total perturbation field to zero.

	   numerator    = 0
	   denominator  = 0
	   perturbation = 0

	   !  Loop over each observation.

	   obs_loop : DO n = 1 , numobs

		!  Every observation has a first guess wind field location that will be specified
		!  as an objective analysis that is circular, elliptical or banana-shaped.

		choice = 'undefine'

		!  To use the observation, it must be far enough inside our domain
		!  (using the cross/dot configuration for the staggered variables).  If
		!  it is not, we simply zoom back to the top of this loop.

		IF ( ( xob(n) .LE.   1 + REAL(crsdot)/2. ) .OR. & 
			 ( yob(n) .LE.   1 + REAL(crsdot)/2. ) .OR. & 
			 ( xob(n) .GE. iew - REAL(crsdot)/2. ) .OR. & 
			 ( yob(n) .GE. jns - REAL(crsdot)/2. ) ) THEN
		   CYCLE obs_loop
		END IF

		!  Compute the window in which this observation has influence.  For
		!  some points inside this window the observation will have zero 
		!  influence, but this allows a direct method rather than a search.

		iew_min = MAX ( NINT ( (xob(n)-REAL(crsdot)/2.) ) - 2*radius_influence -1 ,            1 )
		jns_min = MAX ( NINT ( (yob(n)-REAL(crsdot)/2.) ) - 2*radius_influence -1 ,            1 )
		iew_max = MIN ( NINT ( (xob(n)-REAL(crsdot)/2.) ) + 2*radius_influence +1 , iew - crsdot ) 
		jns_max = MIN ( NINT ( (yob(n)-REAL(crsdot)/2.) ) + 2*radius_influence +1 , jns - crsdot )

		!  For the grid points in the window, compute this observation's
		!  influence.

		DO j = jns_min , jns_max
		   DO i = iew_min , iew_max

			  !  Distance from each of the windowed grid points to the observation.  The distance is in 
			  !  grid point units.  If this is a horizontal wind field or a moisture field, we can use the
			  !  banana scheme, which elongates the acceptable distance from an observation and also is
			  !  able to handle first-guess wind curvature.

			  IF ( ( name(1:8) .EQ. 'U       ' ) .OR. &
				   ( name(1:8) .EQ. 'V       ' ) .OR. &
				   ( name(1:8) .EQ. 'RH      ' )  ) THEN
				 CALL dist ( i , j , iew , jns , crsdot , name(1:8) , obs_value(n) , xob(n) , yob(n) , &
							 dxd , pressure , u_banana , v_banana , radius_influence , &
							 lat_center , choice , distance )
   
			  ELSE 
				 distance = SQRT ( ( ( xob(n) - REAL(crsdot)/2. - REAL(i) ) * ( xob(n) - REAL(crsdot)/2. - REAL(i) ) ) +  &
								   ( ( yob(n) - REAL(crsdot)/2. - REAL(j) ) * ( yob(n) - REAL(crsdot)/2. - REAL(j) ) ) ) 

			  END IF

			  !  If the observation is within a grid point's radius of influence,
			  !  sum its effect.

			  IF ( distance .LT. radius_influence ) THEN
				 weight = ( radius_influence**2 - distance**2 ) / ( radius_influence**2 + distance**2 )
				 numerator(i,j)   = numerator(i,j)   + weight * weight * obs(n)
				 denominator(i,j) = denominator(i,j) + weight
			  END IF

		   END DO
		END DO

	   END DO obs_loop

	   !  We have looped over all of the observations, and have computed the required
	   !  sums for the objective analysis for each grid point.  Compute the grid
	   !  point perturbations from these values.

	   WHERE ( denominator .GT. 0 ) 
		perturbation = numerator / denominator
	   ENDWHERE

	   !  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            ( perturbation , iew, jns, passes , crsdot ) 
	   ELSE IF ( smooth_type .EQ. 2 ) THEN
		CALL smoother_desmoother ( perturbation , 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 + perturbation
	   ELSE
		gridded = perturbation
	   END IF

	END SUBROUTINE cressman

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

	!  This routine computes the distance from an observation to a grid point.
	!  The Euclidean distance is modified to allow for elongation due to
	!  strong stream-wise flow, as well as for curvature.

	SUBROUTINE dist ( i , j , iew , jns , crsdot , name , obs_value , xob , yob , dxd , &
					pressure , u_banana , v_banana , radius_influence , &
					lat_center , choice , distance )
  
	   IMPLICIT NONE
 
	   !  Input variables.

	   INTEGER , INTENT(IN)                      :: i        , &
												  j        , &
												  iew      , &
												  jns      , &
												  crsdot   , &
												  radius_influence
	   REAL    , INTENT(IN)                      :: obs_value , &
												  xob      , &
												  yob      , &
												  dxd      , &
												  pressure , &
												  lat_center
	   REAL    , INTENT(IN) , DIMENSION(iew,jns) :: u_banana , &
												  v_banana
	   CHARACTER (LEN=8), INTENT(IN)             :: name

	   !  Output variables.

	   REAL    , INTENT(OUT)                     :: distance
	   CHARACTER (LEN=8) , INTENT(INOUT)         :: choice

	   !  Local variables.

	   REAL :: radius_of_curvature , k
	   REAL :: u_ob , dudx , speed_ob , dydx
	   REAL :: v_ob , dvdx
	   REAL :: numerator , denominator
	   REAL :: x1 , x2 , y1 , y2 , d12
	   INTEGER :: i1 , i2 , j1 , j2
	   INTEGER :: ia, ja, ib, jb, ic, jc, id, jd
	   REAL :: ua, va, ub, vb, uc, vc, ud, vd, vor
   

	   INTEGER :: i_center , j_center

	   REAL :: theta_ob , theta_ij
   
	   REAL :: rij
   
	   REAL :: vcrit , elongation

	   REAL :: x , y , theta_diff

	   SAVE

	   IF ( choice .EQ. 'undefine' ) THEN

		!  The critical speed to decide whether or not to do a simple
		!  circular Cressman or one of the anisotropic schemes.  If the 
		!  observation speed is < the critical speed, do a circular
		!  Cressman distance computation.  At 1000 hPa, the critical speed
		!  is 5 m/s; this goes linearly to 15 m/s at 500 hPa.
   
		IF ( pressure < 500. ) THEN
		   vcrit = 15.
		ELSE
		   vcrit = 25. - 0.02 * pressure
		END IF
   
		!  Compute the radius of curvature for the streamline passing through
		!  this observation location.  This is the information available from 
		!  eqn 4.16 in Manning and Haagenson, 1994 (MH94).  Note that MH94 eqn 4.16
		!  is different from that computation given below.
   
		!  dy/dx = slope of streamline through grid point, which is v/u
		!  d^2y/dx^2 = rate of change of slope of streamline through grid point, which
		!              is d(v/u)/dx if taken along stream line trajectory, which is
		!              u*dv/dx - v*du/dx
		!              -----------------
		!                    u^2
		!  k = curvature = abs(dy/dx) / ( 1 + (d^2y/dx^2)^2 ) ^(3/2)
		!  radius of curvature = 1 / k
   
		!  Values of u and v nearest the observation point in the first-guess data and the
		!  respective speed at this point.
   
		u_ob = u_banana(NINT(xob),NINT(yob))
		v_ob = v_banana(NINT(xob),NINT(yob))
		speed_ob = SQRT ( u_ob**2 + v_ob**2 ) 
   
		!  We can see if we need to do any fancy computations by just looking at the speed
		!  of the observation.  If the first-guess speed is less than the cutoff, we do a 
		!  simple circular Cressman and vamoose.
   
		IF ( speed_ob .GE. vcrit ) THEN
      
		   !  The dy/dx has to be defined, so u can't be = zero.  If it is, then we end up
		   !  doing an ellipse, not a banana, anyways.
   
		   IF ( ABS(u_ob) .GE. 1.E-5 ) THEN
			  dydx = v_ob / u_ob
		   ELSE
			  dydx = 1000.
		   END IF
   
		   !  Differences are computed from approximately 3 grid points on either side of the 
		   !  observation location, normal and tanget to the flow at the observation location.
    
		   x1 = xob + 3. * ( u_ob/speed_ob )
		   y1 = yob + 3. * ( v_ob/speed_ob )
		   x2 = xob - 3. * ( u_ob/speed_ob )
		   y2 = yob - 3. * ( v_ob/speed_ob )
   
		   !  These values need to be inside the domain since we will use them as indices
		   !  for the u and v first-guess winds.
   
		   i1 = MAX ( MIN ( NINT(x1) , iew ) , 1 )
		   j1 = MAX ( MIN ( NINT(y1) , jns ) , 1 )
		   i2 = MAX ( MIN ( NINT(x2) , iew ) , 1 )
		   j2 = MAX ( MIN ( NINT(y2) , jns ) , 1 )
   
		   !  And the distance between these two points (in grid point units) is 
   
		   d12 = SQRT( (REAL(i1) - REAL(i2))**2 + (REAL(j1) - REAL(j2))**2 )
   
		   !  So we can now compute the du/dx and the dv/dx, where we are assumming that
		   !  we are on a streamline through this observation point.  We extend the 
		   !  difference calculation in approximately 3 grid points to either side of
		   !  the observation point on the tangent to the streamline.
   
		   dudx = ( u_banana(i1,j1) - u_banana(i2,j2) ) / ( d12 * dxd ) 
		   dvdx = ( v_banana(i1,j1) - v_banana(i2,j2) ) / ( d12 * dxd ) 
      
		   !  Build the numerator and denominator of the curvature equation (k from above).
		   !  Since the radius of curvature is the inverse of the curvature, there is no
		   !  real need to do a temporary computation and then an inverse, so we just do it
		   !  directly (denominator / numerator).  The radius of curvature is now in units
		   !  of (m).
      
		   numerator = ABS ( u_ob * dvdx - v_ob * dudx ) / u_ob**2
		   denominator = ( SQRT ( 1. + dydx**2 ) ) ** 3
   
		   IF ( ( ABS ( u_ob) .LT. 1.E-5 ) .OR. ( ABS ( v_ob) .LT. 1.E-5 ) .OR. &
				( ABS ( numerator ) .LT. 1.E-5 ) ) THEN
			  radius_of_curvature = 1000.
		   ELSE
	   !        k = numerator / denominator
	   !        radius_of_curvature = 1. / k
			  radius_of_curvature = denominator / numerator
		   END IF
   
		END IF

	   END IF

	   !  Circular Cressman implies a Euclidean distance without elongation.  This is used
	   !  when the speed of the observation is less than the critical value, for this 
	   !  pressure level.

	   IF ( speed_ob .LT. vcrit ) THEN

		distance = SQRT ( ( ( xob - REAL(crsdot)/2. - REAL(i) ) * ( xob - REAL(crsdot)/2. - REAL(i) ) ) +  &
						  ( ( yob - REAL(crsdot)/2. - REAL(j) ) * ( yob - REAL(crsdot)/2. - REAL(j) ) ) ) 
   
		choice = 'circle  '

	   !  If the wind speed is large enough and the radius of curvature is small enough, then
	   !  the banana shaped ellipse is used.

	   ELSE IF ( ABS(radius_of_curvature) .LT. ( 3.* radius_influence ) * dxd ) THEN

		IF ( choice .EQ. 'undefine' ) THEN

		   !  The radius of curvature needs to have a sign associated with it so that
		   !  from the given point on the radius, we can find the center.  From the 
		   !  observation location, find the four surrounding points that are three
		   !  grid units away.  However, as before, they need to be in the domain.
   
		   ia = MAX ( NINT ( xob ) - 3 ,   1 ) 
		   ja =              yob             
		   ib =              xob              
		   jb = MAX ( NINT ( yob ) - 3 ,   1 ) 
		   ic = MIN ( NINT ( xob ) + 3 , iew ) 
		   jc =              yob             
		   id =              xob              
		   jd = MIN ( NINT ( yob ) + 3 , jns ) 
   
		   !  With these four grid points, get the first-guess u and v at these points.
   
		   ua = u_banana(ia,ja)
		   va = v_banana(ia,ja)
		   ub = u_banana(ib,jb)
		   vb = v_banana(ib,jb)
		   uc = u_banana(ic,jc)
		   vc = v_banana(ic,jc)
		   ud = u_banana(id,jd)
		   vd = v_banana(id,jd)
   
		   !  Get the SIGN of the relative vorticity at the observation point with
		   !  these surrounding values.
      
		   vor = ( vc - va ) - ( ud - ub )
   
		   !  Modify the SIGN of the radius_of_curvature with the SIGN of the relative
		   !  vorticity.
   
		   radius_of_curvature = radius_of_curvature * SIGN ( 1.,vor )
   
		   !  Compute the location of the center of curvature for the streamline
		   !  passing through this observation location.  This is the same as eqn 4.20a
		   !  and 4.20b in MH94, except that i0, ic, j0, jc are not in the reversed
		   !  model directions (i is east west, j is north south).
   
		   i_center = NINT(xob) - radius_of_curvature * v_ob / speed_ob
		   j_center = NINT(yob) + radius_of_curvature * u_ob / speed_ob
      
		END IF