!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_fddaobs_rtfdda 
(docs) 2
! This obs-nudging FDDA module (RTFDDA) is developed by the
! NCAR/RAL/NSAP (National Security Application Programs), under the
! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is
! acknowledged for releasing this capability for WRF community
! research applications.
!
! The NCAR/RAL RTFDDA module was adapted, and significantly modified
! from the obs-nudging module in the standard MM5V3.1 which was originally
! developed by PSU (Stauffer and Seaman, 1994).
!
! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module
! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
! Nov. 2006
!
! References:
!
! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
! implementation of obs-nudging-based FDDA into WRF for supporting
! ATEC test operations. 2005 WRF user workshop. Paper 10.7.
!
! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update
! on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE
! and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7.
!
! Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data
! assimilation. J. Appl. Meteor., 33, 416-434.
!
! http://www.rap.ucar.edu/projects/armyrange/references.html
!
! Modification History:
! 03212007 Modified fddaobs_init to compute Lambert cone factor. -Al Bourgeois
CONTAINS
!------------------------------------------------------------------------------
SUBROUTINE fddaobs_init (docs) (obs_nudge_opt, maxdom, inest, parid, & 1,16
idynin, dtramp, fdaend, restart, &
obs_twindo_cg, obs_twindo, itimestep, &
xlat, xlong, &
#if ( EM_CORE == 1 )
fdob, &
#endif
iprt, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte)
!-----------------------------------------------------------------------
! This routine does initialization for real time fdda obs-nudging.
!
!-----------------------------------------------------------------------
USE module_domain
USE module_dm ! for wrf_dm_min_real
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!=======================================================================
! Definitions
!-----------
INTEGER, intent(in) :: maxdom
INTEGER, intent(in) :: obs_nudge_opt(maxdom)
INTEGER, intent(in) :: ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte
INTEGER, intent(in) :: inest
INTEGER, intent(in) :: parid(maxdom)
INTEGER, intent(in) :: idynin ! flag for dynamic initialization
REAL, intent(in) :: dtramp ! time period for idynin ramping (min)
REAL, intent(in) :: fdaend(maxdom) ! nudging end time for domain (min)
LOGICAL, intent(in) :: restart
REAL, intent(in) :: obs_twindo_cg ! time window on coarse grid
REAL, intent(in) :: obs_twindo
INTEGER, intent(in) :: itimestep
REAL, DIMENSION( ims:ime, jms:jme ), &
INTENT(IN) :: xlat, xlong ! lat/long locations on mass point grid
#if ( EM_CORE == 1 )
TYPE(fdob_type), intent(inout) :: fdob
#endif
LOGICAL, intent(in) :: iprt ! Flag enabling printing warning messages
! Local variables
logical :: nudge_flag ! nudging flag for this nest
integer :: ktau ! current timestep
integer :: nest ! loop counter
integer :: idom ! domain id
integer :: parent ! parent domain
real :: conv ! 180/pi
real :: tl1 ! Lambert standard parallel 1
real :: tl2 ! Lambert standard parallel 2
real :: xn1
real :: known_lat ! Latitude of domain point (i,j)=(1,1)
real :: known_lon ! Longitude of domain point (i,j)=(1,1)
character(len=200) :: msg ! Argument to wrf_message
#if ( EM_CORE == 1 )
! Set flag for nudging on pressure (not sigma) surfaces.
fdob%iwtsig = 0
!**************************************************************************
! *** Initialize datend for dynamic initialization (ajb added 08132008) ***
!**************************************************************************
! Set ending nudging date (used with dynamic ramp-down) to zero.
fdob%datend = 0.
fdob%ieodi = 0
! Check for dynamic initialization flag
if(idynin.eq.1)then
! Set datend to time in minutes after which data are assumed to have ended.
if(dtramp.gt.0.)then
fdob%datend = fdaend(inest) - dtramp
else
fdob%datend = fdaend(inest)
endif
if(iprt) then
call wrf_message("")
write(msg,'(a,i3,a)') &
' *** DYNAMIC-INITIALIZATION OPTION FOR INEST = ', inest, ' ***'
call wrf_message(msg)
write(msg,*) ' FDAEND,DATEND,DTRAMP: ',fdaend(inest),fdob%datend,dtramp
call wrf_message(msg)
call wrf_message("")
endif
endif
! *** end dynamic initialization section ***
!**************************************************************************
! Set time window.
fdob%window = obs_twindo
if(inest.eq.1) then
if(obs_twindo .eq. 0.) then
if(iprt) then
call wrf_message("")
write(msg,*) '*** WARNING: TWINDO=0 on the coarse domain.'
call wrf_message(msg)
write(msg,*) '*** Did you forget to set twindo in the fdda namelist?'
call wrf_message(msg)
call wrf_message("")
endif
endif
else ! nest
if(obs_twindo .eq. 0.) then
fdob%window = obs_twindo_cg
if(iprt) then
call wrf_message("")
write(msg,'(a,i2)') 'WARNING: TWINDO=0. for nest ',inest
call wrf_message(msg)
write(msg,'(a,f12.5,a)') 'Default to coarse-grid value of ', obs_twindo_cg,' hrs'
call wrf_message(msg)
call wrf_message("")
endif
endif
endif
! Initialize flags.
fdob%domain_tot=0
do nest=1,maxdom
fdob%domain_tot = fdob%domain_tot + obs_nudge_opt(nest)
end do
! Set parameters.
fdob%pfree = 50.0
fdob%rinfmn = 1.0
fdob%rinfmx = 2.0
fdob%dpsmx = 7.5
fdob%dcon = 1.0/fdob%dpsmx
! Get known lat and lon and store these on all processors in fdob structure, for
! later use in projection x-formation to map (lat,lon) to (x,y) for each obs.
IF (its .eq. 1 .AND. jts .eq. 1) then
known_lat = xlat(1,1)
known_lon = xlong(1,1)
ELSE
known_lat = 9999.
known_lon = 9999.
END IF
fdob%known_lat = wrf_dm_min_real(known_lat)
fdob%known_lon = wrf_dm_min_real(known_lon)
! Calculate the nest levels, levidn. Note that each nest
! must know the nest levels levidn(maxdom) of each domain.
do nest=1,maxdom
! Initialize nest level for each domain.
if (nest .eq. 1) then
fdob%levidn(nest) = 0 ! Mother domain has nest level 0
else
fdob%levidn(nest) = 1 ! All other domains start with 1
endif
idom = nest
100 parent = parid(idom) ! Go up the parent tree
if (parent .gt. 1) then ! If not up to mother domain
fdob%levidn(nest) = fdob%levidn(nest) + 1
idom = parent
goto 100
endif
enddo
! Check to see if the nudging flag has been set. If not,
! simply RETURN.
nudge_flag = (obs_nudge_opt(inest) .eq. 1)
if (.not. nudge_flag) return
ktau = itimestep
if(restart) then
fdob%ktaur = ktau
else
fdob%ktaur = 0
endif
RETURN
#endif
END SUBROUTINE fddaobs_init
#if ( EM_CORE == 1 )
!-----------------------------------------------------------------------
SUBROUTINE errob (docs) (inest, ub, vb, tb, t0, qvb, pbase, pp, rovcp, & 1,13
phb, ph, g, &
uratx, vratx, tratx, &
nndgv, nerrf, niobf, maxdom, &
levidn, parid, nstat, nstaw, &
iswind, istemp, ismois, ispstr, &
timeob, rio, rjo, rko, &
varobs, errf, ktau, xtime, &
iratio, npfi, &
prt_max, prt_freq, iprt, &
obs_prt, lat_prt, lon_prt, &
mlat_prt, mlon_prt, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
!-----------------------------------------------------------------------
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
USE module_dm, ONLY : get_full_obs_vector
#endif
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
! PURPOSE: THIS SUBROUTINE CALCULATES THE DIFFERENCE BETWEEN THE
! OBSERVED VALUES AND THE FORECAST VALUES AT THE OBSERVATION
! POINTS. THE OBSERVED VALUES CLOSEST TO THE CURRENT
! FORECAST TIME (XTIME) WERE DETERMINED IN SUBROUTINE
! IN4DOB AND STORED IN ARRAY VAROBS. THE DIFFERENCES
! CALCULATED BY SUBROUTINE ERROB WILL BE STORED IN ARRAY
! ERRF FOR THE NSTA OBSERVATION LOCATIONS. MISSING
! OBSERVATIONS ARE DENOTED BY THE DUMMY VALUE 99999.
!
! HISTORY: Original author: MM5 version???
! 02/04/2004 - Creation of WRF version. Al Bourgeois
! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
!------------------------------------------------------------------------------
! THE STORAGE ORDER IN ERRF IS AS FOLLOWS:
! IVAR VARIABLE TYPE(TAU-1)
! ---- --------------------
! 1 U error
! 2 V error
! 3 T error
! 4 Q error
! 5 Surface press error at T points (not used)
! 6 Model surface press at T-points
! 7 Model surface press at U-points
! 8 Model surface press at V-points
! 9 RKO at U-points
!-----------------------------------------------------------------------
!
! Description of input arguments.
!
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: inest ! Domain index.
INTEGER, INTENT(IN) :: nndgv ! Number of nudge variables.
INTEGER, INTENT(IN) :: nerrf ! Number of error fields.
INTEGER, INTENT(IN) :: niobf ! Number of observations.
INTEGER, INTENT(IN) :: maxdom ! Maximum number of domains.
INTEGER, INTENT(IN) :: levidn(maxdom) ! Level of nest.
INTEGER, INTENT(IN) :: parid(maxdom) ! Id of parent grid.
INTEGER, INTENT(IN) :: ktau ! Model time step index
REAL, INTENT(IN) :: xtime ! Forecast time in minutes
INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
INTEGER, INTENT(IN) :: npfi ! Coarse-grid diagnostics freq.
INTEGER, INTENT(IN) :: prt_max ! Max number of obs to print.
INTEGER, INTENT(IN) :: prt_freq ! Frequency of obs to print.
LOGICAL, INTENT(IN) :: iprt ! Print flag
INTEGER, INTENT(INOUT) :: obs_prt(prt_max) ! Print obs indices
REAL, INTENT(INOUT) :: lat_prt(prt_max) ! Print obs latitude
REAL, INTENT(INOUT) :: lon_prt(prt_max) ! Print obs longitude
REAL, INTENT(INOUT) :: mlat_prt(prt_max) ! Print model lat at obs loc
REAL, INTENT(INOUT) :: mlon_prt(prt_max) ! Print model lon at obs loc
INTEGER, INTENT(IN) :: nstat ! # stations held for use
INTEGER, INTENT(IN) :: nstaw ! # stations in current window
INTEGER, intent(in) :: iswind
INTEGER, intent(in) :: istemp
INTEGER, intent(in) :: ismois
INTEGER, intent(in) :: ispstr
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
REAL, INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
REAL, INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
REAL, INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
REAL, INTENT(IN) :: t0
REAL, INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
REAL, INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme )
REAL, INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme ) ! Press. perturbation (Pa)
REAL, INTENT(IN) :: rovcp
REAL, INTENT(IN) :: phb( ims:ime, kms:kme, jms:jme ) ! geopotential (base)
REAL, INTENT(IN) :: ph( ims:ime, kms:kme, jms:jme ) ! geopotential (perturbation)
REAL, INTENT(IN) :: g ! gravity constant
REAL, INTENT(IN) :: uratx( ims:ime, jms:jme ) ! U to U10 ratio on mass points.
REAL, INTENT(IN) :: vratx( ims:ime, jms:jme ) ! V to V10 ratio on mass points.
REAL, INTENT(IN) :: tratx( ims:ime, jms:jme ) ! T to TH2 ratio on mass points.
REAL, INTENT(IN) :: timeob(niobf) ! Obs time (hrs)
REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
REAL, INTENT(INOUT) :: rko(niobf) ! Obs bottom-top coordinate
REAL, INTENT(INOUT) :: varobs(nndgv, niobf)
REAL, INTENT(INOUT) :: errf(nerrf, niobf)
! Local variables
INTEGER :: iobmg(niobf) ! Obs i-coord on mass grid
INTEGER :: jobmg(niobf) ! Obs j-coord on mass grid
INTEGER :: ia(niobf)
INTEGER :: ib(niobf)
INTEGER :: ic(niobf)
REAL :: pbbo(kds:kde) ! column base pressure (cb) at obs loc.
REAL :: ppbo(kds:kde) ! column pressure perturbation (cb) at obs loc.
REAL :: ra(niobf)
REAL :: rb(niobf)
REAL :: rc(niobf)
REAL :: dxobmg(niobf) ! Interp. fraction (x dir) referenced to mass-grid
REAL :: dyobmg(niobf) ! Interp. fraction (y dir) referenced to mass-grid
INTEGER MM(MAXDOM)
INTEGER NNL
real :: uratio( ims:ime, jms:jme ) ! U to U10 ratio on momentum points.
real :: vratio( ims:ime, jms:jme ) ! V to V10 ratio on momentum points.
real :: pug1,pug2,pvg1,pvg2
character(len=200) :: msg ! Argument to wrf_message
! Define staggers for U, V, and T grids, referenced from non-staggered grid.
real, parameter :: gridx_t = 0.5 ! Mass-point x stagger
real, parameter :: gridy_t = 0.5 ! Mass-point y stagger
real, parameter :: gridx_u = 0.0 ! U-point x stagger
real, parameter :: gridy_u = 0.5 ! U-point y stagger
real, parameter :: gridx_v = 0.5 ! V-point x stagger
real, parameter :: gridy_v = 0.0 ! V-point y stagger
real :: dummy = 99999.
real :: pbhi, pphi
real :: press,ttemp !ajb scratch variables
! real model_temp,pot_temp !ajb scratch variables
!*** DECLARATIONS FOR IMPLICIT NONE
integer nsta,ivar,n,ityp
integer iob,job,kob,iob_ms,job_ms
integer k,kbot,nml,nlb,nle
integer iobm,jobm,iobp,jobp,kobp,inpf,i,j
integer i_start,i_end,j_start,j_end ! loop ranges for uratio,vratio calc.
integer k_start,k_end
integer ips ! For printing obs information
integer pnx ! obs index for printout
real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
real pob
real hob
real grfacx,grfacy,uratiob,vratiob,tratiob,tratxob,fnpf
real stagx ! For x correction to mass-point stagger
real stagy ! For y correction to mass-point stagger
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
LOGICAL MP_LOCAL_DUMMASK(NIOBF) ! Mask for work to be done on this processor
LOGICAL MP_LOCAL_UOBMASK(NIOBF) ! Dot-point mask
LOGICAL MP_LOCAL_VOBMASK(NIOBF) ! Dot-point mask
LOGICAL MP_LOCAL_COBMASK(NIOBF) ! Cross-point mask
#endif
! LOGICAL, EXTERNAL :: TILE_MASK
NSTA=NSTAT
! FIRST, DETERMINE THE GRID TYPE CORRECTION FOR U-momentum, V-momentum,
! AND MASS POINTS, AND WHEN INEST=2, CONVERT THE STORED COARSE MESH INDICES
! TO THE FINE MESH INDEX EQUIVALENTS
! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
if (iprt) then
write(msg,'(a,i5,a,i2,a,i5,a)') '++++++CALL ERROB AT KTAU = ', &
KTAU,' AND INEST = ',INEST,': NSTA = ',NSTAW,' ++++++'
call wrf_message(msg)
endif
ERRF = 0.0 ! Zero out errf array
! Set up loop bounds for this grid's boundary conditions
i_start = max( its-1,ids )
i_end = min( ite+1,ide-1 )
j_start = max( jts-1,jds )
j_end = min( jte+1,jde-1 )
k_start = kts
k_end = min( kte, kde-1 )
DO ITYP=1,3 ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP
! Set grid stagger
IF(ITYP.EQ.1) THEN ! U-POINT CASE
GRIDX = gridx_u
GRIDY = gridy_u
ELSE IF(ITYP.EQ.2) THEN ! V-POINT CASE
GRIDX = gridx_v
GRIDY = gridy_v
ELSE ! MASS-POINT CASE
GRIDX = gridx_t
GRIDY = gridy_t
ENDIF
! Compute URATIO and VRATIO fields on momentum (u,v) points.
IF(ityp.eq.1)THEN
call upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, uratx, uratio)
ELSE IF (ityp.eq.2) THEN
call vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, vratx, vratio)
ENDIF
IF(INEST.EQ.1) THEN ! COARSE MESH CASE...
DO N=1,NSTA
RA(N)=RIO(N)-GRIDX
RB(N)=RJO(N)-GRIDY
IA(N)=RA(N)
IB(N)=RB(N)
IOB=MAX0(1,IA(N))
IOB=MIN0(IOB,ide-1)
JOB=MAX0(1,IB(N))
JOB=MIN0(JOB,jde-1)
DXOB=RA(N)-FLOAT(IA(N))
DYOB=RB(N)-FLOAT(IB(N))
! Save mass-point arrays for computing rko for all var types
if(ityp.eq.1) then
iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
endif
iob_ms = iobmg(n)
job_ms = jobmg(n)
dxob_ms = dxobmg(n)
dyob_ms = dyobmg(n)
!if(n.eq.1 .and. iprt) then
! write(6,*) 'ERROB - COARSE MESH:'
! write(6,'(a,i1,a,i1,4(a,f5.2),2(a,i3),2(a,f6.3))') 'OBS= ',n, &
! ' ityp= ',ityp, &
! ' ra= ',ra(n),' rb= ',rb(n), &
! ' rio= ',rio(n),' rjo= ',rjo(n), &
! ' iob= ',iob,' job= ',job, &
! ' dxob= ',dxob,' dyob= ',dyob
! write(6,'(a,i3,a,i3,a,f5.2,a,f5.2)') &
! ' iob_ms= ',iob_ms,' job_ms= ',job_ms, &
! ' dxob_ms= ',dxob_ms,' dyob_ms= ',dyob_ms
!endif
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! Set mask for obs to be handled by this processor
MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
IF ( MP_LOCAL_DUMMASK(N) ) THEN
#endif
! Interpolate pressure to obs location column and convert from Pa to cb.
do k = kds, kde
pbbo(k) = .001*( &
(1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) + &
DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
ppbo(k) = .001*( &
(1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) + &
DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
! write(6,'(a,i2,2(a,f9.3)') ' k= ',k,' pbbo= ',pbbo(k),' ppbo= ',ppbo(k)
enddo
!ajb 20040119: Note based on bugfix for dot/cross points split across processors,
!ajb which was actually a serial code fix: The ityp=2 (v-points) and
!ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
!ajb case rko! This is necessary for bit-for-bit reproducability
!ajb with the parallel run. (coarse mesh)
if(abs(rko(n)+99).lt.1.)then
pob = varobs(5,n)
if(pob .eq.-888888.) then
hob = varobs(6,n)
if(hob .gt. -800000. ) then
pob = ht_to_p( hob, ppbo, pbbo, phb, ph, g, iob_ms, job_ms, &
dxob_ms, dyob_ms, k_start, k_end, kds,kde, &
ims,ime, jms,jme, kms,kme )
endif
endif
if(pob .gt.-800000.)then
do k=k_end-1,1,-1
kbot = k
if(pob .le. pbbo(k)+ppbo(k)) then
goto 199
endif
enddo
199 continue
pphi = ppbo(kbot+1)
pbhi = pbbo(kbot+1)
rko(n) = real(kbot+1)- &
( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
rko(n)=max(rko(n),1.0)
endif
endif
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
#endif
RC(N)=RKO(N)
ENDDO ! END COARSE MESH LOOP OVER NSTA
ELSE ! FINE MESH CASE
!**********************************************************************
!ajb_07012008: Conversion of obs coordinates to the fine mesh here
!ajb is no longer necessary, since implementation of the WRF map pro-
!ajb jections (to each nest directly) is implemented in sub in4dob.
!**********************************************************************
!ajb
!ajb GET (I,J,K) OF OBSERVATIONS ON FINE MESH VALUES.
DO N=1,NSTA
! write(6,*) 'UPDATE: ra(n) = ',ra(n),' rb(n) = ',rb(n)
RA(N)=RIO(N)-GRIDX ! ajb added 07012008
RB(N)=RJO(N)-GRIDY ! ajb added 07012008
IA(N)=RA(N)
IB(N)=RB(N)
IOB=MAX0(1,IA(N))
IOB=MIN0(IOB,ide-1)
JOB=MAX0(1,IB(N))
JOB=MIN0(JOB,jde-1)
DXOB=RA(N)-FLOAT(IA(N))
DYOB=RB(N)-FLOAT(IB(N))
! Save mass-point arrays for computing rko for all var types
if(ityp.eq.1) then
stagx = grfacx - gridx_t !Correct x stagger to mass-point
stagy = grfacy - gridy_t !Correct y stagger to mass-point
iobmg(n) = MIN0(MAX0(1,int(RA(n)+stagx)),ide-1)
jobmg(n) = MIN0(MAX0(1,int(RB(n)+stagy)),jde-1)
dxobmg(n) = RA(N)+stagx-FLOAT(int(RA(N)+stagx))
dyobmg(n) = RB(N)+stagy-FLOAT(int(RB(N)+stagy))
endif
iob_ms = iobmg(n)
job_ms = jobmg(n)
dxob_ms = dxobmg(n)
dyob_ms = dyobmg(n)
!if(n.eq.1) then
! write(6,*) 'ERROB - FINE MESH:'
! write(6,*) 'RA = ',ra(n),' RB = ',rb(n)
! write(6,'(a,i1,a,i1,4(a,f5.2),2(a,i3),2(a,f6.3))') 'OBS= ',n, &
! ' ityp= ',ityp, &
! ' ra= ',ra(n),' rb= ',rb(n), &
! ' rio= ',rio(n),' rjo= ',rjo(n), &
! ' iob= ',iob,' job= ',job, &
! ' dxob= ',dxob,' dyob= ',dyob
! write(6,'(a,i3,a,i3,a,f5.2,a,f5.2)') &
! ' iob_ms= ',iob_ms,' job_ms= ',job_ms, &
! ' dxob_ms= ',dxob_ms,' dyob_ms= ',dyob_ms
!endif
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! Set mask for obs to be handled by this processor
MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
IF ( MP_LOCAL_DUMMASK(N) ) THEN
#endif
! Interpolate pressure to obs location column and convert from Pa to cb.
do k = kds, kde
pbbo(k) = .001*( &
(1.-DYOB_MS)*( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS) + &
DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
ppbo(k) = .001*( &
(1.-DYOB_MS)*( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS) + &
DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
! write(6,'(a,i2,2(a,f9.3)') ' k= ',k,' pbbo= ',pbbo(k),' ppbo= ',ppbo(k)
enddo
!ajb 20040119: Note based on bugfix for dot/cross points split across processors,
!ajb which was actually a serial code fix: The ityp=2 (v-points) and
!ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
!ajb case) rko! This is necessary for bit-for-bit reproducability
!ajb with parallel run. (fine mesh)
if(abs(rko(n)+99).lt.1.)then
pob = varobs(5,n)
if(pob .eq.-888888.) then
hob = varobs(6,n)
if(hob .gt. -800000. ) then
pob = ht_to_p( hob, ppbo, pbbo, phb, ph, g, iob_ms, job_ms, &
dxob_ms, dyob_ms, k_start, k_end, kds,kde, &
ims,ime, jms,jme, kms,kme )
endif
endif
if(pob .gt.-800000.)then
do k=k_end-1,1,-1
kbot = k
if(pob .le. pbbo(k)+ppbo(k)) then
goto 198
endif
enddo
198 continue
pphi = ppbo(kbot+1)
pbhi = pbbo(kbot+1)
rko(n) = real(kbot+1)- &
( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
rko(n)=max(rko(n),1.0)
endif
endif
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
#endif
RC(N)=RKO(N)
ENDDO ! END FINE MESH LOOP OVER NSTA
ENDIF ! end if(inest.eq.1)
! Print obs information.
if(ityp.eq.3) then !mass-point case
CALL print_obs_info(iprt,inest,niobf,rio,rjo,rko, &
prt_max,prt_freq,obs_prt,lat_prt,lon_prt, &
mlat_prt,mlon_prt,timeob,xtime)
endif
! INITIALIZE THE ARRAY OF DIFFERENCES BETWEEN THE OBSERVATIONS
! AND THE LOCAL FORECAST VALUES TO ZERO. FOR THE FINE MESH
! ONLY, SET THE DIFFERENCE TO A LARGE DUMMY VALUE IF THE
! OBSERVATION IS OUTSIDE THE FINE MESH DOMAIN.
! SET DIFFERENCE VALUE TO A DUMMY VALUE FOR OBS POINTS OUTSIDE
! CURRENT DOMAIN
IF(ITYP.EQ.1) THEN
NLB=1
NLE=1
ELSE IF(ITYP.EQ.2) THEN
NLB=2
NLE=2
ELSE
NLB=3
NLE=5
ENDIF
DO IVAR=NLB,NLE
DO N=1,NSTA
IF((RA(N)-1.).LT.0)THEN
ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
ENDIF
IF((RB(N)-1.).LT.0)THEN
ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
ENDIF
IF((FLOAT(ide)-2.0*GRIDX-RA(N)-1.E-10).LT.0)THEN
ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
ENDIF
IF((FLOAT(jde)-2.0*GRIDY-RB(N)-1.E-10).LT.0)THEN
ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
ENDIF
if(rc(n).lt.1.)errf(ivar,n)=errf(ivar,n)+dummy
ENDDO
ENDDO
! NOW FIND THE EXACT OFFSET OF EACH OBSERVATION FROM THE
! GRID POINT TOWARD THE LOWER LEFT
DO N=1,NSTA
IA(N)=RA(N)
IB(N)=RB(N)
IC(N)=RC(N)
ENDDO
DO N=1,NSTA
RA(N)=RA(N)-FLOAT(IA(N))
RB(N)=RB(N)-FLOAT(IB(N))
RC(N)=RC(N)-FLOAT(IC(N))
ENDDO
! PERFORM A TRILINEAR EIGHT-POINT (3-D) INTERPOLATION
! TO FIND THE FORECAST VALUE AT THE EXACT OBSERVATION
! POINTS FOR U, V, T, AND Q.
! Compute local masks for dot and cross points.
if(ityp.eq.1) then
DO N=1,NSTA
IOB=MAX0(1,IA(N))
IOB=MIN0(IOB,ide-1)
JOB=MAX0(1,IB(N))
JOB=MIN0(JOB,jde-1)
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! Set mask for U-momemtum points to be handled by this processor
MP_LOCAL_UOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
#endif
ENDDO
endif
if(ityp.eq.2) then
DO N=1,NSTA
IOB=MAX0(1,IA(N))
IOB=MIN0(IOB,ide-1)
JOB=MAX0(1,IB(N))
JOB=MIN0(JOB,jde-1)
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! Set mask for V-momentum points to be handled by this processor
MP_LOCAL_VOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
#endif
ENDDO
endif
if(ityp.eq.3) then
DO N=1,NSTA
IOB=MAX0(1,IA(N))
IOB=MIN0(IOB,ide-1)
JOB=MAX0(1,IB(N))
JOB=MIN0(JOB,jde-1)
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! Set mask for cross (mass) points to be handled by this processor
MP_LOCAL_COBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
#endif
ENDDO
endif
!**********************************************************
! PROCESS U VARIABLE (IVAR=1)
!**********************************************************
IF(ITYP.EQ.1) THEN
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
DO N=1,NSTA
IF(MP_LOCAL_UOBMASK(N)) THEN
ERRF(9,N)=rko(n) !RKO is needed by neighboring processors !ajb
ENDIF
ENDDO
#endif
IF(ISWIND.EQ.1) THEN
DO N=1,NSTA
IOB=MAX0(2,IA(N))
IOB=MIN0(IOB,ide-1)
IOBM=MAX0(1,IOB-1)
IOBP=MIN0(ide-1,IOB+1)
JOB=MAX0(2,IB(N))
JOB=MIN0(JOB,jde-1)
JOBM=MAX0(1,JOB-1)
JOBP=MIN0(jde-1,JOB+1)
KOB=MIN0(K_END,IC(N))
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
IF(MP_LOCAL_UOBMASK(N))THEN ! Do if obs on this processor
#endif
KOBP=MIN0(KOB+1,k_end)
DXOB=RA(N)
DYOB=RB(N)
DZOB=RC(N)
! Compute surface pressure values at surrounding U and V points
PUG1 = .5*( pbase(IOBM,1,JOB) + pbase(IOB,1,JOB) )
PUG2 = .5*( pbase(IOB,1,JOB) + pbase(IOBP,1,JOB) )
! This is to correct obs to model sigma level using reverse similarity theory
if(rko(n).eq.1.0)then
uratiob=((1.-DYOB)*((1.-DXOB)*uratio(IOB,JOB)+ &
DXOB*uratio(IOBP,JOB) &
)+DYOB*((1.-DXOB)*uratio(IOB,JOBP)+ &
DXOB*uratio(IOBP,JOBP)))
else
uratiob=1.
endif
!YLIU Some PBL scheme do not define the vratio/uratio
if(abs(uratiob).lt.1.0e-3) then
uratiob=1.
endif
! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
! OUTSIDE THE CURRENT DOMAIN
! U COMPONENT WIND ERROR
ERRF(1,N)=ERRF(1,N)+uratiob*VAROBS(1,N)-((1.-DZOB)* &
((1.-DyOB)*((1.- &
DxOB)*UB(IOB,KOB,JOB)+DxOB*UB(IOB+1,KOB,JOB) &
)+DyOB*((1.-DxOB)*UB(IOB,KOB,JOB+1)+DxOB* &
UB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
*UB(IOB,KOBP,JOB)+DxOB*UB(IOB+1,KOBP,JOB))+ &
DyOB*((1.-DxOB)*UB(IOB,KOBP,JOB+1)+DxOB* &
UB(IOB+1,KOBP,JOB+1))))
! if(n.le.10) then
! write(6,*)
! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF1 at ',iob,job,kob, &
! ' N = ',n,' inest = ',inest
! write(6,*) 'VAROBS(1,N) = ',varobs(1,n)
! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
! write(6,*) 'UB(IOB,KOB,JOB) = ',UB(IOB,KOB,JOB)
! write(6,*) 'UB(IOB+1,KOB,JOB) = ',UB(IOB+1,KOB,JOB)
! write(6,*) 'UB(IOB,KOB,JOB+1) = ',UB(IOB,KOB,JOB+1)
! write(6,*) 'UB(IOB+1,KOB,JOB+1) = ',UB(IOB+1,KOB,JOB+1)
! write(6,*) 'UB(IOB,KOBP,JOB) = ',UB(IOB,KOBP,JOB)
! write(6,*) 'UB(IOB+1,KOBP,JOB) = ',UB(IOB+1,KOBP,JOB)
! write(6,*) 'UB(IOB,KOBP,JOB+1) = ',UB(IOB,KOBP,JOB+1)
! write(6,*) 'UB(IOB+1,KOBP,JOB+1) = ',UB(IOB+1,KOBP,JOB+1)
! write(6,*) 'uratiob = ',uratiob
! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
! write(6,*) 'ERRF(1,N) = ',errf(1,n)
! write(6,*)
! endif
! Store model surface pressure (not the error!) at U point.
ERRF(7,N)=.001*( (1.-DXOB)*PUG1 + DXOB*PUG2 )
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
ENDIF ! end IF( MP_LOCAL_UOBMASK(N) )
#endif
ENDDO ! END U-POINT LOOP OVER OBS
ENDIF ! end if(iswind.eq.1)
ENDIF ! ITYP=1: PROCESS U
!**********************************************************
! PROCESS V VARIABLE (IVAR=2)
!**********************************************************
IF(ITYP.EQ.2) THEN
IF(ISWIND.EQ.1) THEN
DO N=1,NSTA
IOB=MAX0(2,IA(N))
IOB=MIN0(IOB,ide-1)
IOBM=MAX0(1,IOB-1)
IOBP=MIN0(ide-1,IOB+1)
JOB=MAX0(2,IB(N))
JOB=MIN0(JOB,jde-1)
JOBM=MAX0(1,JOB-1)
JOBP=MIN0(jde-1,JOB+1)
KOB=MIN0(K_END,IC(N))
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
IF(MP_LOCAL_VOBMASK(N))THEN ! Do if obs on this processor
#endif
KOBP=MIN0(KOB+1,k_end)
DXOB=RA(N)
DYOB=RB(N)
DZOB=RC(N)
! Compute surface pressure values at surrounding U and V points
PVG1 = .5*( pbase(IOB,1,JOBM) + pbase(IOB,1,JOB) )
PVG2 = .5*( pbase(IOB,1,JOB) + pbase(IOB,1,JOBP) )
! This is to correct obs to model sigma level using reverse similarity theory
if(rko(n).eq.1.0)then
vratiob=((1.-DYOB)*((1.-DXOB)*vratio(IOB,JOB)+ &
DXOB*vratio(IOBP,JOB) &
)+DYOB*((1.-DXOB)*vratio(IOB,JOBP)+ &
DXOB*vratio(IOBP,JOBP)))
else
vratiob=1.
endif
!YLIU Some PBL scheme do not define the vratio/uratio
if(abs(vratiob).lt.1.0e-3) then
vratiob=1.
endif
! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
! OUTSIDE THE CURRENT DOMAIN
! V COMPONENT WIND ERROR
ERRF(2,N)=ERRF(2,N)+vratiob*VAROBS(2,N)-((1.-DZOB)* &
((1.-DyOB)*((1.- &
DxOB)*VB(IOB,KOB,JOB)+DxOB*VB(IOB+1,KOB,JOB) &
)+DyOB*((1.-DxOB)*VB(IOB,KOB,JOB+1)+DxOB* &
VB(IOB+1,KOB,JOB+1)))+DZOB*((1.-DyOB)*((1.-DxOB) &
*VB(IOB,KOBP,JOB)+DxOB*VB(IOB+1,KOBP,JOB))+ &
DyOB*((1.-DxOB)*VB(IOB,KOBP,JOB+1)+DxOB* &
VB(IOB+1,KOBP,JOB+1))))
! if(n.le.10) then
! write(6,*)
! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF2 at ',iob,job,kob, &
! ' N = ',n,' inest = ',inest
! write(6,*) 'VAROBS(2,N) = ',varobs(2,n)
! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
! write(6,*) 'VB(IOB,KOB,JOB) = ',VB(IOB,KOB,JOB)
! write(6,*) 'ERRF(2,N) = ',errf(2,n)
! write(6,*) 'vratiob = ',vratiob
! write(6,*)
! endif
! Store model surface pressure (not the error!) at V point.
ERRF(8,N)=.001*( (1.-DYOB)*PVG1 + DYOB*PVG2 )
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
ENDIF ! end IF( MP_LOCAL_VOBMASK(N) )
#endif
ENDDO ! END V-POINT LOOP OVER OBS
ENDIF ! end if(iswind.eq.1)
ENDIF ! ITYP=1: PROCESS V
!**********************************************************
! PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
!**********************************************************
IF(ITYP.EQ.3) THEN
IF(ISTEMP.EQ.1 .OR. ISMOIS.EQ.1) THEN
DO N=1,NSTA
IOB=MAX0(1,IA(N))
IOB=MIN0(IOB,ide-1)
JOB=MAX0(1,IB(N))
JOB=MIN0(JOB,jde-1)
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
#endif
KOB=MIN0(k_end,IC(N))
KOBP=MIN0(KOB+1,K_END)
DXOB=RA(N)
DYOB=RB(N)
DZOB=RC(N)
! This is to correct obs to model sigma level using reverse similarity theory
if(rko(n).eq.1.0)then
tratxob=((1.-DYOB)*((1.-DXOB)*tratx(IOB,JOB)+ &
DXOB*tratx(IOB+1,JOB) &
)+DYOB*((1.-DXOB)*tratx(IOB,JOB+1)+ &
DXOB*tratx(IOB+1,JOB+1)))
else
tratxob=1.
endif
!yliu
if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
!ajb testing only
! if(iprt .and. n.eq.81) then
! write(6,*) 'POTENTIAL TEMP FOR N=81:'
! write(6,*)
! write(6,*) ' K THETA TEMPERATURE', &
! ' PBASE'
! write(6,*)
! do k=k_end,1,-1
! press = pbase(iob,k,job)+pp(iob,k,job)
! ttemp = (t0+tb(iob,k,job)) / (100000./press)**.2857143
! write(6,*) k,t0+TB(IOB,k,JOB),ttemp,pbase(iob,k,job)
! enddo
! endif
!ajb end testing only
! TEMPERATURE ERROR
! if(n.le.10) then
! write(6,*) 'before: errf(3,n) = ',errf(3,n)
! endif
ERRF(3,N)=ERRF(3,N)+tratxob*VAROBS(3,N)-((1.-DZOB)* &
((1.-DyOB)*((1.- &
DxOB)*(TB(IOB,KOB,JOB))+DxOB* &
(TB(IOB+1,KOB,JOB)))+DyOB*((1.-DxOB)* &
(TB(IOB,KOB,JOB+1))+DxOB* &
(TB(IOB+1,KOB,JOB+1))))+DZOB*((1.- &
DyOB)*((1.-DxOB)*(TB(IOB,KOBP,JOB))+DxOB* &
(TB(IOB+1,KOBP,JOB)))+DyOB*((1.-DxOB)* &
(TB(IOB,KOBP,JOB+1))+DxOB* &
(TB(IOB+1,KOBP,JOB+1)))))
! if(n.le.10) then
! write(6,*)
! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF3 at ',iob,job,kob, &
! ' N = ',n,' inest = ',inest
! write(6,*) 'VAROBS(3,N) = ',varobs(3,n)
! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
! write(6,*) 'TB(IOB,KOB,JOB) = ',TB(IOB,KOB,JOB)
! write(6,*) 'TB(IOB+1,KOB,JOB) = ',TB(IOB+1,KOB,JOB)
! write(6,*) 'TB(IOB,KOB,JOB+1) = ',TB(IOB,KOB,JOB+1)
! write(6,*) 'TB(IOB+1,KOB,JOB+1) = ',TB(IOB+1,KOB,JOB+1)
! write(6,*) 'TB(IOB,KOBP,JOB) = ',TB(IOB,KOBP,JOB)
! write(6,*) 'TB(IOB+1,KOBP,JOB) = ',TB(IOB+1,KOBP,JOB)
! write(6,*) 'TB(IOB,KOBP,JOB+1) = ',TB(IOB,KOBP,JOB+1)
! write(6,*) 'TB(IOB+1,KOBP,JOB+1) = ',TB(IOB+1,KOBP,JOB+1)
! write(6,*) 'tratxob = ',tratxob
! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
! write(6,*) 'ERRF(3,N) = ',errf(3,n)
! write(6,*)
! endif
! MOISTURE ERROR
ERRF(4,N)=ERRF(4,N)+VAROBS(4,N)-((1.-DZOB)*((1.-DyOB)*((1.- &
DxOB)*QVB(IOB,KOB,JOB)+DxOB* &
QVB(IOB+1,KOB,JOB))+DyOB*((1.-DxOB)* &
QVB(IOB,KOB,JOB+1)+DxOB* &
QVB(IOB+1,KOB,JOB+1)))+DZOB*((1.- &
DyOB)*((1.-DxOB)*QVB(IOB,KOBP,JOB)+DxOB &
*QVB(IOB+1,KOBP,JOB))+DyOB*((1.-DxOB &
)*QVB(IOB,KOBP,JOB+1)+DxOB* &
QVB(IOB+1,KOBP,JOB+1))))
! Store model surface pressure (not the error!) at T-point
ERRF(6,N)= .001* &
((1.-DyOB)*((1.-DxOB)*pbase(IOB,1,JOB)+DxOB* &
pbase(IOB+1,1,JOB))+DyOB*((1.-DxOB)* &
pbase(IOB,1,JOB+1)+DxOB*pbase(IOB+1,1,JOB+1) ))
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
#endif
ENDDO ! END T and QV LOOP OVER OBS
ENDIF !end if(istemp.eq.1 .or. ismois.eq.1)
!ajb 01062008 - ISPSTR is not used
!!**********************************************************
!! PROCESS SURFACE PRESSURE CROSS-POINT FIELD, IVAR=5,
!! USING BILINEAR FOUR-POINT 2-D INTERPOLATION
!!**********************************************************
! IF(ISPSTR.EQ.1) THEN
! DO N=1,NSTA
! IOB=MAX0(1,IA(N))
! IOB=MIN0(IOB,ide-1)
! JOB=MAX0(1,IB(N))
! JOB=MIN0(JOB,jde-1)
!#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
!#endif
! DXOB=RA(N)
! DYOB=RB(N)
!!ajb fix this (put in correct pressure calc for IOB,JOB here)
! ERRF(5,N)=ERRF(5,N)+VAROBS(5,N)-((1.-DyOB)*((1.-DxOB)* &
! pbase(IOB,1,JOB)+DxOB*pbase(IOB+1,1,JOB))+DyOB* &
! ((1.-DxOB)*pbase(IOB,1,JOB+1)+DxOB* &
! pbase(IOB+1,1,JOB+1)))
!
!#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
!#endif
! ENDDO
! ENDIF ! end if(ispstr.eq.1)
!!**********************************************************
ENDIF ! end if(ityp.eq.3)
ENDDO ! END BIG LOOP
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
! Gather the errf values calculated by the processors with the obs.
CALL get_full_obs_vector(nsta, nerrf, niobf, mp_local_uobmask, &
mp_local_vobmask, mp_local_cobmask, errf)
#endif
! DIFFERENCE BETWEEN OBS AND FCST IS COMPLETED
IF(INEST.EQ.1)THEN
INPF=NPFI
ELSE
FNPF=IRATIO**LEVIDN(INEST)
INPF=FNPF*NPFI
ENDIF
! Gross error check for temperature. Set all vars bad.
do n=1,nsta
if((abs(errf(3,n)).gt.20.).and. &
(errf(3,n).gt.-800000.))then
errf(1,n)=-888888.
errf(2,n)=-888888.
errf(3,n)=-888888.
errf(4,n)=-888888.
varobs(1,n)=-888888.
varobs(2,n)=-888888.
varobs(3,n)=-888888.
varobs(4,n)=-888888.
endif
enddo
! For printout
! IF(MOD(KTAU,INPF).NE.0) THEN
! RETURN
! ENDIF
RETURN
END SUBROUTINE errob
SUBROUTINE upoint (docs) (i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, & 1
arrin, arrout)
!------------------------------------------------------------------------------
! PURPOSE: This subroutine interpolates a real 2D array defined over mass
! coordinate points, to U (momentum) points.
!
!------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
INTEGER, INTENT(IN) :: ids ! Starting i index for entire model domain
INTEGER, INTENT(IN) :: ide ! Ending i index for entire model domain
INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on U points
! Local variables
integer :: i, j
! Do domain interior first
do j = j_start, j_end
do i = max(2,i_start), i_end
arrout(i,j) = 0.5*(arrin(i,j)+arrin(i-1,j))
enddo
enddo
! Do west-east boundaries
if(i_start .eq. ids) then
do j = j_start, j_end
arrout(i_start,j) = arrin(i_start,j)
enddo
endif
if(i_end .eq. ide-1) then
do j = j_start, j_end
arrout(i_end+1,j) = arrin(i_end,j)
enddo
endif
RETURN
END SUBROUTINE upoint
SUBROUTINE vpoint (docs) (i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, & 1
arrin, arrout)
!------------------------------------------------------------------------------
! PURPOSE: This subroutine interpolates a real 2D array defined over mass
! coordinate points, to V (momentum) points.
!
!------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
INTEGER, INTENT(IN) :: jds ! Starting j index for entire model domain
INTEGER, INTENT(IN) :: jde ! Ending j index for entire model domain
INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on V points
! Local variables
integer :: i, j
! Do domain interior first
do j = max(2,j_start), j_end
do i = i_start, i_end
arrout(i,j) = 0.5*(arrin(i,j)+arrin(i,j-1))
enddo
enddo
! Do south-north boundaries
if(j_start .eq. jds) then
do i = i_start, i_end
arrout(i,j_start) = arrin(i,j_start)
enddo
endif
if(j_end .eq. jde-1) then
do i = i_start, i_end
arrout(i,j_end+1) = arrin(i,j_end)
enddo
endif
RETURN
END SUBROUTINE vpoint
LOGICAL FUNCTION TILE_MASK (docs) (iloc, jloc, its, ite, jts, jte) 5
!------------------------------------------------------------------------------
! PURPOSE: Check to see if an i, j grid coordinate is in the tile index range.
!
! Returns: TRUE if the grid coordinate (ILOC,JLOC) is in the tile defined by
! tile-range indices (its,jts) and (ite,jte)
! FALSE otherwise.
!
!------------------------------------------------------------------------------
IMPLICIT NONE
INTEGER, INTENT(IN) :: iloc
INTEGER, INTENT(IN) :: jloc
INTEGER, INTENT(IN) :: its
INTEGER, INTENT(IN) :: ite
INTEGER, INTENT(IN) :: jts
INTEGER, INTENT(IN) :: jte
! Local variables
LOGICAL :: retval
TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND. &
jloc .LE. jte .AND. jloc .GE. jts )
RETURN
END FUNCTION TILE_MASK
!-----------------------------------------------------------------------
SUBROUTINE nudob (docs) (j, ivar, aten, inest, ifrest, ktau, ktaur, & 4,12
xtime, mu, msfx, msfy, nndgv, nerrf, niobf, maxdom, &
npfi, ionf, rinxy, twindo, levidn, &
parid, nstat, &
fdob, lev_in_ob, plfo, nlevs_ob, &
iratio, dx, dtmin, rio, rjo, rko, &
timeob, varobs, errf, pbase, ptop, pp, &
iswind, istemp, ismois, giv, git, giq, &
savwt, kpblt, nscan, &
iprt, &
ids,ide, jds,jde, kds,kde, & ! domain dims
ims,ime, jms,jme, kms,kme, & ! memory dims
its,ite, jts,jte, kts,kte ) ! tile dims
!-----------------------------------------------------------------------
USE module_model_constants
USE module_domain
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
!
! PURPOSE: THIS SUBROUTINE GENERATES NUDGING TENDENCIES FOR THE J-TH
! VERTICAL SLICE (I-K PLANE) FOR FOUR-DIMENSIONAL DATA
! ASSIMILATION FROM INDIVIDUAL OBSERVATIONS. THE NUDGING
! TENDENCIES ARE FOUND FROM A ONE-PASS CALCULATION OF
! WEIGHTING FACTORS SIMILAR TO THE BENJAMIN-SEAMAN OBJECTIVE
! ANALYSIS. THIS SUBROUTINE IS DESIGNED FOR RAPID EXECUTION
! AND MINIMAL STORAGE REQUIREMENTS. ALGORITHMS SHOULD BE
! VECTORIZED WHEREVER POSSIBLE.
!
! HISTORY: Original author: MM5 version???
! 02/04/2004 - Creation of WRF version. Al Bourgeois
! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
!------------------------------------------------------------------------------
!
! NOTE: This routine was originally designed for MM5, which uses
! a nonstandard (I,J) coordinate system. For WRF, I is the
! east-west running coordinate, and J is the south-north
! running coordinate. So a "J-slab" here is west-east in
! extent, not south-north as for MM5. -ajb 06/10/2004
!
! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
!
! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
! TYPES OF FACTORS:
! 1) TIME WEIGHTING - ONLY OBSERVATIONS WITHIN A SELECTED
! TIME WINDOW (TWINDO) CENTERED AT THE CURRENT FORECAST
! TIME (XTIME) ARE USED. OBSERVATIONS CLOSEST TO
! XTIME ARE TIME-WEIGHTED MOST HEAVILY (TIMEWT)
! 2) VERTICAL WEIGHTING - NON-ZERO WEIGHTS (WTSIG) ARE
! CALCULATED WITHIN A VERTICAL REGION OF INFLUENCE
! (RINSIG).
! 3) HORIZONTAL WEIGHTING - NON-ZERO WEIGHTS (WTIJ) ARE
! CALCULATED WITHIN A RADIUS OF INFLUENCE (RINXY). THE
! VALUE OF RIN IS DEFINED IN KILOMETERS, AND CONVERTED
! TO GRID LENGTHS FOR THE APPROPRIATE MESH SIZE.
!
! THE FIVE FORECAST VARIABLES ARE PROCESSED BY CHANGING THE
! VALUE OF IVAR AS FOLLOWS:
! IVAR VARIABLE(TAU-1)
! ---- ---------------
! 1 U
! 2 V
! 3 T
! 4 QV
! 5 PSB(CROSS) REMOVED IN V3
! (6) PSB(DOT)
!
!-----------------------------------------------------------------------
!
! Description of input arguments.
!
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
INTEGER, INTENT(IN) :: j ! south-north running coordinate.
INTEGER, INTENT(IN) :: ivar
INTEGER, INTENT(IN) :: inest ! domain index
LOGICAL, INTENT(IN) :: ifrest
INTEGER, INTENT(IN) :: ktau
INTEGER, INTENT(IN) :: ktaur
REAL, INTENT(IN) :: xtime ! forecast time in minutes
INTEGER, INTENT(IN) :: nndgv ! number of nudge variables
INTEGER, INTENT(IN) :: nerrf ! number of error fields
INTEGER, INTENT(IN) :: niobf ! number of observations
INTEGER, INTENT(IN) :: maxdom ! maximum number of domains
INTEGER, INTENT(IN) :: npfi
INTEGER, INTENT(IN) :: ionf
REAL, INTENT(IN) :: rinxy
REAL, INTENT(IN) :: twindo
INTEGER, INTENT(IN) :: levidn(maxdom) ! level of nest
INTEGER, INTENT(IN) :: parid(maxdom) ! parent domain id
INTEGER, INTENT(IN) :: nstat ! number of obs stations
TYPE(fdob_type), intent(inout) :: fdob
REAL, INTENT(IN) :: lev_in_ob(niobf) ! Level in sounding-type obs.
REAL, intent(IN) :: plfo(niobf)
REAL, INTENT(IN) :: nlevs_ob(niobf) ! Number of levels in sounding.
INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
REAL, INTENT(IN) :: dx ! This domain grid cell-size (m)
REAL, INTENT(IN) :: dtmin
REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
REAL, INTENT(INOUT) :: rko(niobf) ! Obs vertical coordinate.
REAL, INTENT(IN) :: timeob(niobf)
REAL, INTENT(IN) :: varobs(nndgv,niobf)
REAL, INTENT(IN) :: errf(nerrf, niobf)
REAL, INTENT(IN) :: pbase( ims:ime, kms:kme ) ! Base pressure.
REAL, INTENT(IN) :: ptop
REAL, INTENT(IN) :: pp( ims:ime, kms:kme ) ! Pressure perturbation (Pa)
REAL, INTENT(IN) :: mu(ims:ime) ! Air mass on u, v, or mass-grid
REAL, INTENT(IN) :: msfx(ims:ime) ! Map scale (only used for vars u & v)
REAL, INTENT(IN) :: msfy(ims:ime) ! Map scale (only used for vars u & v)
INTEGER, intent(in) :: iswind ! Nudge flag for wind
INTEGER, intent(in) :: istemp ! Nudge flag for temperature
INTEGER, intent(in) :: ismois ! Nudge flag for moisture
REAL, intent(in) :: giv ! Coefficient for wind
REAL, intent(in) :: git ! Coefficient for temperature
REAL, intent(in) :: giq ! Coefficient for moisture
REAL, INTENT(INOUT) :: aten( ims:ime, kms:kme)
REAL, INTENT(INOUT) :: savwt( nndgv, ims:ime, kms:kme )
INTEGER, INTENT(IN) :: kpblt(its:ite)
INTEGER, INTENT(IN) :: nscan ! number of scans
LOGICAL, INTENT(IN) :: iprt ! print flag
! Local variables
integer :: mm(maxdom)
integer :: kobs ! k-lev below obs (for obs straddling pblt)
real :: ra(niobf)
real :: rb(niobf)
real :: psurf(niobf)
real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
real :: rscale(ims:ime) ! For converting to rho-coupled units.
! real :: tscale(ims:ime,kms:kme) ! For converting to potential temp. units.
real :: reserf(100)
character*40 name
character*3 chr_hr
character(len=200) :: msg ! Argument to wrf_message
!*** DECLARATIONS FOR IMPLICIT NONE
integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
integer :: komin,komax,nn,nhi,nlo,nnjc
integer :: i_s,i_e
integer :: istq
real :: gfactor,rfactor,gridx,gridy,rindx,schnes,ris
real :: grfacx,grfacy
real :: fdtim,tw1,tw2,tconst,timewt,timewt2,ttim,dift,pob
real :: ri,rj,rx,ry,rsq,wtij,pdfac,erfivr,dk,slope,rinfac
real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
real :: scratch
! print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
! if(ivar.ne.4)return
!yliu start -- for multi-scans: NSCAN=0: original
! NSCAN=1: added a scan with a larger Ri and smaller G
! if(NSCAN.ne.0 .and. NSCAN.ne.1) stop
! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
if(NSCAN.ne.0) then
IF (iprt) then
write(msg,*) 'SAVWT must be resized for NSCAN=1'
call wrf_message(msg)
ENDIF
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
endif
IPLO=0 + NSCAN*4
GFACTOR=1. + NSCAN*(-1. + 0.33333)
RFACTOR=1. + NSCAN*(-1. + 3.0)
!yliu end
! jc
! return if too close to j boundary
if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
! write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
! $ ' too close to boundary.'
return
endif
if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
! write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
! $ ' too close to boundary.'
return
endif
! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
ICUT=0
IF(INEST.GT.1)ICUT=1
i_s = max0(2+icut,its)
i_e = min0(ide-1-icut,ite)
IPL=IVAR + IPLO !yliu +IPLO
! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
INPF=(IRATIO**LEVIDN(INEST))*NPFI
INFR=(IRATIO**LEVIDN(INEST))*IONF
GRIDX=0.0
GRIDY=0.0
IGRID=0
IF(IVAR.GE.3)THEN
GRIDX=0.5
GRIDY=0.5
IGRID=1
ELSEIF(IVAR.eq.1) THEN
GRIDY=0.5
GRIDX=0.0
IGRID=1
ELSEIF(IVAR.eq.2) THEN
GRIDX=0.5
GRIDY=0.0
IGRID=1
ENDIF
! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
! KILOMETERS TO GRID LENGTHS, RINDX
RINDX=RINXY*1000./DX * RFACTOR !yliu *RFACTOR
! jc
! make horizontal radius vary per nest
! rindx=rindx/float(inest)
! yliu test1 -- En 3, 4
! rindx=rindx/float(3**(in-1)) !YLIU
! jc
! make horizontal radius vary per nest
! schnes=1/float(inest) !JC
! yliu test1 -- En 3, 4 !YLIU
schnes=1/float(3**(inest-1)) !JC
! reduce the Rinf in the nested grid proportionally
rindx=rindx*schnes
! rinfmn =1., rinfmx=2., pfree=50 in param.F
! yliu test: for upper-air data, use larger influence radii
! Essentially increase the slope -- the same radii
! at 500mb and above as the coarse mesh and the
! same small radii at sfc as that for sfc obs
fdob%rinfmx=2. *1.0 /schnes !YLIU
! rinfmx=1.2*1.0 /schnes !YLIU
! jc
RIS=RINDX*RINDX
IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
IF(MOD(KTAU,INFR).NE.0)GOTO 126
5 CONTINUE
IF (iprt) THEN
IF(J.EQ.10) then
write(msg,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
call wrf_message(msg)
ENDIF
ENDIF
6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,', &
'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2, &
' rindx=',f4.1)
!********************************************************************
!ajb_07012008 Setting ra and rb for the fine mesh her is now simple.
! Values are no longer calculated here based on the
! coarse mesh, since direct use of WRF map projections
! on each nest was implemented in subroutine in4dob.
!********************************************************************
! SET RA AND RB
DO N=1,NSTAT
RA(N)=RIO(N)-GRIDX
RB(N)=RJO(N)-GRIDY
ENDDO
! OUTPUT OBS PER GRID EVERY HOUR
if ( mod(xtime,60.).gt.56. .and. ivar.eq.3 .and. j.eq.10) then
! IF (iprt) print *,'outputting obs number on grid ', &
! inest,' at time=',xtime
! write(chr_hr(1:3),'(i3)')nint(xtime/60.)
! if(chr_hr(1:1).eq.' ')chr_hr(1:1)='0'
! if(chr_hr(2:2).eq.' ')chr_hr(2:2)='0'
! IF (iprt) print *,'chr_hr=',chr_hr(1:3),nint(xtime/60.)
! open(91,file= &
! 'obs_g'//char(inest+ichar('0'))//'_'//chr_hr(1:3), &
! form='FORMATted',status='unknown')
! write(91,911)nstat
! write(msg,911)nstat
! call wrf_message(msg)
!911 FORMAT('total obs=',i8)
nsthis=0
nsmetar=0
nsspeci=0
nsship=0
nssynop=0
nstemp=0
nspilot=0
nssatwnds=0
nssams=0
nsprofs=0
! print *,'ide,jde=',ide,jde
do jjjn=1,nstat
! DETERMINE THE TIME-WEIGHT FACTOR FOR N
FDTIM=XTIME-DTMIN
! CONVERT TWINDO AND TIMEOB FROM HOURS TO MINUTES:
TW1=TWINDO/2.*60.
TW2=TWINDO*60.
TCONST=1./TW1
TIMEWT2=0.0
TTIM=TIMEOB(jjjn)*60.
!***********TTIM=TARGET TIME IN MINUTES
DIFT=ABS(FDTIM-TTIM)
IF(DIFT.LE.TW1)TIMEWT2=1.0
IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
IF(FDTIM.LT.TTIM)TIMEWT2=(FDTIM-(TTIM-TW2))*TCONST
IF(FDTIM.GT.TTIM)TIMEWT2=((TTIM+TW2)-FDTIM)*TCONST
ENDIF
! print *,'timewt2=',timewt2,ttim,fdtim
if (ra(jjjn).ge.1. .and. rb(jjjn).ge.1. &
.and.ra(jjjn).le.real(ide) .and. rb(jjjn).le.real(jde) &
.and.timewt2.gt.0.) then
if(lev_in_ob(jjjn).eq.1.)nsthis=nsthis+1
if(plfo(jjjn).eq.1.)nsmetar=nsmetar+1
if(plfo(jjjn).eq.2.)nsspeci=nsspeci+1
if(plfo(jjjn).eq.3.)nsship=nsship+1
if(plfo(jjjn).eq.4.)nssynop=nssynop+1
if(plfo(jjjn).eq.5..and.lev_in_ob(jjjn).eq.1.) nstemp=nstemp+1
if(plfo(jjjn).eq.6..and.lev_in_ob(jjjn).eq.1.) nspilot=nspilot+1
if(plfo(jjjn).eq.7.)nssatwnds=nssatwnds+1
if(plfo(jjjn).eq.8.)nssams=nssams+1
if(plfo(jjjn).eq.9..and.lev_in_ob(jjjn).eq.1.) nsprofs=nsprofs+1
endif
enddo
! write(91,912)nsthis
! write(msg,912)nsthis
! call wrf_message(msg)
!912 FORMAT('total obs on this grid=',i8)
! write(91,921)nsmetar
! write(msg,921)nsmetar
! call wrf_message(msg)
!921 FORMAT('total metar obs on this grid=',i8)
! write(91,922)nsspeci
! write(msg,922)nsspeci
! call wrf_message(msg)
!922 FORMAT('total special obs on this grid=',i8)
! write(91,923)nsship
! write(msg,923)nsship
! call wrf_message(msg)
!923 FORMAT('total ship obs on this grid=',i8)
! write(91,924)nssynop
! write(msg,924)nssynop
! call wrf_message(msg)
!924 FORMAT('total synop obs on this grid=',i8)
! write(91,925)nstemp
! write(msg,925)nstemp
! call wrf_message(msg)
!925 FORMAT('total temp obs on this grid=',i8)
! write(91,926)nspilot
! write(msg,926)nspilot
! call wrf_message(msg)
!926 FORMAT('total pilot obs on this grid=',i8)
! write(91,927)nssatwnds
! write(msg,927)nssatwnds
! call wrf_message(msg)
!927 FORMAT('total sat-wind obs on this grid=',i8)
! write(91,928)nssams
! write(msg,928)nssams
! call wrf_message(msg)
!928 FORMAT('total sams obs on this grid=',i8)
! write(91,929)nsprofs
! write(msg,929)nsprofs
! call wrf_message(msg)
!929 FORMAT('total profiler obs on this grid=',i8)
! close(91)
endif ! END OUTPUT OBS PER GRID EVERY HOUR
! INITIALIZE WEIGHTING ARRAYS TO ZERO
DO I=its,ite
DO K=1,kte
WT(I,K)=0.0
WT2ERR(I,K)=0.0
ENDDO
ENDDO
! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
!
! COMPUTE P* AT OBS LOCATION (RA,RB). DO THIS AS SEPARATE VECTOR LOOP H
! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
DO N=1,NSTAT
IF(IVAR.GE.3)THEN
PSURF(N)=ERRF(6,N)
ELSE
IF(IVAR.EQ.1)THEN
PSURF(N)=ERRF(7,N) ! U-points
ELSE
PSURF(N)=ERRF(8,N) ! V-points
ENDIF
ENDIF
ENDDO
! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
! J-STRIP
MAXJ=J+IFIX(RINDX*fdob%RINFMX+0.99) !ajb
MINJ=J-IFIX(RINDX*fdob%RINFMX+0.99) !ajb
! jc comment out this? want to use obs beyond the domain?
! MAXJ=MIN0(JL-IGRID,MAXJ) !yliu
! MINJ=MAX0(1,MINJ) !yliu
n=1
!***********************************************************************
DO nnn=1,NSTAT ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
!***********************************************************************
! Soundings are consecutive obs, but they need to be treated as a single
! entity. Thus change the looping to nnn, with n defined separately.
!yliu
! note for sfc data: nsndlev=1 and njcsnd=1
nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
! yliu start -- set together with the other parts
! test: do the sounding levels as individual obs
! nsndlev=1
! yliu end
njcsnd=nsndlev
! set pob here, to be used later
pob=varobs(5,n)
! print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
! print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP. THIS
! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
! ATION.
!yliu: Skip bad obs if it is sfc or single level sounding.
!yliu: Before this (020918), a snd will be skipped if its first
!yliu level has bad data- often true due to elevation
IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
! print *, " bad obs skipped"
ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
! print *, " skipped obs far away from this J-slice"
!----------------------------------------------------------------------
ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
!----------------------------------------------------------------------
! DETERMINE THE TIME-WEIGHT FACTOR FOR N
FDTIM=XTIME-DTMIN
! TWINDO IS IN MINUTES:
TW1=TWINDO/2.*60.
TW2=TWINDO*60.
TCONST=1./TW1
TIMEWT=0.0
TTIM=TIMEOB(N)*60.
!***********TTIM=TARGET TIME IN MINUTES
DIFT=ABS(FDTIM-TTIM)
IF(DIFT.LE.TW1)TIMEWT=1.0
IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
ENDIF
! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
! FOR THE VERTICAL WEIGHTING, WTSIG
! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
!ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
!ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
rko(n) = errf(9,n) !ajb 20021210
#endif
KOB=nint(RKO(N))
KOB=MIN0(kte,KOB)
KOB=MAX0(1,KOB)
! ASSIMILATE SURFACE LAYER DATA ON SIGMA
IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
DO K=1,kte
WTSIG(K)=0.0
ENDDO
! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
! WTSIG(1)=1.0
! WTSIG(2)=0.67
! WTSIG(3)=0.33
! KOMIN=3
! KOMAX=1
! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
! fix this because kpblt at 1 and il is 0
MAXI=IFIX(RA(N)+0.99+RINDX)
MAXI=MIN0(ide-1,MAXI)
MINI=IFIX(RA(N)-RINDX-0.99)
MINI=MAX0(2,MINI)
!yliu start
! use also obs outside of this domain -- surface obs
! if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
! & RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
! print *, " skipped obs far away from this domain"
! currently can use obs within this domain or ones very close to (1/3
! influence of radius in the coarse domain) this
! domain. In later case, use BC column value to approximate the model value
! at obs point -- ERRF need model field in errob.F !!
if ( RA(N).GE.(0.-RINDX/3) &
.and. RA(N).LE.float(ide)+RINDX/3 &
.and. RB(N).GE.(0.-RINDX/3) &
.and. RB(N).LE.float(jde)+RINDX/3) then
! or use obs within this domain only
! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
! print *, " skipped obs far outside of this domain"
! if(j.eq.3 .and. ivar.eq.3) then
! write(6,*) 'N = ',n,' exit 120 3'
! endif
!yliu end
!
! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
RJ=FLOAT(J)
RX=RJ-RB(N)
! WEIGHTS FOR THE 3-D VARIABLES
ERFIVR=ERRF(IVAR,N)
!
!JM I will be local, because it indexes into PDOC, WT, and others
! if((ivar.eq.1 .or. ivar.eq.3) .and. n.le.200) then
! write(6,'(a,i3,a,i3)')'SURF OBS NEAR: N = ',n,' nest = ',inest
! write(6,'(a,f10.3,a,f10.3,a,i3,a,i3,a,i3,a,i2)')
! $ ' RA =',RA(N),' RB =',RB(N),' J =',j,
! $ ' MINI =',MINI,' MAXI =',MAXI,' NEST =',inest
! endif
DO I=max0(its,MINI),min0(ite,MAXI)
RI=FLOAT(I)
RY=RI-RA(N)
RIS=RINDX*RINDX
RSQ=RX*RX+RY*RY
! DPRIM=SQRT(RSQ)
! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
! D=DPRIM+RINDX*DCON*ABS(PSBO(N)-PDOC(I,J))
! DSQ=D*D
! WTIJ=(RIS-DSQ)/(RIS+DSQ)
wtij=(ris-rsq)/(ris+rsq)
scratch = (abs(psurf(n)-.001*pbase(i,1))*fdob%DCON)
pdfac=1.-AMIN1(1.0,scratch)
wtij=wtij*pdfac
WTIJ=AMAX1(0.0,WTIJ)
! try making sfc obs weighting go thru pbl
! jc kpbl is at dot or cross only - need to interpolate?
! wtsig(1)=1.
komax=max0(3,kpblt(i))
! jc arbitrary check here
IF (iprt) THEN
if (kpblt(i).gt.25 .and. ktau.ne.0) then
write(msg,552)inest,i,j,kpblt(i)
call wrf_message(msg)
endif
ENDIF
552 FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
if(kpblt(i).gt.25) komax=3
komin=1
dk=float(komax)
do k=komin,komax
wtsig(k)=float(komax-k+1)/dk
WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ
WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ*WTSIG(K) &
*WTSIG(K)*ERFIVR
! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
!
! write(6,'(a,i2,a,f8.3,a,f8.3,a,f8.3,a,f8.3,a,f8.3)')
! 'Surface obs, after: k = ',k, &
! ' WT2ERR = ',wt2err(i,k), &
! ' TIMEWT = ',timewt, &
! ' WTIJ = ',wtij, &
! ' WSIG = ',wtsig(k), &
! ' ERFIVR = ',erfivr
! endif
enddo
ENDDO
! print *, " Surface "
endif ! end check for obs in domain
! END SURFACE-LAYER U OR V OBS NUDGING
ELSE
! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
!
! print *,'in upper air section'
! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
!ajb SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
slope = (fdob%RINFMN-fdob%RINFMX)/(psurf(n)-fdob%PFREE)
RINFAC=SLOPE*POB+fdob%RINFMX-SLOPE*fdob%pfree
RINFAC=AMAX1(RINFAC,fdob%RINFMN)
RINFAC=AMIN1(RINFAC,fdob%RINFMX)
!yliu: for multilevel upper-air data, take the maximum
! for the I loop.
if(nsndlev.gt.1) RINFAC = fdob%RINFMX
!yliu end
MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
MAXI=MIN0(ide-IGRID,MAXI)
MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
MINI=MAX0(1,MINI)
! yliu start
! use also obs outside of but close to this domain -- upr data
! if( RA(N).LT.(0.-RINFAC*RINDX)
! & .or. RA(N).GT.float(IL)+RINFAC*RINDX
! & .or. RB(N).LT.(0.-RINFAC*RINDX)
! & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then
! print *, " skipped obs far away from this I-range"
! currently can use obs within this domain or ones very close to (1/3
! influence of radius in the coarse domain) this
! domain. In later case, use BC column value to approximate the model value
! at obs point -- ERRF need model field in errob.F !!
!cc if (i.eq.39 .and. j.eq.34) then
!cc write(6,*) 'RA(N) = ',ra(n)
!cc write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
!cc endif
if( RA(N).GE.(0.-RINFAC*RINDX/3) &
.and. RA(N).LE.float(ide)+RINFAC*RINDX/3 &
.and. RB(N).GE.(0.-RINFAC*RINDX/3) &
.and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
! or use obs within this domain only
! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
! print *, " skipped obs far outside of this domain"
! yliu end
! is this 2 needed here - kpbl not used?
! MINI=MAX0(2,MINI)
! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
RJ=FLOAT(J)
RX=RJ-RB(N)
! WEIGHTS FOR THE 3-D VARIABLES
!
ERFIVR=ERRF(IVAR,N)
! jc
nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
! yliu start
! test: do the sounding levels as individual obs
! nsndlev=1
! yliu end
njcsnd=nsndlev
!
DO I=max0(its,MINI),min0(ite,MAXI)
! jc
RI=FLOAT(I)
RY=RI-RA(N)
RIS=RINDX*RINFAC*RINDX*RINFAC
RSQ=RX*RX+RY*RY
! yliu test: for upper-air data, keep D1 influence radii
! RIS=RIS /schnes /schnes
WTIJ=(RIS-RSQ)/(RIS+RSQ)
WTIJ=AMAX1(0.0,WTIJ)
! weight ob in vertical with +- 50 mb
! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
if(nsndlev.eq.1) then
rinprs=7.5
else
rinprs=3.0
endif
! yliu end
!
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
if(nsndlev.eq.1)then
!----------------------------------------------------------------------
! --- HANDLE 1-LEVEL OBSERVATIONS ---
!----------------------------------------------------------------------
! if(I.eq.MINI) print *, " Single snd "
! ERFIVR is the residual (difference) between the ob and the model
! at that point. We can analyze that residual up and down.
! First find komin for ob.
!yliu start -- in the old code, komax and komin were reversed!
do k=kte,1,-1
pijk = .001*(pbase(i,k)+pp(i,k))
! print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
if(pijk.ge.(pob+rinprs)) then
komin=k
go to 325
endif
enddo
komin=1
325 continue
! now find komax for ob
do k=3,kte
pijk = .001*(pbase(i,k)+pp(i,k))
if(pijk.le.(pob-rinprs)) then
komax=k
go to 326
endif
enddo
komax=kte ! yliu 20050706
326 continue
! yliu: single-level upper-air data will act either above or below the PBL top
! ajb: Reset komin or komax. if kobs>kpblt, komin=kpblt+1, else komax=kpblt
if( (kpblt(i).le.komax) .and. (kpblt(i).ge.komin) ) then
kobs = 1
OBS_K: do k = komin, komax
if( pob .gt. .001*(pbase(i,k)+pp(i,k)) ) then
kobs = k
EXIT OBS_K
endif
enddo OBS_K
if(kobs.gt.kpblt(i)) then
komin=max0(kobs, komin) ! kobs here is kpblt(i)+1
else
komax=min0(kpblt(i), komax)
endif
endif
! yliu end
!
! print *,'1 level, komin,komax=',komin,komax
! if(i.eq.MINI) then
! print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
! ERFIVR=0
! endif
do k=1,kte
reserf(k)=0.0
wtsig(k)=0.0
enddo
!yliu end
!cc if (i.eq.39 .and. j.eq.34) then
!cc write(6,*) ' komin = ', komin,' komax = ',komax
!cc endif
do k=komin,komax
pijk = .001*(pbase(i,k)+pp(i,k))
reserf(k)=erfivr
wtsig(k)=1.-abs(pijk-pob)/rinprs
wtsig(k)=amax1(wtsig(k),0.0)
! print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
! Now calculate WT and WT2ERR for each i,j,k point cajb
WT(I,K)=WT(I,K)+TIMEWT*WTIJ*wtsig(k)
WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ* &
reserf(k)*wtsig(k)*wtsig(k)
enddo
else
!----------------------------------------------------------------------
! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
!----------------------------------------------------------------------
!yliu start
! if(I.eq.MINI) print *, " Multi-level snd "
! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
IF (iprt) THEN
write(msg,*) "n = ",n,"nsndlev = ",nsndlev
call wrf_message(msg)
write(msg,*) "nlevs_ob,lev_in_ob", &
nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
call wrf_message(msg)
call wrf_message("in nudobs.F: sounding level messed up, stopping")
ENDIF
call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
endif
!yliu end
! This is for a multi-level observation
! The trick here is that the sounding is "one ob". You don't
! want multiple levels to each be treated like separate
! and independent observations.
! At each i,j want to interpolate sounding to the model levels at that
! particular point.
komin=1
komax=kte-2
! this loop goes to 1501
! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
!yliu start
do k=1,kte
reserf(k)=0.0
wtsig(k)=0.0
enddo
!yliu end
do k=komax,komin,-1
pijk = .001*(pbase(i,k)+pp(i,k))
! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
if(pijk.gt.varobs(5,n)) then
go to 1501
endif
! if sigma level pressure is .lt. than the highest ob level, don't interpolate
if(pijk.le.varobs(5,n+nsndlev-1)) then
go to 1501
endif
! now interpolate sounding to this k
! yliu start-- recalculate WTij for each k-level
!ajb SLOPE = (fdob%RINFMN-fdob%RINFMX)/(pdoc(i,j)+PTOP-fdob%PFREE)
slope = (fdob%RINFMN-fdob%RINFMX)/ (.001*pbase(i,1)-fdob%PFREE)
RINFAC=SLOPE*pijk+fdob%RINFMX-SLOPE*fdob%PFREE
RINFAC=AMAX1(RINFAC,fdob%RINFMN)
RINFAC=AMIN1(RINFAC,fdob%RINFMX)
RIS=RINDX*RINFAC*RINDX*RINFAC
RSQ=RX*RX+RY*RY
! for upper-air data, keep D1 influence radii
! RIS=RIS /schnes /schnes
WTIJ=(RIS-RSQ)/(RIS+RSQ)
WTIJ=AMAX1(0.0,WTIJ)
! yliu end
! this loop goes to 1503
do nn=2,nsndlev
! only set pobhi if varobs(ivar) is ok
pobhi=-888888.
if(varobs(ivar,n+nn-1).gt.-800000. &
.and. varobs(5,n+nn-1).gt.-800000.) then
pobhi=varobs(5,n+nn-1)
nhi=n+nn-1
if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.20.) then
go to 1502 ! within 200mb of obs height
endif
endif
enddo
! did not find any ob above within 100 mb, so jump out
go to 1501
1502 continue
nlo=nhi-1
do nnjc=nhi-1,n,-1
if(varobs(ivar,nnjc).gt.-800000. &
.and. varobs(5,nnjc).gt.-800000.) then
poblo=varobs(5,nnjc)
nlo=nnjc
if(poblo.gt.pijk .and. abs(poblo-pijk).lt.20.) then
go to 1505 ! within 200mb of obs height
endif
endif
enddo
!yliu end --
! did not find any ob below within 200 mb, so jump out
go to 1501
1505 continue
! interpolate to model level
pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
reserf(k)=errf(ivar,nlo)+ &
(errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
wtsig(k)=1.
1501 continue
! now calculate WT and WT2ERR for each i,j,k point cajb
WT(I,K)=WT(I,K)+TIMEWT*WTIJ*wtsig(k)
WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ*WTIJ* &
reserf(k)*wtsig(k)*wtsig(k)
! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
!
! if(wt(i,k) .ne. 0.0) then
! scratch = WT2ERR(I,K)/WT(I,K)
! else
! scratch = 999.
! endif
!
! write(6,'(a,i2,a,f8.3,a,f4.2,a,f7.4,a,f4.2,a,f5.3,a,f7.4)')
! $ 'Multi-level obs: k = ',k,
! $ ' WT2ERR = ',wt2err(i,k),
! $ ' WTIJ = ',wtij,
! $ ' RSF = ',reserf(k),
! $ ' WSIG = ',wtsig(k),
! $ ' WT = ',wt(i,k),
! $ ' W2EOWT = ',scratch
! endif
! end do k
enddo ! enddo k levels
! end multi-levels
endif ! end if(nsndlev.eq.1)
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
!
ENDDO ! END DO MINI,MAXI LOOP
endif ! check for obs in domain
! END OF NUDGING TO OBS ON PRESSURE LEVELS
ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
!----------------------------------------------------------------------
ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
!----------------------------------------------------------------------
! n=n+1
n=n+njcsnd
!yliu 1202 continue
if(n.gt.nstat)then
! print *,'n,nstat=',n,nstat,ivar,j
go to 1203
endif
! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
!***********************************************************************
ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
!***********************************************************************
1203 continue
! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
! THE ATEN ARRAY
! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
! THEY ARE USED BELOW IN THE DENOMINATOR.
DO K=kts,kte
DO I=its,ite
IF(WT(I,K).EQ.0)THEN
WT2ERR(I,K)=0.0
ENDIF
IF(WT(I,K).EQ.0)THEN
WT(I,K)=1.0
ENDIF
ENDDO
ENDDO
126 CONTINUE
IF(IVAR.GE.3)GOTO 170
! this is for u,v
! 3-D DOT POINT TENDENCIES
! Calculate scales for converting nudge factor from u (v)
! to rho_u (or rho_v) units.
IF (IVAR == 1) THEN
call calc_rcouple_scales(mu,msfy,rscale,ims,ime,its,ite)
ELSE IF (IVAR == 2) THEN
call calc_rcouple_scales(mu,msfx,rscale,ims,ime,its,ite)
END IF
DO K=1,kte
DO I=i_s,i_e
IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
W2EOWT=WT2ERR(I,K)/WT(I,K)
ELSE
W2EOWT=SAVWT(IPL,I,K)
ENDIF
! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
! scratch = GIV*RSCALE(I)*W2EOWT*fdob%TFACI*ISWIND*GFACTOR
! write(6,*) 'ATEN calc: k = ',k
! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
! $ ' W2EOWT = ',w2eowt
! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
! $ ' GFACTOR = ',gfactor
! endif
!
! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
! scratch = GIV*RSCALE(I)*W2EOWT*fdob%TFACI*ISWIND*GFACTOR
! write(6,*) 'ATEN calc: k = ',k
! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
! $ ' W2EOWT = ',w2eowt
! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
! $ ' GFACTOR = ',gfactor
! endif
ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
*W2EOWT*fdob%TFACI &
*ISWIND *GFACTOR !yliu *GFACTOR
! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
! endif
! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
! endif
ENDDO
ENDDO
IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
DO K=1,kte
DO I=its,ite
SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
ENDDO
ENDDO
ENDIF
RETURN
170 CONTINUE
! 3-D CROSS-POINT TENDENCIES
! this is for t (ivar=3) and q (ivsr=4)
IF(3-IVAR.LT.0)THEN
GITQ=GIQ
ELSE
GITQ=GIT
ENDIF
IF(3-IVAR.LT.0)THEN
ISTQ=ISMOIS
ELSE
ISTQ=ISTEMP
ENDIF
DO K=1,kte
DO I=i_s,i_e
IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
W2EOWT=WT2ERR(I,K)/WT(I,K)
ELSE
W2EOWT=SAVWT(IPL,I,K)
ENDIF
! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
! scratch = GITQ*MU(I)*W2EOWT*fdob%TFACI*ISTQ*GFACTOR
! write(6,*) 'ATEN calc: k = ',k
! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
! $ ' W2EOWT = ',w2eowt
! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
! $ ' GFACTOR = ',gfactor
! endif
!
! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
! scratch = GITQ*MU(I)*W2EOWT*fdob%TFACI*ISTQ*GFACTOR
! write(6,*) 'ATEN calc: k = ',k
! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
! $ ' W2EOWT = ',w2eowt
! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
! $ ' GFACTOR = ',gfactor
! endif
ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
*W2EOWT*fdob%TFACI*ISTQ *GFACTOR !yliu *GFACTOR
! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
! endif
! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
! endif
ENDDO
ENDDO
IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
DO K=1,kte
DO I=its,ite
SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
ENDDO
ENDDO
ENDIF
RETURN
END SUBROUTINE nudob
SUBROUTINE calc_rcouple_scales (docs) (a, msf, rscale, ims,ime, its,ite) 2
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
! Local variables
integer :: i
! Calculate scales to be used for producing rho-coupled nudging factors.
do i = its,ite
rscale(i) = a(i)/msf(i)
enddo
RETURN
END SUBROUTINE calc_rcouple_scales
!ajb: Not used
SUBROUTINE set_real_array (docs) (rscale, value, ims,ime, its,ite)
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
REAL, INTENT(IN) :: value ! Constant array value
REAL, INTENT(OUT) :: rscale( ims:ime ) ! Output array
! Local variables
integer :: i
! Set array to constant value
do i = its,ite
rscale(i) = value
enddo
RETURN
END SUBROUTINE set_real_array
!ajb: Not used
SUBROUTINE calc_pottemp_scales (docs) (ivar, rcp, pb, p, tscale, &
ims,ime, its,ite, &
kms,kme, kts,kte)
!-----------------------------------------------------------------------
IMPLICIT NONE
!-----------------------------------------------------------------------
INTEGER, INTENT(IN) :: ims,ime, kms,kme ! Memory dimensions
INTEGER, INTENT(IN) :: its,ite, kts,kte ! Tile dimensions
INTEGER, INTENT(IN) :: ivar ! Variable identifier
REAL, INTENT(IN) :: rcp ! Constant (2./7.)
REAL, INTENT(IN) :: pb(ims:ime, kms:kme) ! Base pressure (Pa) array
REAL, INTENT(IN) :: p(ims:ime, kms:kme) ! Pressure pert. (Pa) array
REAL, INTENT(OUT) :: tscale(ims:ime, kms:kme) ! Scales for pot. temp.
! Local variables
integer :: i,k
if(ivar.eq.3) then
! Calculate scales to be used for producing potential temperature nudging factors.
do k = kts,kte
do i = its,ite
tscale(i,k) = ( 1000000. / ( pb(i,k)+p(i,k)) )**rcp
enddo
enddo
else
! Set to 1. for moisture scaling.
do k = kts,kte
do i = its,ite
tscale(i,k) = 1.0
enddo
enddo
endif
RETURN
END SUBROUTINE calc_pottemp_scales
SUBROUTINE print_obs_info (docs) (iprt,inest,niobf,rio,rjo,rko, & 1,7
prt_max,prt_freq,obs,lat,lon, &
mlat,mlon,timeob,xtime)
!*************************************************************************
! Purpose: Print obs information.
!*************************************************************************
IMPLICIT NONE
LOGICAL, intent(in) :: iprt ! Print flag
INTEGER, intent(in) :: inest ! Nest level
INTEGER, intent(in) :: niobf ! Maximum number of observations
REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
REAL, intent(in) :: rko(niobf) ! Bottom-top north coord (non-stagger)
INTEGER, intent(in) :: prt_max ! Max no. of obs for diagnostic printout
INTEGER, intent(in) :: prt_freq ! Frequency for diagnostic printout
INTEGER, intent(inout) :: obs(prt_max) ! Saved obs indices to print
REAL, intent(inout) :: lat(prt_max) ! Saved latitudes
REAL, intent(inout) :: lon(prt_max) ! Saved longitudes
REAL, intent(inout) :: mlat(prt_max) ! Saved model latitudes
REAL, intent(inout) :: mlon(prt_max) ! Saved longitudes
REAL, intent(in) :: timeob(niobf) ! Times of each observation (hours)
REAL, intent(in) :: xtime ! Model time in minutes
! Local variables
integer :: n ! Loop counter over obs
integer :: pnx ! Obs index for printout
character(len=200) :: msg ! Argument to wrf_message
if(iprt) then
if(prt_max.gt.0) then
if(obs(1).ne.-999) then
call wrf_message("")
write(msg,'(a,i4,a,f8.1,a)') 'REPORTING OBS MASS-PT LOCS FOR NEST ', &
inest,' AT XTIME=',xtime,' MINUTES'
call wrf_message(msg)
write(msg,'(a,i4,a,i5,a)') 'FREQ=',prt_freq,', MAX=',prt_max, &
' LOCS, NEWLY READ OBS ONLY, -999 => OBS OFF PROC'
call wrf_message(msg)
call wrf_message("")
write(msg,'(a,a)') ' OBS# I J K OBS LAT', &
' OBS LON XLAT(I,J) XLONG(I,J) TIME(hrs)'
call wrf_message(msg)
endif
endif
! Note: rio and rjo are referenced to non-staggered grid (not mass-point!)
! Hence subtract .5 from each to get mass-point coords.
do n=1,prt_max
pnx = obs(n)
if(pnx.ne.-999) then
write(msg,'(2x,i7,3f8.3,2f9.3,2x,f9.3,2x,f9.3,3x,f6.2)') &
pnx,rio(pnx)-.5,rjo(pnx)-.5,rko(pnx),lat(n),lon(n), &
mlat(n),mlon(n),timeob(pnx)
call wrf_message(msg)
endif
enddo
if(obs(1).ne.-999) call wrf_message("")
endif
END SUBROUTINE print_obs_info
REAL FUNCTION ht_to_p (docs) ( h, pbbc, ppbc, phb, ph, g, ic, jc, dx, dy, & 2,1
k_start, k_end, kds,kde, ims,ime, jms,jme, kms,kme )
!*************************************************************************
! Purpose: Interpolate pressure at a specific x, y, height coordinate.
! The input pressure column (base and perturbation) is already
! horizontally interpolated to the x, y position. The input
! geopotential field (base: phb, pert: ph) is 3d, and must be
! horizontally interpolated (to iob, job) to get a height column (z_at_p)
! before p is vertically interpolated.
!*************************************************************************
IMPLICIT NONE
REAL, INTENT(IN) :: h ! height value (m)
INTEGER, INTENT(IN) :: k_start, k_end ! loop bounds
INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
REAL, INTENT(IN) :: pbbc(kds:kde) ! column base pressure (cb)
REAL, INTENT(IN) :: ppbc(kds:kde) ! column pressure perturbation (cb)
REAL, INTENT(IN) :: phb( ims:ime, kms:kme, jms:jme ) ! geopotential (base)
REAL, INTENT(IN) :: ph( ims:ime, kms:kme, jms:jme ) ! geopotential (perturbation)
REAL, INTENT(IN) :: g ! gravity constant
INTEGER, INTENT(IN) :: ic ! i-coord of desired p
INTEGER, INTENT(IN) :: jc ! j-coord of desired p
REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
! Local variables
INTEGER :: k ! loop counter
INTEGER :: klo ! lower k bound
REAL :: zlo ! lower z bound for h
REAL :: zhi ! upper z bound for h
REAL :: p ! interpolated pressure value
REAL :: ln_p ! log p
REAL :: ln_plo ! log plo
REAL :: ln_phi ! log phi
REAL :: z_at_p( kms:kme ) ! height at p levels
! Calculate z at p levels.
call get_height_column(phb, ph, g, ic, jc, dx, dy, z_at_p, &
k_start, k_end, kds,kde, &
ims,ime, jms,jme, kms,kme )
! Now we have pbbc, ppbc, z_at_p, so compute p at h. First, find
! bounding layers klo and khi so that z_at_p(klo) <= h <= z_at_p(khi)
ZLEVS: do k = k_start+1, k_end
klo = k-1
if(h .le. z_at_p(k)) then
EXIT ZLEVS
endif
enddo ZLEVS
zlo = z_at_p(klo)
zhi = z_at_p(klo+1)
! Interpolate natural log of pressure
ln_plo = log( pbbc(klo+1) + ppbc(klo+1) )
ln_phi = log( pbbc(klo) + ppbc(klo) )
if(h.le.zlo) then
ln_p = ln_phi ! set to k=1 pressure
else if (h.ge.zhi) then
ln_p = ln_plo ! set to k=k_end pressure
else
ln_p = ln_plo + (ln_phi-ln_plo)*((zhi-h)/(zhi-zlo))
endif
! Return pressure
p = exp(ln_p)
ht_to_p = p
RETURN
END FUNCTION ht_to_p
SUBROUTINE get_height_column (docs) ( phb, ph, g, ic, jc, dx, dy, z_at_p, & 1
k_start, k_end, kds,kde, &
ims,ime, jms,jme, kms,kme )
!*************************************************************************
! Purpose: Compute column height at ic, jc location on p levels
!*************************************************************************
IMPLICIT NONE
INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
REAL, INTENT(IN) :: phb( ims:ime, kms:kme, jms:jme ) ! geopotential (base)
REAL, INTENT(IN) :: ph( ims:ime, kms:kme, jms:jme ) ! geopotential (perturbation)
REAL, INTENT(IN) :: g ! gravity constant
INTEGER, INTENT(IN) :: ic ! i-coord of desired p
INTEGER, INTENT(IN) :: jc ! j-coord of desired p
REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
! Local variables
INTEGER :: k ! loop counter
REAL :: phbc(kds:kde) ! column geopotential (base)
REAL :: phc(kds:kde) ! column geopotential (perturbation)
REAL :: z_at_w(kds:kde) ! column height at w levels
do k = kds, kde
phbc(k) = &
(1.-DY)*( (1.-DX)*phb(IC,K,JC) + &
DX *phb(IC+1,K,JC) ) + &
DY* ( (1.-DX)*phb(IC,K,JC+1) + &
DX *phb(IC+1,K,JC+1) )
phc(k) = &
(1.-DY)*( (1.-DX)*ph(IC,K,JC) + &
DX *ph(IC+1,K,JC) ) + &
DY* ( (1.-DX)*ph(IC,K,JC+1) + &
DX *ph(IC+1,K,JC+1) )
z_at_w(k) = (phbc(k)+phc(k))/g
enddo
do k = k_start, k_end
z_at_p(k) = 0.5*(z_at_w(k) +z_at_w(k+1))
enddo
END SUBROUTINE get_height_column
#endif
END MODULE module_fddaobs_rtfdda