Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_fddaobs_rtfdda.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 89.6 KB

Rev	Line
[2759]	1	!WRF:MODEL_LAYER:PHYSICS
	2	!
	3	MODULE module_fddaobs_rtfdda
	4
	5	! This obs-nudging FDDA module (RTFDDA) is developed by the
	6	! NCAR/RAL/NSAP (National Security Application Programs), under the
	7	! sponsorship of ATEC (Army Test and Evaluation Commands). ATEC is
	8	! acknowledged for releasing this capability for WRF community
	9	! research applications.
	10	!
	11	! The NCAR/RAL RTFDDA module was adapted, and significantly modified
	12	! from the obs-nudging module in the standard MM5V3.1 which was originally
	13	! developed by PSU (Stauffer and Seaman, 1994).
	14	!
	15	! Yubao Liu (NCAR/RAL): lead developer of the RTFDDA module
	16	! Al Bourgeois (NCAR/RAL): lead engineer implementing RTFDDA into WRF-ARW
	17	! Nov. 2006
	18	!
	19	! References:
	20	!
	21	! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and J. Hacker, 2005: An
	22	! implementation of obs-nudging-based FDDA into WRF for supporting
	23	! ATEC test operations. 2005 WRF user workshop. Paper 10.7.
	24	!
	25	! Liu, Y., A. Bourgeois, T. Warner, S. Swerdlin and W. Yu, 2006: An update
	26	! on "obs-nudging"-based FDDA for WRF-ARW: Verification using OSSE
	27	! and performance of real-time forecasts. 2006 WRF user workshop. Paper 4.7.
	28
	29	!
	30	! Stauffer, D.R., and N.L. Seaman, 1994: Multi-scale four-dimensional data
	31	! assimilation. J. Appl. Meteor., 33, 416-434.
	32	!
	33	! http://www.rap.ucar.edu/projects/armyrange/references.html
	34	!
	35	! Modification History:
	36	! 03212007 Modified fddaobs_init to compute Lambert cone factor. -Al Bourgeois
	37
	38	CONTAINS
	39
	40	!------------------------------------------------------------------------------
	41	SUBROUTINE fddaobs_init(obs_nudge_opt, maxdom, inest, parid, &
	42	idynin, dtramp, fdaend, restart, &
	43	obs_twindo_cg, obs_twindo, itimestep, &
	44	cen_lat, cen_lon, stand_lon, &
	45	true_lat1, true_lat2, map_proj, &
	46	xlat, xlong, &
	47	e_sn, s_sn_cg, e_sn_cg, s_we_cg, e_we_cg, &
	48	#if ( EM_CORE == 1 )
	49	fdob, &
	50	#endif
	51	iprt, &
	52	ids,ide, jds,jde, kds,kde, &
	53	ims,ime, jms,jme, kms,kme, &
	54	its,ite, jts,jte, kts,kte)
	55	!-----------------------------------------------------------------------
	56	! This routine does initialization for real time fdda obs-nudging.
	57	!
	58	!-----------------------------------------------------------------------
	59	USE module_domain
	60	USE module_dm ! for wrf_dm_min_real
	61	!-----------------------------------------------------------------------
	62	IMPLICIT NONE
	63	!-----------------------------------------------------------------------
	64
	65	!=======================================================================
	66	! Definitions
	67	!-----------
	68	INTEGER, intent(in) :: maxdom
	69	INTEGER, intent(in) :: obs_nudge_opt(maxdom)
	70	INTEGER, intent(in) :: ids,ide, jds,jde, kds,kde, &
	71	ims,ime, jms,jme, kms,kme, &
	72	its,ite, jts,jte, kts,kte
	73	INTEGER, intent(in) :: inest
	74	INTEGER, intent(in) :: parid(maxdom)
	75	INTEGER, intent(in) :: idynin ! flag for dynamic initialization
	76	REAL, intent(in) :: dtramp ! time period for idynin ramping (min)
	77	REAL, intent(in) :: fdaend(maxdom) ! nudging end time for domain (min)
	78	LOGICAL, intent(in) :: restart
	79	REAL, intent(in) :: obs_twindo_cg ! time window on coarse grid
	80	REAL, intent(in) :: obs_twindo
	81	INTEGER, intent(in) :: itimestep
	82	REAL , INTENT(IN) :: cen_lat ! domain center latitude (for map proj)
	83	REAL , INTENT(IN) :: cen_lon ! domain center longitude (for map proj)
	84	REAL , INTENT(IN) :: stand_lon ! domain standard longitude (for map proj)
	85	REAL, intent(in) :: true_lat1 ! domain truelat1 (for map proj)
	86	REAL, intent(in) :: true_lat2 ! domain truelat2 (for map proj)
	87	INTEGER, intent(in) :: map_proj ! map projection index
	88	REAL, DIMENSION( ims:ime, jms:jme ), &
	89	INTENT(IN) :: xlat, xlong ! lat/long locations on mass point grid
	90	INTEGER, intent(in) :: e_sn ! ending south-north grid index
	91	INTEGER, intent(in) :: s_sn_cg ! starting south-north coarse-grid index
	92	INTEGER, intent(in) :: e_sn_cg ! ending south-north coarse-grid index
	93	INTEGER, intent(in) :: s_we_cg ! starting west-east coarse-grid index
	94	INTEGER, intent(in) :: e_we_cg ! ending west-east coarse-grid index
	95	#if ( EM_CORE == 1 )
	96	TYPE(fdob_type), intent(inout) :: fdob
	97	#endif
	98	LOGICAL, intent(in) :: iprt ! Flag enabling printing warning messages
	99
	100	! Local variables
	101	logical :: nudge_flag ! nudging flag for this nest
	102	integer :: ktau ! current timestep
	103	integer :: nest ! loop counter
	104	integer :: idom ! domain id
	105	integer :: parent ! parent domain
	106	real :: conv ! 180/pi
	107	real :: tl1 ! Lambert standard parallel 1
	108	real :: tl2 ! Lambert standard parallel 2
	109	real :: xn1
	110	real :: known_lat ! Latitude of domain point (i,j)=(1,1)
	111	real :: known_lon ! Longitude of domain point (i,j)=(1,1)
	112
	113	#if ( EM_CORE == 1 )
	114
	115	! Set flag for nudging on pressure (not sigma) surfaces.
	116	fdob%iwtsig = 0
	117
	118	!**************************************************************************
	119	! * Initialize datend for dynamic initialization (ajb added 08132008) *
	120	!**************************************************************************
	121	! Set ending nudging date (used with dynamic ramp-down) to zero.
	122	fdob%datend = 0.
	123	fdob%ieodi = 0
	124
	125	! Check for dynamic initialization flag
	126	if(idynin.eq.1)then
	127	! Set datend to time in minutes after which data are assumed to have ended.
	128	if(dtramp.gt.0.)then
	129	fdob%datend = fdaend(inest) - dtramp
	130	else
	131	fdob%datend = fdaend(inest)
	132	endif
	133	if(iprt) then
	134	print*,' '
	135	print,' DYNAMIC-INITIALIZATION OPTION FOR INEST = ',inest, ' *'
	136	print*,' FDAEND,DATEND,DTRAMP: ',fdaend(inest),fdob%datend,dtramp
	137	print*,' '
	138	endif
	139	endif
	140	! * end dynamic initialization section *
	141	!**************************************************************************
	142
	143	! Set time window.
	144	fdob%window = obs_twindo
	145
	146	if(inest.eq.1) then
	147	if(obs_twindo .eq. 0.) then
	148	if(iprt) then
	149	write(6,'(/a)') '*** WARNING: TWINDO=0 on the coarse domain.'
	150	write(6,'(a/)') '*** Did you forget to set twindo in the fdda namelist?'
	151	endif
	152	endif
	153	else ! nest
	154	if(obs_twindo .eq. 0.) then
	155	fdob%window = obs_twindo_cg
	156	if(iprt) then
	157	write(6,'(/a,i2)') 'WARNING: TWINDO=0. for nest ',inest
	158	write(6,'(a,f12.5,a/)') 'Default to coarse-grid value of ', obs_twindo_cg,' hrs'
	159	endif
	160	endif
	161	endif
	162
	163	! Initialize flags.
	164
	165	fdob%domain_tot=0
	166	do nest=1,maxdom
	167	fdob%domain_tot = fdob%domain_tot + obs_nudge_opt(nest)
	168	end do
	169
	170	! Set parameters.
	171
	172	fdob%pfree = 50.0
	173	fdob%rinfmn = 1.0
	174	fdob%rinfmx = 2.0
	175	fdob%dpsmx = 7.5
	176	fdob%dcon = 1.0/fdob%dpsmx
	177
	178	! Get known lat and lon and store these on all processors in fdob structure, for
	179	! later use in projection x-formation to map (lat,lon) to (x,y) for each obs.
	180	IF (its .eq. 1 .AND. jts .eq. 1) then
	181	known_lat = xlat(1,1)
	182	known_lon = xlong(1,1)
	183	ELSE
	184	known_lat = 9999.
	185	known_lon = 9999.
	186	END IF
	187	fdob%known_lat = wrf_dm_min_real(known_lat)
	188	fdob%known_lon = wrf_dm_min_real(known_lon)
	189
	190	fdob%sn_maxcg = e_sn_cg - s_sn_cg + 1 ! coarse domain grid dimension in N-S
	191	fdob%we_maxcg = e_we_cg - s_we_cg + 1 ! coarse domain grid dimension in W-E
	192	fdob%sn_end = e_sn - 1 ! ending S-N grid coordinate
	193
	194	! Calculate the nest levels, levidn. Note that each nest
	195	! must know the nest levels levidn(maxdom) of each domain.
	196	do nest=1,maxdom
	197
	198	! Initialize nest level for each domain.
	199	if (nest .eq. 1) then
	200	fdob%levidn(nest) = 0 ! Mother domain has nest level 0
	201	else
	202	fdob%levidn(nest) = 1 ! All other domains start with 1
	203	endif
	204	idom = nest
	205	100 parent = parid(idom) ! Go up the parent tree
	206
	207	if (parent .gt. 1) then ! If not up to mother domain
	208	fdob%levidn(nest) = fdob%levidn(nest) + 1
	209	idom = parent
	210	goto 100
	211	endif
	212	enddo
	213
	214	! Check to see if the nudging flag has been set. If not,
	215	! simply RETURN.
	216	nudge_flag = (obs_nudge_opt(inest) .eq. 1)
	217	if (.not. nudge_flag) return
	218
	219	ktau = itimestep
	220	if(restart) then
	221	fdob%ktaur = ktau
	222	else
	223	fdob%ktaur = 0
	224	endif
	225
	226	RETURN
	227	#endif
	228	END SUBROUTINE fddaobs_init
	229
	230	#if ( EM_CORE == 1 )
	231	!-----------------------------------------------------------------------
	232	SUBROUTINE errob(inest, ub, vb, tb, t0, qvb, pbase, pp, rovcp, &
	233	uratx, vratx, tratx, nndgv, &
	234	nerrf, niobf, maxdom, levidn, parid, nstat, &
	235	nstaw, iswind, &
	236	istemp, ismois, ispstr, rio, rjo, rko, varobs, &
	237	errf, i_parent_start, j_parent_start, &
	238	ktau, iratio, npfi, nobs_prt, iprt, &
	239	ids,ide, jds,jde, kds,kde, &
	240	ims,ime, jms,jme, kms,kme, &
	241	its,ite, jts,jte, kts,kte )
	242
	243	!-----------------------------------------------------------------------
	244	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	245	USE module_dm, ONLY : get_full_obs_vector
	246	#endif
	247
	248	!-----------------------------------------------------------------------
	249	IMPLICIT NONE
	250	!-----------------------------------------------------------------------
	251	!
	252	! PURPOSE: THIS SUBROUTINE CALCULATES THE DIFFERENCE BETWEEN THE
	253	! OBSERVED VALUES AND THE FORECAST VALUES AT THE OBSERVATION
	254	! POINTS. THE OBSERVED VALUES CLOSEST TO THE CURRENT
	255	! FORECAST TIME (XTIME) WERE DETERMINED IN SUBROUTINE
	256	! IN4DOB AND STORED IN ARRAY VAROBS. THE DIFFERENCES
	257	! CALCULATED BY SUBROUTINE ERROB WILL BE STORED IN ARRAY
	258	! ERRF FOR THE NSTA OBSERVATION LOCATIONS. MISSING
	259	! OBSERVATIONS ARE DENOTED BY THE DUMMY VALUE 99999.
	260	!
	261	! HISTORY: Original author: MM5 version???
	262	! 02/04/2004 - Creation of WRF version. Al Bourgeois
	263	! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
	264	!------------------------------------------------------------------------------
	265
	266	! THE STORAGE ORDER IN VAROBS AND ERRF IS AS FOLLOWS:
	267	! IVAR VARIABLE TYPE(TAU-1)
	268	! ---- --------------------
	269	! 1 U error
	270	! 2 V error
	271	! 3 T error
	272	! 4 Q error
	273	! 5 Surface press error at T points (not used)
	274	! 6 Model surface press at T-points
	275	! 7 Model surface press at U-points
	276	! 8 Model surface press at V-points
	277	! 9 RKO at U-points
	278
	279	!-----------------------------------------------------------------------
	280	!
	281	! Description of input arguments.
	282	!
	283	!-----------------------------------------------------------------------
	284
	285	INTEGER, INTENT(IN) :: inest ! Domain index.
	286	INTEGER, INTENT(IN) :: nndgv ! Number of nudge variables.
	287	INTEGER, INTENT(IN) :: nerrf ! Number of error fields.
	288	INTEGER, INTENT(IN) :: niobf ! Number of observations.
	289	INTEGER, INTENT(IN) :: maxdom ! Maximum number of domains.
	290	INTEGER, INTENT(IN) :: levidn(maxdom) ! Level of nest.
	291	INTEGER, INTENT(IN) :: parid(maxdom) ! Id of parent grid.
	292	INTEGER, INTENT(IN) :: i_parent_start(maxdom) ! Start i index in parent domain.
	293	INTEGER, INTENT(IN) :: j_parent_start(maxdom) ! Start j index in parent domain.
	294	INTEGER, INTENT(IN) :: ktau
	295	INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
	296	INTEGER, INTENT(IN) :: npfi ! Coarse-grid diagnostics freq.
	297	INTEGER, INTENT(IN) :: nobs_prt ! Number of current obs to print info.
	298	LOGICAL, INTENT(IN) :: iprt ! Print flag
	299	INTEGER, INTENT(IN) :: nstat ! # stations held for use
	300	INTEGER, INTENT(IN) :: nstaw ! # stations in current window
	301	INTEGER, intent(in) :: iswind
	302	INTEGER, intent(in) :: istemp
	303	INTEGER, intent(in) :: ismois
	304	INTEGER, intent(in) :: ispstr
	305	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
	306	INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
	307	INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
	308
	309	REAL, INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
	310	REAL, INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
	311	REAL, INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
	312	REAL, INTENT(IN) :: t0
	313	REAL, INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
	314	REAL, INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme )
	315	REAL, INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme ) ! Press. perturbation (Pa)
	316	REAL, INTENT(IN) :: rovcp
	317	REAL, INTENT(IN) :: uratx( ims:ime, jms:jme ) ! U to U10 ratio on mass points.
	318	REAL, INTENT(IN) :: vratx( ims:ime, jms:jme ) ! V to V10 ratio on mass points.
	319	REAL, INTENT(IN) :: tratx( ims:ime, jms:jme ) ! T to TH2 ratio on mass points.
	320	REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
	321	REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
	322	REAL, INTENT(INOUT) :: rko(niobf)
	323	REAL, INTENT(INOUT) :: varobs(nndgv, niobf)
	324	REAL, INTENT(INOUT) :: errf(nerrf, niobf)
	325
	326	! Local variables
	327	INTEGER :: iobmg(niobf) ! Obs i-coord on mass grid
	328	INTEGER :: jobmg(niobf) ! Obs j-coord on mass grid
	329	INTEGER :: ia(niobf)
	330	INTEGER :: ib(niobf)
	331	INTEGER :: ic(niobf)
	332	REAL :: pbbo(kds:kde) ! column base pressure (cb) at obs loc.
	333	REAL :: ppbo(kds:kde) ! column pressure perturbation (cb) at obs loc.
	334
	335	REAL :: ra(niobf)
	336	REAL :: rb(niobf)
	337	REAL :: rc(niobf)
	338	REAL :: dxobmg(niobf) ! Interp. fraction (x dir) referenced to mass-grid
	339	REAL :: dyobmg(niobf) ! Interp. fraction (y dir) referenced to mass-grid
	340	INTEGER MM(MAXDOM)
	341	INTEGER NNL
	342	real :: uratio( ims:ime, jms:jme ) ! U to U10 ratio on momentum points.
	343	real :: vratio( ims:ime, jms:jme ) ! V to V10 ratio on momentum points.
	344	real :: pug1,pug2,pvg1,pvg2
	345
	346	! Define staggers for U, V, and T grids, referenced from non-staggered grid.
	347	real, parameter :: gridx_t = 0.5 ! Mass-point x stagger
	348	real, parameter :: gridy_t = 0.5 ! Mass-point y stagger
	349	real, parameter :: gridx_u = 0.0 ! U-point x stagger
	350	real, parameter :: gridy_u = 0.5 ! U-point y stagger
	351	real, parameter :: gridx_v = 0.5 ! V-point x stagger
	352	real, parameter :: gridy_v = 0.0 ! V-point y stagger
	353
	354	real :: dummy = 99999.
	355
	356	real :: pbhi, pphi
	357	real :: press,ttemp !ajb scratch variables
	358	! real model_temp,pot_temp !ajb scratch variables
	359
	360	!*** DECLARATIONS FOR IMPLICIT NONE
	361	integer nsta,ivar,n,ityp
	362	integer iob,job,kob,iob_ms,job_ms
	363	integer k,kbot,nml,nlb,nle
	364	integer iobm,jobm,iobp,jobp,kobp,inpf,i,j
	365	integer i_start,i_end,j_start,j_end ! loop ranges for uratio,vratio calc.
	366	integer k_start,k_end
	367	integer ips ! For printing obs information
	368
	369	real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
	370	real pob
	371	real grfacx,grfacy,uratiob,vratiob,tratiob,tratxob,fnpf
	372	real stagx ! For x correction to mass-point stagger
	373	real stagy ! For y correction to mass-point stagger
	374
	375	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	376	LOGICAL MP_LOCAL_DUMMASK(NIOBF) ! Mask for work to be done on this processor
	377	LOGICAL MP_LOCAL_UOBMASK(NIOBF) ! Dot-point mask
	378	LOGICAL MP_LOCAL_VOBMASK(NIOBF) ! Dot-point mask
	379	LOGICAL MP_LOCAL_COBMASK(NIOBF) ! Cross-point mask
	380	#endif
	381	! LOGICAL, EXTERNAL :: TILE_MASK
	382
	383	NSTA=NSTAT
	384
	385	! FIRST, DETERMINE THE GRID TYPE CORRECTION FOR U-momentum, V-momentum,
	386	! AND MASS POINTS, AND WHEN INEST=2, CONVERT THE STORED COARSE MESH INDICES
	387	! TO THE FINE MESH INDEX EQUIVALENTS
	388
	389	! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
	390
	391	if (iprt) then
	392	write(6,'(a,i5,a,i2,a,i5,a)') '++++++CALL ERROB AT KTAU = ', &
	393	KTAU,' AND INEST = ',INEST,': NSTA = ',NSTAW,' ++++++'
	394	endif
	395
	396	ERRF = 0.0 ! Zero out errf array
	397
	398	! Set up loop bounds for this grid's boundary conditions
	399	i_start = max( its-1,ids )
	400	i_end = min( ite+1,ide-1 )
	401	j_start = max( jts-1,jds )
	402	j_end = min( jte+1,jde-1 )
	403	k_start = kts
	404	k_end = min( kte, kde-1 )
	405
	406	DO ITYP=1,3 ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP
	407
	408	! Set grid stagger
	409	IF(ITYP.EQ.1) THEN ! U-POINT CASE
	410	GRIDX = gridx_u
	411	GRIDY = gridy_u
	412	ELSE IF(ITYP.EQ.2) THEN ! V-POINT CASE
	413	GRIDX = gridx_v
	414	GRIDY = gridy_v
	415	ELSE ! MASS-POINT CASE
	416	GRIDX = gridx_t
	417	GRIDY = gridy_t
	418	ENDIF
	419
	420	! Compute URATIO and VRATIO fields on momentum (u,v) points.
	421	IF(ityp.eq.1)THEN
	422	call upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, uratx, uratio)
	423	ELSE IF (ityp.eq.2) THEN
	424	call vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, vratx, vratio)
	425	ENDIF
	426
	427	IF(INEST.EQ.1) THEN ! COARSE MESH CASE...
	428	DO N=1,NSTA
	429	RA(N)=RIO(N)-GRIDX
	430	RB(N)=RJO(N)-GRIDY
	431	IA(N)=RA(N)
	432	IB(N)=RB(N)
	433	IOB=MAX0(1,IA(N))
	434	IOB=MIN0(IOB,ide-1)
	435	JOB=MAX0(1,IB(N))
	436	JOB=MIN0(JOB,jde-1)
	437	DXOB=RA(N)-FLOAT(IA(N))
	438	DYOB=RB(N)-FLOAT(IB(N))
	439
	440	! Save mass-point arrays for computing rko for all var types
	441	if(ityp.eq.1) then
	442	iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
	443	jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
	444	dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
	445	dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
	446	endif
	447	iob_ms = iobmg(n)
	448	job_ms = jobmg(n)
	449	dxob_ms = dxobmg(n)
	450	dyob_ms = dyobmg(n)
	451
	452
	453	!if(n.eq.1 .and. iprt) then
	454	! write(6,*) 'ERROB - COARSE MESH:'
	455	! write(6,'(a,i1,a,i1,4(a,f5.2),2(a,i3),2(a,f6.3))') 'OBS= ',n, &
	456	! ' ityp= ',ityp, &
	457	! ' ra= ',ra(n),' rb= ',rb(n), &
	458	! ' rio= ',rio(n),' rjo= ',rjo(n), &
	459	! ' iob= ',iob,' job= ',job, &
	460	! ' dxob= ',dxob,' dyob= ',dyob
	461	! write(6,'(a,i3,a,i3,a,f5.2,a,f5.2)') &
	462	! ' iob_ms= ',iob_ms,' job_ms= ',job_ms, &
	463	! ' dxob_ms= ',dxob_ms,' dyob_ms= ',dyob_ms
	464	!endif
	465
	466	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	467	! Set mask for obs to be handled by this processor
	468	MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
	469
	470	IF ( MP_LOCAL_DUMMASK(N) ) THEN
	471	#endif
	472
	473	! Interpolate pressure to obs location column and convert from Pa to cb.
	474
	475	do k = kds, kde
	476	pbbo(k) = .001*( &
	477	(1.-DYOB_MS)( (1.-DXOB_MS)pbase(IOB_MS,K,JOB_MS) + &
	478	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
	479	DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
	480	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
	481	ppbo(k) = .001*( &
	482	(1.-DYOB_MS)( (1.-DXOB_MS)pp(IOB_MS,K,JOB_MS) + &
	483	DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
	484	DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
	485	DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
	486
	487	! write(6,'(a,i2,2(a,f9.3)') ' k= ',k,' pbbo= ',pbbo(k),' ppbo= ',ppbo(k)
	488	enddo
	489
	490	!ajb 20040119: Note based on bugfix for dot/cross points split across processors,
	491	!ajb which was actually a serial code fix: The ityp=2 (v-points) and
	492	!ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
	493	!ajb case rko! This is necessary for bit-for-bit reproducability
	494	!ajb with the parallel run. (coarse mesh)
	495
	496
	497	if(abs(rko(n)+99).lt.1.)then
	498	pob = varobs(5,n)
	499
	500	if(pob .gt.-800000.)then
	501	do k=k_end-1,1,-1
	502	kbot = k
	503	if(pob .le. pbbo(k)+ppbo(k)) then
	504	goto 199
	505	endif
	506	enddo
	507	199 continue
	508
	509	pphi = ppbo(kbot+1)
	510	pbhi = pbbo(kbot+1)
	511
	512	rko(n) = real(kbot+1)- &
	513	( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
	514
	515	rko(n)=max(rko(n),1.0)
	516	endif
	517	endif
	518
	519	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	520	ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
	521	#endif
	522
	523	RC(N)=RKO(N)
	524
	525	ENDDO ! END COARSE MESH LOOP OVER NSTA
	526
	527	ELSE ! FINE MESH CASE
	528
	529	!**********************************************************************
	530	!ajb_07012008: Conversion of obs coordinates to the fine mesh here
	531	!ajb is no longer necessary, since implementation of the WRF map pro-
	532	!ajb jections (to each nest directly) is implemented in sub in4dob.
	533	!**********************************************************************
	534	!ajb
	535	!ajb GET (I,J,K) OF OBSERVATIONS ON FINE MESH VALUES.
	536	DO N=1,NSTA
	537
	538	! write(6,*) 'UPDATE: ra(n) = ',ra(n),' rb(n) = ',rb(n)
	539
	540	RA(N)=RIO(N)-GRIDX ! ajb added 07012008
	541	RB(N)=RJO(N)-GRIDY ! ajb added 07012008
	542	IA(N)=RA(N)
	543	IB(N)=RB(N)
	544	IOB=MAX0(1,IA(N))
	545	IOB=MIN0(IOB,ide-1)
	546	JOB=MAX0(1,IB(N))
	547	JOB=MIN0(JOB,jde-1)
	548	DXOB=RA(N)-FLOAT(IA(N))
	549	DYOB=RB(N)-FLOAT(IB(N))
	550
	551	! Save mass-point arrays for computing rko for all var types
	552	if(ityp.eq.1) then
	553	stagx = grfacx - gridx_t !Correct x stagger to mass-point
	554	stagy = grfacy - gridy_t !Correct y stagger to mass-point
	555	iobmg(n) = MIN0(MAX0(1,int(RA(n)+stagx)),ide-1)
	556	jobmg(n) = MIN0(MAX0(1,int(RB(n)+stagy)),jde-1)
	557	dxobmg(n) = RA(N)+stagx-FLOAT(int(RA(N)+stagx))
	558	dyobmg(n) = RB(N)+stagy-FLOAT(int(RB(N)+stagy))
	559	endif
	560	iob_ms = iobmg(n)
	561	job_ms = jobmg(n)
	562	dxob_ms = dxobmg(n)
	563	dyob_ms = dyobmg(n)
	564
	565	!if(n.eq.1) then
	566	! write(6,*) 'ERROB - FINE MESH:'
	567	! write(6,*) 'RA = ',ra(n),' RB = ',rb(n)
	568	! write(6,'(a,i1,a,i1,4(a,f5.2),2(a,i3),2(a,f6.3))') 'OBS= ',n, &
	569	! ' ityp= ',ityp, &
	570	! ' ra= ',ra(n),' rb= ',rb(n), &
	571	! ' rio= ',rio(n),' rjo= ',rjo(n), &
	572	! ' iob= ',iob,' job= ',job, &
	573	! ' dxob= ',dxob,' dyob= ',dyob
	574	! write(6,'(a,i3,a,i3,a,f5.2,a,f5.2)') &
	575	! ' iob_ms= ',iob_ms,' job_ms= ',job_ms, &
	576	! ' dxob_ms= ',dxob_ms,' dyob_ms= ',dyob_ms
	577	!endif
	578
	579	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	580	! Set mask for obs to be handled by this processor
	581	MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
	582
	583	IF ( MP_LOCAL_DUMMASK(N) ) THEN
	584	#endif
	585
	586	! Interpolate pressure to obs location column and convert from Pa to cb.
	587
	588	do k = kds, kde
	589	pbbo(k) = .001*( &
	590	(1.-DYOB_MS)( (1.-DXOB_MS)pbase(IOB_MS,K,JOB_MS) + &
	591	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
	592	DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
	593	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
	594	ppbo(k) = .001*( &
	595	(1.-DYOB_MS)( (1.-DXOB_MS)pp(IOB_MS,K,JOB_MS) + &
	596	DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
	597	DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
	598	DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
	599
	600	! write(6,'(a,i2,2(a,f9.3)') ' k= ',k,' pbbo= ',pbbo(k),' ppbo= ',ppbo(k)
	601	enddo
	602
	603	!ajb 20040119: Note based on bugfix for dot/cross points split across processors,
	604	!ajb which was actually a serial code fix: The ityp=2 (v-points) and
	605	!ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
	606	!ajb case) rko! This is necessary for bit-for-bit reproducability
	607	!ajb with parallel run. (fine mesh)
	608
	609	if(abs(rko(n)+99).lt.1.)then
	610	pob = varobs(5,n)
	611
	612	if(pob .gt.-800000.)then
	613	do k=k_end-1,1,-1
	614	kbot = k
	615	if(pob .le. pbbo(k)+ppbo(k)) then
	616	goto 198
	617	endif
	618	enddo
	619	198 continue
	620
	621	pphi = ppbo(kbot+1)
	622	pbhi = pbbo(kbot+1)
	623
	624	rko(n) = real(kbot+1)- &
	625	( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
	626	rko(n)=max(rko(n),1.0)
	627	endif
	628	endif
	629
	630	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	631	ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
	632	#endif
	633
	634	RC(N)=RKO(N)
	635
	636	ENDDO ! END FINE MESH LOOP OVER NSTA
	637
	638	ENDIF ! end if(inest.eq.1)
	639
	640	! Print obs information.
	641	if (iprt) then
	642	if(ityp.eq.3) then !mass-point case
	643	ips = min(nstaw,nobs_prt)
	644	if(ips.gt.0) then
	645	write(6,*) ' OBS# I J K'
	646	endif
	647	do n=1,ips
	648	write(6,'(2x,i4,3f8.3)') n,ra(n),rb(n),rc(n)
	649	enddo
	650	endif
	651	endif
	652
	653	! INITIALIZE THE ARRAY OF DIFFERENCES BETWEEN THE OBSERVATIONS
	654	! AND THE LOCAL FORECAST VALUES TO ZERO. FOR THE FINE MESH
	655	! ONLY, SET THE DIFFERENCE TO A LARGE DUMMY VALUE IF THE
	656	! OBSERVATION IS OUTSIDE THE FINE MESH DOMAIN.
	657
	658	! SET DIFFERENCE VALUE TO A DUMMY VALUE FOR OBS POINTS OUTSIDE
	659	! CURRENT DOMAIN
	660	IF(ITYP.EQ.1) THEN
	661	NLB=1
	662	NLE=1
	663	ELSE IF(ITYP.EQ.2) THEN
	664	NLB=2
	665	NLE=2
	666	ELSE
	667	NLB=3
	668	NLE=5
	669	ENDIF
	670	DO IVAR=NLB,NLE
	671	DO N=1,NSTA
	672	IF((RA(N)-1.).LT.0)THEN
	673	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
	674	ENDIF
	675	IF((RB(N)-1.).LT.0)THEN
	676	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
	677	ENDIF
	678	IF((FLOAT(ide)-2.0*GRIDX-RA(N)-1.E-10).LT.0)THEN
	679	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
	680	ENDIF
	681	IF((FLOAT(jde)-2.0*GRIDY-RB(N)-1.E-10).LT.0)THEN
	682	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
	683	ENDIF
	684	if(rc(n).lt.1.)errf(ivar,n)=errf(ivar,n)+dummy
	685	ENDDO
	686	ENDDO
	687
	688	! NOW FIND THE EXACT OFFSET OF EACH OBSERVATION FROM THE
	689	! GRID POINT TOWARD THE LOWER LEFT
	690	DO N=1,NSTA
	691	IA(N)=RA(N)
	692	IB(N)=RB(N)
	693	IC(N)=RC(N)
	694	ENDDO
	695	DO N=1,NSTA
	696	RA(N)=RA(N)-FLOAT(IA(N))
	697	RB(N)=RB(N)-FLOAT(IB(N))
	698	RC(N)=RC(N)-FLOAT(IC(N))
	699	ENDDO
	700	! PERFORM A TRILINEAR EIGHT-POINT (3-D) INTERPOLATION
	701	! TO FIND THE FORECAST VALUE AT THE EXACT OBSERVATION
	702	! POINTS FOR U, V, T, AND Q.
	703
	704	! Compute local masks for dot and cross points.
	705	if(ityp.eq.1) then
	706	DO N=1,NSTA
	707	IOB=MAX0(1,IA(N))
	708	IOB=MIN0(IOB,ide-1)
	709	JOB=MAX0(1,IB(N))
	710	JOB=MIN0(JOB,jde-1)
	711	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	712	! Set mask for U-momemtum points to be handled by this processor
	713	MP_LOCAL_UOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
	714	#endif
	715	ENDDO
	716	endif
	717	if(ityp.eq.2) then
	718	DO N=1,NSTA
	719	IOB=MAX0(1,IA(N))
	720	IOB=MIN0(IOB,ide-1)
	721	JOB=MAX0(1,IB(N))
	722	JOB=MIN0(JOB,jde-1)
	723	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	724	! Set mask for V-momentum points to be handled by this processor
	725	MP_LOCAL_VOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
	726	#endif
	727	ENDDO
	728	endif
	729	if(ityp.eq.3) then
	730	DO N=1,NSTA
	731	IOB=MAX0(1,IA(N))
	732	IOB=MIN0(IOB,ide-1)
	733	JOB=MAX0(1,IB(N))
	734	JOB=MIN0(JOB,jde-1)
	735	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	736	! Set mask for cross (mass) points to be handled by this processor
	737	MP_LOCAL_COBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
	738	#endif
	739	ENDDO
	740	endif
	741
	742	!**********************************************************
	743	! PROCESS U VARIABLE (IVAR=1)
	744	!**********************************************************
	745	IF(ITYP.EQ.1) THEN
	746	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	747	DO N=1,NSTA
	748	IF(MP_LOCAL_UOBMASK(N)) THEN
	749	ERRF(9,N)=rko(n) !RKO is needed by neighboring processors !ajb
	750	ENDIF
	751	ENDDO
	752	#endif
	753	IF(ISWIND.EQ.1) THEN
	754	DO N=1,NSTA
	755	IOB=MAX0(2,IA(N))
	756	IOB=MIN0(IOB,ide-1)
	757	IOBM=MAX0(1,IOB-1)
	758	IOBP=MIN0(ide-1,IOB+1)
	759	JOB=MAX0(2,IB(N))
	760	JOB=MIN0(JOB,jde-1)
	761	JOBM=MAX0(1,JOB-1)
	762	JOBP=MIN0(jde-1,JOB+1)
	763	KOB=MIN0(K_END,IC(N))
	764
	765	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	766	IF(MP_LOCAL_UOBMASK(N))THEN ! Do if obs on this processor
	767	#endif
	768	KOBP=MIN0(KOB+1,k_end)
	769	DXOB=RA(N)
	770	DYOB=RB(N)
	771	DZOB=RC(N)
	772
	773	! Compute surface pressure values at surrounding U and V points
	774	PUG1 = .5*( pbase(IOBM,1,JOB) + pbase(IOB,1,JOB) )
	775	PUG2 = .5*( pbase(IOB,1,JOB) + pbase(IOBP,1,JOB) )
	776
	777	! This is to correct obs to model sigma level using reverse similarity theory
	778	if(rko(n).eq.1.0)then
	779	uratiob=((1.-DYOB)((1.-DXOB)uratio(IOB,JOB)+ &
	780	DXOB*uratio(IOBP,JOB) &
	781	)+DYOB((1.-DXOB)uratio(IOB,JOBP)+ &
	782	DXOB*uratio(IOBP,JOBP)))
	783	else
	784	uratiob=1.
	785	endif
	786	!YLIU Some PBL scheme do not define the vratio/uratio
	787	if(abs(uratiob).lt.1.0e-3) then
	788	uratiob=1.
	789	endif
	790
	791	! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
	792	! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
	793	! OUTSIDE THE CURRENT DOMAIN
	794
	795	! U COMPONENT WIND ERROR
	796	ERRF(1,N)=ERRF(1,N)+uratiobVAROBS(1,N)-((1.-DZOB) &
	797	((1.-DyOB)*((1.- &
	798	DxOB)UB(IOB,KOB,JOB)+DxOBUB(IOB+1,KOB,JOB) &
	799	)+DyOB((1.-DxOB)UB(IOB,KOB,JOB+1)+DxOB* &
	800	UB(IOB+1,KOB,JOB+1)))+DZOB((1.-DyOB)((1.-DxOB) &
	801	UB(IOB,KOBP,JOB)+DxOBUB(IOB+1,KOBP,JOB))+ &
	802	DyOB((1.-DxOB)UB(IOB,KOBP,JOB+1)+DxOB* &
	803	UB(IOB+1,KOBP,JOB+1))))
	804
	805	! if(n.le.10) then
	806	! write(6,*)
	807	! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF1 at ',iob,job,kob, &
	808	! ' N = ',n,' inest = ',inest
	809	! write(6,*) 'VAROBS(1,N) = ',varobs(1,n)
	810	! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
	811	! write(6,*) 'UB(IOB,KOB,JOB) = ',UB(IOB,KOB,JOB)
	812	! write(6,*) 'UB(IOB+1,KOB,JOB) = ',UB(IOB+1,KOB,JOB)
	813	! write(6,*) 'UB(IOB,KOB,JOB+1) = ',UB(IOB,KOB,JOB+1)
	814	! write(6,*) 'UB(IOB+1,KOB,JOB+1) = ',UB(IOB+1,KOB,JOB+1)
	815	! write(6,*) 'UB(IOB,KOBP,JOB) = ',UB(IOB,KOBP,JOB)
	816	! write(6,*) 'UB(IOB+1,KOBP,JOB) = ',UB(IOB+1,KOBP,JOB)
	817	! write(6,*) 'UB(IOB,KOBP,JOB+1) = ',UB(IOB,KOBP,JOB+1)
	818	! write(6,*) 'UB(IOB+1,KOBP,JOB+1) = ',UB(IOB+1,KOBP,JOB+1)
	819	! write(6,*) 'uratiob = ',uratiob
	820	! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
	821	! write(6,*) 'ERRF(1,N) = ',errf(1,n)
	822	! write(6,*)
	823	! endif
	824
	825
	826	! Store model surface pressure (not the error!) at U point.
	827	ERRF(7,N)=.001( (1.-DXOB)PUG1 + DXOB*PUG2 )
	828
	829	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	830	ENDIF ! end IF( MP_LOCAL_UOBMASK(N) )
	831	#endif
	832	ENDDO ! END U-POINT LOOP OVER OBS
	833
	834	ENDIF ! end if(iswind.eq.1)
	835
	836	ENDIF ! ITYP=1: PROCESS U
	837
	838	!**********************************************************
	839	! PROCESS V VARIABLE (IVAR=2)
	840	!**********************************************************
	841	IF(ITYP.EQ.2) THEN
	842
	843	IF(ISWIND.EQ.1) THEN
	844	DO N=1,NSTA
	845	IOB=MAX0(2,IA(N))
	846	IOB=MIN0(IOB,ide-1)
	847	IOBM=MAX0(1,IOB-1)
	848	IOBP=MIN0(ide-1,IOB+1)
	849	JOB=MAX0(2,IB(N))
	850	JOB=MIN0(JOB,jde-1)
	851	JOBM=MAX0(1,JOB-1)
	852	JOBP=MIN0(jde-1,JOB+1)
	853	KOB=MIN0(K_END,IC(N))
	854
	855	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	856	IF(MP_LOCAL_VOBMASK(N))THEN ! Do if obs on this processor
	857	#endif
	858	KOBP=MIN0(KOB+1,k_end)
	859	DXOB=RA(N)
	860	DYOB=RB(N)
	861	DZOB=RC(N)
	862
	863	! Compute surface pressure values at surrounding U and V points
	864	PVG1 = .5*( pbase(IOB,1,JOBM) + pbase(IOB,1,JOB) )
	865	PVG2 = .5*( pbase(IOB,1,JOB) + pbase(IOB,1,JOBP) )
	866
	867	! This is to correct obs to model sigma level using reverse similarity theory
	868	if(rko(n).eq.1.0)then
	869	vratiob=((1.-DYOB)((1.-DXOB)vratio(IOB,JOB)+ &
	870	DXOB*vratio(IOBP,JOB) &
	871	)+DYOB((1.-DXOB)vratio(IOB,JOBP)+ &
	872	DXOB*vratio(IOBP,JOBP)))
	873	else
	874	vratiob=1.
	875	endif
	876	!YLIU Some PBL scheme do not define the vratio/uratio
	877	if(abs(vratiob).lt.1.0e-3) then
	878	vratiob=1.
	879	endif
	880
	881	! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
	882	! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
	883	! OUTSIDE THE CURRENT DOMAIN
	884
	885	! V COMPONENT WIND ERROR
	886	ERRF(2,N)=ERRF(2,N)+vratiobVAROBS(2,N)-((1.-DZOB) &
	887	((1.-DyOB)*((1.- &
	888	DxOB)VB(IOB,KOB,JOB)+DxOBVB(IOB+1,KOB,JOB) &
	889	)+DyOB((1.-DxOB)VB(IOB,KOB,JOB+1)+DxOB* &
	890	VB(IOB+1,KOB,JOB+1)))+DZOB((1.-DyOB)((1.-DxOB) &
	891	VB(IOB,KOBP,JOB)+DxOBVB(IOB+1,KOBP,JOB))+ &
	892	DyOB((1.-DxOB)VB(IOB,KOBP,JOB+1)+DxOB* &
	893	VB(IOB+1,KOBP,JOB+1))))
	894
	895	! if(n.le.10) then
	896	! write(6,*)
	897	! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF2 at ',iob,job,kob, &
	898	! ' N = ',n,' inest = ',inest
	899	! write(6,*) 'VAROBS(2,N) = ',varobs(2,n)
	900	! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
	901	! write(6,*) 'VB(IOB,KOB,JOB) = ',VB(IOB,KOB,JOB)
	902	! write(6,*) 'ERRF(2,N) = ',errf(2,n)
	903	! write(6,*) 'vratiob = ',vratiob
	904	! write(6,*)
	905	! endif
	906
	907
	908	! Store model surface pressure (not the error!) at V point.
	909	ERRF(8,N)=.001( (1.-DYOB)PVG1 + DYOB*PVG2 )
	910
	911	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	912	ENDIF ! end IF( MP_LOCAL_VOBMASK(N) )
	913	#endif
	914	ENDDO ! END V-POINT LOOP OVER OBS
	915
	916	ENDIF ! end if(iswind.eq.1)
	917
	918	ENDIF ! ITYP=1: PROCESS V
	919
	920	!**********************************************************
	921	! PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
	922	!**********************************************************
	923	IF(ITYP.EQ.3) THEN
	924
	925	IF(ISTEMP.EQ.1 .OR. ISMOIS.EQ.1) THEN
	926	DO N=1,NSTA
	927	IOB=MAX0(1,IA(N))
	928	IOB=MIN0(IOB,ide-1)
	929	JOB=MAX0(1,IB(N))
	930	JOB=MIN0(JOB,jde-1)
	931	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	932	IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
	933	#endif
	934	KOB=MIN0(k_end,IC(N))
	935	KOBP=MIN0(KOB+1,K_END)
	936	DXOB=RA(N)
	937	DYOB=RB(N)
	938	DZOB=RC(N)
	939
	940	! This is to correct obs to model sigma level using reverse similarity theory
	941	if(rko(n).eq.1.0)then
	942	tratxob=((1.-DYOB)((1.-DXOB)tratx(IOB,JOB)+ &
	943	DXOB*tratx(IOB+1,JOB) &
	944	)+DYOB((1.-DXOB)tratx(IOB,JOB+1)+ &
	945	DXOB*tratx(IOB+1,JOB+1)))
	946	else
	947	tratxob=1.
	948	endif
	949
	950	!yliu
	951	if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
	952
	953	!ajb testing only
	954	if(iprt .and. n.eq.81) then
	955	write(6,*) 'POTENTIAL TEMP FOR N=81:'
	956	write(6,*)
	957	write(6,*) ' K THETA TEMPERATURE', &
	958	' PBASE'
	959	write(6,*)
	960	do k=k_end,1,-1
	961	press = pbase(iob,k,job)+pp(iob,k,job)
	962	ttemp = exp ( alog(300.+TB(IOB,k,JOB)) - &
	963	.2857143*alog(100000./press) )
	964	write(6,*) k,300.+TB(IOB,k,JOB),ttemp,pbase(iob,k,job)
	965	enddo
	966	endif
	967	!ajb end testing only
	968
	969	! TEMPERATURE ERROR
	970	! if(n.le.10) then
	971	! write(6,*) 'before: errf(3,n) = ',errf(3,n)
	972	! endif
	973	ERRF(3,N)=ERRF(3,N)+tratxobVAROBS(3,N)-((1.-DZOB) &
	974	((1.-DyOB)*((1.- &
	975	DxOB)(TB(IOB,KOB,JOB))+DxOB &
	976	(TB(IOB+1,KOB,JOB)))+DyOB((1.-DxOB) &
	977	(TB(IOB,KOB,JOB+1))+DxOB* &
	978	(TB(IOB+1,KOB,JOB+1))))+DZOB*((1.- &
	979	DyOB)((1.-DxOB)(TB(IOB,KOBP,JOB))+DxOB* &
	980	(TB(IOB+1,KOBP,JOB)))+DyOB((1.-DxOB) &
	981	(TB(IOB,KOBP,JOB+1))+DxOB* &
	982	(TB(IOB+1,KOBP,JOB+1)))))
	983
	984	! if(n.le.10) then
	985	! write(6,*)
	986	! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF3 at ',iob,job,kob, &
	987	! ' N = ',n,' inest = ',inest
	988	! write(6,*) 'VAROBS(3,N) = ',varobs(3,n)
	989	! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
	990	! write(6,*) 'TB(IOB,KOB,JOB) = ',TB(IOB,KOB,JOB)
	991	! write(6,*) 'TB(IOB+1,KOB,JOB) = ',TB(IOB+1,KOB,JOB)
	992	! write(6,*) 'TB(IOB,KOB,JOB+1) = ',TB(IOB,KOB,JOB+1)
	993	! write(6,*) 'TB(IOB+1,KOB,JOB+1) = ',TB(IOB+1,KOB,JOB+1)
	994	! write(6,*) 'TB(IOB,KOBP,JOB) = ',TB(IOB,KOBP,JOB)
	995	! write(6,*) 'TB(IOB+1,KOBP,JOB) = ',TB(IOB+1,KOBP,JOB)
	996	! write(6,*) 'TB(IOB,KOBP,JOB+1) = ',TB(IOB,KOBP,JOB+1)
	997	! write(6,*) 'TB(IOB+1,KOBP,JOB+1) = ',TB(IOB+1,KOBP,JOB+1)
	998	! write(6,*) 'tratxob = ',tratxob
	999	! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
	1000	! write(6,*) 'ERRF(3,N) = ',errf(3,n)
	1001	! write(6,*)
	1002	! endif
	1003
	1004
	1005	! MOISTURE ERROR
	1006	ERRF(4,N)=ERRF(4,N)+VAROBS(4,N)-((1.-DZOB)((1.-DyOB)((1.- &
	1007	DxOB)QVB(IOB,KOB,JOB)+DxOB &
	1008	QVB(IOB+1,KOB,JOB))+DyOB((1.-DxOB) &
	1009	QVB(IOB,KOB,JOB+1)+DxOB* &
	1010	QVB(IOB+1,KOB,JOB+1)))+DZOB*((1.- &
	1011	DyOB)((1.-DxOB)QVB(IOB,KOBP,JOB)+DxOB &
	1012	QVB(IOB+1,KOBP,JOB))+DyOB((1.-DxOB &
	1013	)QVB(IOB,KOBP,JOB+1)+DxOB &
	1014	QVB(IOB+1,KOBP,JOB+1))))
	1015
	1016	! Store model surface pressure (not the error!) at T-point
	1017	ERRF(6,N)= .001* &
	1018	((1.-DyOB)((1.-DxOB)pbase(IOB,1,JOB)+DxOB* &
	1019	pbase(IOB+1,1,JOB))+DyOB((1.-DxOB) &
	1020	pbase(IOB,1,JOB+1)+DxOB*pbase(IOB+1,1,JOB+1) ))
	1021
	1022	if(iprt .and. n.eq.81) then
	1023	write(6,*) 'ERRF(6,81) calculation:'
	1024	write(6,*) 'iob,job = ',iob,job
	1025	write(6,*) 'pbase(iob,1,job) = ',pbase(iob,1,job)
	1026	write(6,*) 'pbase(iob+1,1,job) = ',pbase(iob+1,1,job)
	1027	write(6,*) 'pbase(iob,1,job+1) = ',pbase(iob,1,job+1)
	1028	write(6,*) 'pbase(iob+1,1,job+1) = ',pbase(iob+1,1,job+1)
	1029	write(6,*) 'ERRF(6,81) = ',errf(6,n)
	1030	! call flush(6)
	1031	endif
	1032
	1033	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	1034	ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
	1035	#endif
	1036	ENDDO ! END T and QV LOOP OVER OBS
	1037
	1038	ENDIF !end if(istemp.eq.1 .or. ismois.eq.1)
	1039
	1040	!**********************************************************
	1041	! PROCESS SURFACE PRESSURE CROSS-POINT FIELD, IVAR=5,
	1042	! USING BILINEAR FOUR-POINT 2-D INTERPOLATION
	1043	!**********************************************************
	1044	IF(ISPSTR.EQ.1) THEN
	1045	DO N=1,NSTA
	1046	IOB=MAX0(1,IA(N))
	1047	IOB=MIN0(IOB,ide-1)
	1048	JOB=MAX0(1,IB(N))
	1049	JOB=MIN0(JOB,jde-1)
	1050	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	1051	IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
	1052	#endif
	1053	DXOB=RA(N)
	1054	DYOB=RB(N)
	1055	!ajb fix this (put in correct pressure calc for IOB,JOB here)
	1056	ERRF(5,N)=ERRF(5,N)+VAROBS(5,N)-((1.-DyOB)((1.-DxOB) &
	1057	pbase(IOB,1,JOB)+DxOBpbase(IOB+1,1,JOB))+DyOB &
	1058	((1.-DxOB)pbase(IOB,1,JOB+1)+DxOB &
	1059	pbase(IOB+1,1,JOB+1)))
	1060
	1061	if(n.eq.81) then
	1062	write(6,*) 'ERRF(5,81) calculation:'
	1063	write(6,*) 'iob,job = ',iob,job
	1064	write(6,*) 'pbase(iob,1,job) = ',pbase(iob,1,job)
	1065	write(6,*) 'pbase(iob+1,1,job) = ',pbase(iob+1,1,job)
	1066	write(6,*) 'pbase(iob,1,job+1) = ',pbase(iob,1,job+1)
	1067	write(6,*) 'pbase(iob+1,1,job+1) = ',pbase(iob+1,1,job+1)
	1068	write(6,*) 'errf(5,81) = ',errf(5,n)
	1069	! call flush(6)
	1070	endif
	1071
	1072	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	1073	ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
	1074	#endif
	1075
	1076	ENDDO
	1077
	1078	ENDIF ! end if(ispstr.eq.1)
	1079
	1080	ENDIF ! end if(ityp.eq.3)
	1081
	1082	ENDDO ! END BIG LOOP
	1083
	1084	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	1085	! Gather the errf values calculated by the processors with the obs.
	1086	CALL get_full_obs_vector(nsta, nerrf, niobf, mp_local_uobmask, &
	1087	mp_local_vobmask, mp_local_cobmask, errf)
	1088	#endif
	1089
	1090	! DIFFERENCE BETWEEN OBS AND FCST IS COMPLETED
	1091	IF(INEST.EQ.1)THEN
	1092	INPF=NPFI
	1093	ELSE
	1094	FNPF=IRATIO**LEVIDN(INEST)
	1095	INPF=FNPF*NPFI
	1096	ENDIF
	1097	! Gross error check for temperature. Set all vars bad.
	1098	do n=1,nsta
	1099	if((abs(errf(3,n)).gt.20.).and. &
	1100	(errf(3,n).gt.-800000.))then
	1101
	1102	errf(1,n)=-888888.
	1103	errf(2,n)=-888888.
	1104	errf(3,n)=-888888.
	1105	errf(4,n)=-888888.
	1106	varobs(1,n)=-888888.
	1107	varobs(2,n)=-888888.
	1108	varobs(3,n)=-888888.
	1109	varobs(4,n)=-888888.
	1110	endif
	1111	enddo
	1112
	1113	! For printout
	1114	! IF(MOD(KTAU,INPF).NE.0) THEN
	1115	! RETURN
	1116	! ENDIF
	1117
	1118	RETURN
	1119	END SUBROUTINE errob
	1120
	1121	SUBROUTINE upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, &
	1122	arrin, arrout)
	1123	!------------------------------------------------------------------------------
	1124	! PURPOSE: This subroutine interpolates a real 2D array defined over mass
	1125	! coordinate points, to U (momentum) points.
	1126	!
	1127	!------------------------------------------------------------------------------
	1128	IMPLICIT NONE
	1129
	1130	INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
	1131	INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
	1132	INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
	1133	INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
	1134	INTEGER, INTENT(IN) :: ids ! Starting i index for entire model domain
	1135	INTEGER, INTENT(IN) :: ide ! Ending i index for entire model domain
	1136	INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
	1137	INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
	1138	INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
	1139	INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
	1140	REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
	1141	REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on U points
	1142
	1143	! Local variables
	1144	integer :: i, j
	1145
	1146	! Do domain interior first
	1147	do j = j_start, j_end
	1148	do i = max(2,i_start), i_end
	1149	arrout(i,j) = 0.5*(arrin(i,j)+arrin(i-1,j))
	1150	enddo
	1151	enddo
	1152
	1153	! Do west-east boundaries
	1154	if(i_start .eq. ids) then
	1155	do j = j_start, j_end
	1156	arrout(i_start,j) = arrin(i_start,j)
	1157	enddo
	1158	endif
	1159	if(i_end .eq. ide-1) then
	1160	do j = j_start, j_end
	1161	arrout(i_end+1,j) = arrin(i_end,j)
	1162	enddo
	1163	endif
	1164
	1165	RETURN
	1166	END SUBROUTINE upoint
	1167
	1168	SUBROUTINE vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, &
	1169	arrin, arrout)
	1170	!------------------------------------------------------------------------------
	1171	! PURPOSE: This subroutine interpolates a real 2D array defined over mass
	1172	! coordinate points, to V (momentum) points.
	1173	!
	1174	!------------------------------------------------------------------------------
	1175	IMPLICIT NONE
	1176
	1177	INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
	1178	INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
	1179	INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
	1180	INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
	1181	INTEGER, INTENT(IN) :: jds ! Starting j index for entire model domain
	1182	INTEGER, INTENT(IN) :: jde ! Ending j index for entire model domain
	1183	INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
	1184	INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
	1185	INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
	1186	INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
	1187	REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
	1188	REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on V points
	1189
	1190	! Local variables
	1191	integer :: i, j
	1192
	1193	! Do domain interior first
	1194	do j = max(2,j_start), j_end
	1195	do i = i_start, i_end
	1196	arrout(i,j) = 0.5*(arrin(i,j)+arrin(i,j-1))
	1197	enddo
	1198	enddo
	1199
	1200	! Do south-north boundaries
	1201	if(j_start .eq. jds) then
	1202	do i = i_start, i_end
	1203	arrout(i,j_start) = arrin(i,j_start)
	1204	enddo
	1205	endif
	1206	if(j_end .eq. jde-1) then
	1207	do i = i_start, i_end
	1208	arrout(i,j_end+1) = arrin(i,j_end)
	1209	enddo
	1210	endif
	1211
	1212	RETURN
	1213	END SUBROUTINE vpoint
	1214
	1215	LOGICAL FUNCTION TILE_MASK(iloc, jloc, its, ite, jts, jte)
	1216	!------------------------------------------------------------------------------
	1217	! PURPOSE: Check to see if an i, j grid coordinate is in the tile index range.
	1218	!
	1219	! Returns: TRUE if the grid coordinate (ILOC,JLOC) is in the tile defined by
	1220	! tile-range indices (its,jts) and (ite,jte)
	1221	! FALSE otherwise.
	1222	!
	1223	!------------------------------------------------------------------------------
	1224	IMPLICIT NONE
	1225
	1226	INTEGER, INTENT(IN) :: iloc
	1227	INTEGER, INTENT(IN) :: jloc
	1228	INTEGER, INTENT(IN) :: its
	1229	INTEGER, INTENT(IN) :: ite
	1230	INTEGER, INTENT(IN) :: jts
	1231	INTEGER, INTENT(IN) :: jte
	1232
	1233	! Local variables
	1234	LOGICAL :: retval
	1235
	1236	TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND. &
	1237	jloc .LE. jte .AND. jloc .GE. jts )
	1238
	1239	RETURN
	1240	END FUNCTION TILE_MASK
	1241
	1242	!-----------------------------------------------------------------------
	1243	SUBROUTINE nudob(j, ivar, aten, inest, ifrest, ktau, ktaur, &
	1244	xtime, mu, msfx, msfy, nndgv, nerrf, niobf, maxdom, &
	1245	npfi, ionf, rinxy, twindo, levidn, &
	1246	parid, nstat, i_parent_start, j_parent_start, &
	1247	fdob, lev_in_ob, plfo, nlevs_ob, &
	1248	iratio, dx, dtmin, rio, rjo, rko, &
	1249	timeob, varobs, errf, pbase, ptop, pp, &
	1250	iswind, istemp, ismois, giv, git, giq, &
	1251	savwt, kpblt, nscan, &
	1252	iprt, &
	1253	ids,ide, jds,jde, kds,kde, & ! domain dims
	1254	ims,ime, jms,jme, kms,kme, & ! memory dims
	1255	its,ite, jts,jte, kts,kte ) ! tile dims
	1256
	1257	!-----------------------------------------------------------------------
	1258	USE module_model_constants
	1259	USE module_domain
	1260	!-----------------------------------------------------------------------
	1261	IMPLICIT NONE
	1262	!-----------------------------------------------------------------------
	1263	!
	1264	! PURPOSE: THIS SUBROUTINE GENERATES NUDGING TENDENCIES FOR THE J-TH
	1265	! VERTICAL SLICE (I-K PLANE) FOR FOUR-DIMENSIONAL DATA
	1266	! ASSIMILATION FROM INDIVIDUAL OBSERVATIONS. THE NUDGING
	1267	! TENDENCIES ARE FOUND FROM A ONE-PASS CALCULATION OF
	1268	! WEIGHTING FACTORS SIMILAR TO THE BENJAMIN-SEAMAN OBJECTIVE
	1269	! ANALYSIS. THIS SUBROUTINE IS DESIGNED FOR RAPID EXECUTION
	1270	! AND MINIMAL STORAGE REQUIREMENTS. ALGORITHMS SHOULD BE
	1271	! VECTORIZED WHEREVER POSSIBLE.
	1272	!
	1273	! HISTORY: Original author: MM5 version???
	1274	! 02/04/2004 - Creation of WRF version. Al Bourgeois
	1275	! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
	1276	!------------------------------------------------------------------------------
	1277	!
	1278	! NOTE: This routine was originally designed for MM5, which uses
	1279	! a nonstandard (I,J) coordinate system. For WRF, I is the
	1280	! east-west running coordinate, and J is the south-north
	1281	! running coordinate. So a "J-slab" here is west-east in
	1282	! extent, not south-north as for MM5. -ajb 06/10/2004
	1283	!
	1284	! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
	1285	! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
	1286	!
	1287	! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
	1288	! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
	1289	! TYPES OF FACTORS:
	1290	! 1) TIME WEIGHTING - ONLY OBSERVATIONS WITHIN A SELECTED
	1291	! TIME WINDOW (TWINDO) CENTERED AT THE CURRENT FORECAST
	1292	! TIME (XTIME) ARE USED. OBSERVATIONS CLOSEST TO
	1293	! XTIME ARE TIME-WEIGHTED MOST HEAVILY (TIMEWT)
	1294	! 2) VERTICAL WEIGHTING - NON-ZERO WEIGHTS (WTSIG) ARE
	1295	! CALCULATED WITHIN A VERTICAL REGION OF INFLUENCE
	1296	! (RINSIG).
	1297	! 3) HORIZONTAL WEIGHTING - NON-ZERO WEIGHTS (WTIJ) ARE
	1298	! CALCULATED WITHIN A RADIUS OF INFLUENCE (RINXY). THE
	1299	! VALUE OF RIN IS DEFINED IN KILOMETERS, AND CONVERTED
	1300	! TO GRID LENGTHS FOR THE APPROPRIATE MESH SIZE.
	1301	!
	1302	! THE FIVE FORECAST VARIABLES ARE PROCESSED BY CHANGING THE
	1303	! VALUE OF IVAR AS FOLLOWS:
	1304	! IVAR VARIABLE(TAU-1)
	1305	! ---- ---------------
	1306	! 1 U
	1307	! 2 V
	1308	! 3 T
	1309	! 4 QV
	1310	! 5 PSB(CROSS) REMOVED IN V3
	1311	! (6) PSB(DOT)
	1312	!
	1313	!-----------------------------------------------------------------------
	1314	!
	1315	! Description of input arguments.
	1316	!
	1317	!-----------------------------------------------------------------------
	1318
	1319	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
	1320	INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
	1321	INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
	1322	INTEGER, INTENT(IN) :: j ! south-north running coordinate.
	1323	INTEGER, INTENT(IN) :: ivar
	1324	INTEGER, INTENT(IN) :: inest ! domain index
	1325	LOGICAL, INTENT(IN) :: ifrest
	1326	INTEGER, INTENT(IN) :: ktau
	1327	INTEGER, INTENT(IN) :: ktaur
	1328	REAL, INTENT(IN) :: xtime ! forecast time in minutes
	1329	INTEGER, INTENT(IN) :: nndgv ! number of nudge variables
	1330	INTEGER, INTENT(IN) :: nerrf ! number of error fields
	1331	INTEGER, INTENT(IN) :: niobf ! number of observations
	1332	INTEGER, INTENT(IN) :: maxdom ! maximum number of domains
	1333	INTEGER, INTENT(IN) :: npfi
	1334	INTEGER, INTENT(IN) :: ionf
	1335	REAL, INTENT(IN) :: rinxy
	1336	REAL, INTENT(IN) :: twindo
	1337	INTEGER, INTENT(IN) :: levidn(maxdom) ! level of nest
	1338	INTEGER, INTENT(IN) :: parid(maxdom) ! parent domain id
	1339	INTEGER, INTENT(IN) :: nstat ! number of obs stations
	1340	INTEGER, INTENT(IN) :: i_parent_start(maxdom) ! Start i index in parent domain.
	1341	INTEGER, INTENT(IN) :: j_parent_start(maxdom) ! Start j index in parent domain.
	1342	TYPE(fdob_type), intent(inout) :: fdob
	1343	REAL, INTENT(IN) :: lev_in_ob(niobf) ! Level in sounding-type obs.
	1344	REAL, intent(IN) :: plfo(niobf)
	1345	REAL, INTENT(IN) :: nlevs_ob(niobf) ! Number of levels in sounding.
	1346	INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
	1347	REAL, INTENT(IN) :: dx ! This domain grid cell-size (m)
	1348	REAL, INTENT(IN) :: dtmin
	1349	REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
	1350	REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
	1351	REAL, INTENT(INOUT) :: rko(niobf) ! Obs vertical coordinate.
	1352	REAL, INTENT(IN) :: timeob(niobf)
	1353	REAL, INTENT(IN) :: varobs(nndgv,niobf)
	1354	REAL, INTENT(IN) :: errf(nerrf, niobf)
	1355	REAL, INTENT(IN) :: pbase( ims:ime, kms:kme ) ! Base pressure.
	1356	REAL, INTENT(IN) :: ptop
	1357	REAL, INTENT(IN) :: pp( ims:ime, kms:kme ) ! Pressure perturbation (Pa)
	1358	REAL, INTENT(IN) :: mu(ims:ime) ! Air mass on u, v, or mass-grid
	1359	REAL, INTENT(IN) :: msfx(ims:ime) ! Map scale (only used for vars u & v)
	1360	REAL, INTENT(IN) :: msfy(ims:ime) ! Map scale (only used for vars u & v)
	1361	INTEGER, intent(in) :: iswind ! Nudge flag for wind
	1362	INTEGER, intent(in) :: istemp ! Nudge flag for temperature
	1363	INTEGER, intent(in) :: ismois ! Nudge flag for moisture
	1364	REAL, intent(in) :: giv ! Coefficient for wind
	1365	REAL, intent(in) :: git ! Coefficient for temperature
	1366	REAL, intent(in) :: giq ! Coefficient for moisture
	1367	REAL, INTENT(INOUT) :: aten( ims:ime, kms:kme)
	1368	REAL, INTENT(INOUT) :: savwt( nndgv, ims:ime, kms:kme )
	1369	INTEGER, INTENT(IN) :: kpblt(its:ite)
	1370	INTEGER, INTENT(IN) :: nscan ! number of scans
	1371	LOGICAL, INTENT(IN) :: iprt ! print flag
	1372
	1373	! Local variables
	1374	integer :: mm(maxdom)
	1375	integer :: kobs ! k-lev below obs (for obs straddling pblt)
	1376	real :: ra(niobf)
	1377	real :: rb(niobf)
	1378	real :: psurf(niobf)
	1379	real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
	1380	real :: rscale(ims:ime) ! For converting to rho-coupled units.
	1381	! real :: tscale(ims:ime,kms:kme) ! For converting to potential temp. units.
	1382	real :: reserf(100)
	1383	character*40 name
	1384	character*3 chr_hr
	1385
	1386	!*** DECLARATIONS FOR IMPLICIT NONE
	1387	integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
	1388	integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
	1389	integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
	1390	integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
	1391	integer :: komin,komax,nn,nhi,nlo,nnjc
	1392	integer :: i_s,i_e
	1393	integer :: istq
	1394	real :: gfactor,rfactor,gridx,gridy,rindx,schnes,ris
	1395	real :: grfacx,grfacy
	1396	real :: fdtim,tw1,tw2,tconst,timewt,timewt2,ttim,dift,pob
	1397	real :: ri,rj,rx,ry,rsq,wtij,pdfac,erfivr,dk,slope,rinfac
	1398	real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
	1399
	1400	real :: scratch
	1401
	1402	! print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
	1403	! if(ivar.ne.4)return
	1404	!yliu start -- for multi-scans: NSCAN=0: original
	1405	! NSCAN=1: added a scan with a larger Ri and smaller G
	1406	! if(NSCAN.ne.0 .and. NSCAN.ne.1) stop
	1407	! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
	1408	if(NSCAN.ne.0) then
	1409	IF (iprt) write(6,*) 'SAVWT must be resized for NSCAN=1'
	1410	stop
	1411	endif
	1412	IPLO=0 + NSCAN*4
	1413	GFACTOR=1. + NSCAN*(-1. + 0.33333)
	1414	RFACTOR=1. + NSCAN*(-1. + 3.0)
	1415	!yliu end
	1416	! jc
	1417
	1418	! return if too close to j boundary
	1419	if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
	1420	! write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
	1421	! $ ' too close to boundary.'
	1422	return
	1423	endif
	1424	if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
	1425	! write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
	1426	! $ ' too close to boundary.'
	1427	return
	1428	endif
	1429
	1430	! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
	1431	ICUT=0
	1432	IF(INEST.GT.1)ICUT=1
	1433	i_s = max0(2+icut,its)
	1434	i_e = min0(ide-1-icut,ite)
	1435
	1436	IPL=IVAR + IPLO !yliu +IPLO
	1437
	1438	! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
	1439
	1440	INPF=(IRATIO*LEVIDN(INEST))NPFI
	1441	INFR=(IRATIO*LEVIDN(INEST))IONF
	1442
	1443	GRIDX=0.0
	1444	GRIDY=0.0
	1445	IGRID=0
	1446	IF(IVAR.GE.3)THEN
	1447	GRIDX=0.5
	1448	GRIDY=0.5
	1449	IGRID=1
	1450	ELSEIF(IVAR.eq.1) THEN
	1451	GRIDY=0.5
	1452	GRIDX=0.0
	1453	IGRID=1
	1454	ELSEIF(IVAR.eq.2) THEN
	1455	GRIDX=0.5
	1456	GRIDY=0.0
	1457	IGRID=1
	1458	ENDIF
	1459
	1460	! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
	1461	! KILOMETERS TO GRID LENGTHS, RINDX
	1462
	1463	RINDX=RINXY1000./DX RFACTOR !yliu *RFACTOR
	1464
	1465	! jc
	1466	! make horizontal radius vary per nest
	1467	! rindx=rindx/float(inest)
	1468	! yliu test1 -- En 3, 4
	1469	! rindx=rindx/float(3**(in-1)) !YLIU
	1470	! jc
	1471	! make horizontal radius vary per nest
	1472	! schnes=1/float(inest) !JC
	1473	! yliu test1 -- En 3, 4 !YLIU
	1474	schnes=1/float(3**(inest-1)) !JC
	1475	! reduce the Rinf in the nested grid proportionally
	1476	rindx=rindx*schnes
	1477	! rinfmn =1., rinfmx=2., pfree=50 in param.F
	1478	! yliu test: for upper-air data, use larger influence radii
	1479	! Essentially increase the slope -- the same radii
	1480	! at 500mb and above as the coarse mesh and the
	1481	! same small radii at sfc as that for sfc obs
	1482	fdob%rinfmx=2. *1.0 /schnes !YLIU
	1483	! rinfmx=1.2*1.0 /schnes !YLIU
	1484	! jc
	1485	RIS=RINDX*RINDX
	1486	IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
	1487	IF(MOD(KTAU,INFR).NE.0)GOTO 126
	1488	5 CONTINUE
	1489	IF (iprt) THEN
	1490	IF(J.EQ.10) write(6,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
	1491	ENDIF
	1492	6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,', &
	1493	'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2, &
	1494	' rindx=',f4.1)
	1495
	1496	!********************************************************************
	1497	!ajb_07012008 Setting ra and rb for the fine mesh her is now simple.
	1498	! Values are no longer calculated here based on the
	1499	! coarse mesh, since direct use of WRF map projections
	1500	! on each nest was implemented in subroutine in4dob.
	1501	!********************************************************************
	1502	! SET RA AND RB
	1503	DO N=1,NSTAT
	1504	RA(N)=RIO(N)-GRIDX
	1505	RB(N)=RJO(N)-GRIDY
	1506	ENDDO
	1507
	1508	! OUTPUT OBS PER GRID EVERY HOUR
	1509	if ( mod(xtime,60.).gt.56. .and. ivar.eq.3 .and. j.eq.10) then
	1510	IF (iprt) print *,'outputting obs number on grid ', &
	1511	inest,' at time=',xtime
	1512	write(chr_hr(1:3),'(i3)')nint(xtime/60.)
	1513	if(chr_hr(1:1).eq.' ')chr_hr(1:1)='0'
	1514	if(chr_hr(2:2).eq.' ')chr_hr(2:2)='0'
	1515	IF (iprt) print *,'chr_hr=',chr_hr(1:3),nint(xtime/60.)
	1516	open(91,file= &
	1517	'obs_g'//char(inest+ichar('0'))//'_'//chr_hr(1:3), &
	1518	form='FORMATted',status='unknown')
	1519	write(91,911)nstat
	1520	write(6,911)nstat
	1521	911 FORMAT('total obs=',i8)
	1522	nsthis=0
	1523	nsmetar=0
	1524	nsspeci=0
	1525	nsship=0
	1526	nssynop=0
	1527	nstemp=0
	1528	nspilot=0
	1529	nssatwnds=0
	1530	nssams=0
	1531	nsprofs=0
	1532	! print *,'ide,jde=',ide,jde
	1533	do jjjn=1,nstat
	1534	! DETERMINE THE TIME-WEIGHT FACTOR FOR N
	1535	FDTIM=XTIME-DTMIN
	1536	! CONVERT TWINDO AND TIMEOB FROM HOURS TO MINUTES:
	1537	TW1=TWINDO/2.*60.
	1538	TW2=TWINDO*60.
	1539	TCONST=1./TW1
	1540	TIMEWT2=0.0
	1541	TTIM=TIMEOB(jjjn)*60.
	1542	!***********TTIM=TARGET TIME IN MINUTES
	1543	DIFT=ABS(FDTIM-TTIM)
	1544	IF(DIFT.LE.TW1)TIMEWT2=1.0
	1545
	1546	IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
	1547	IF(FDTIM.LT.TTIM)TIMEWT2=(FDTIM-(TTIM-TW2))*TCONST
	1548	IF(FDTIM.GT.TTIM)TIMEWT2=((TTIM+TW2)-FDTIM)*TCONST
	1549	ENDIF
	1550
	1551	! print *,'timewt2=',timewt2,ttim,fdtim
	1552	if (ra(jjjn).ge.1. .and. rb(jjjn).ge.1. &
	1553	.and.ra(jjjn).le.real(ide) .and. rb(jjjn).le.real(jde) &
	1554	.and.timewt2.gt.0.) then
	1555	if(lev_in_ob(jjjn).eq.1.)nsthis=nsthis+1
	1556	if(plfo(jjjn).eq.1.)nsmetar=nsmetar+1
	1557	if(plfo(jjjn).eq.2.)nsspeci=nsspeci+1
	1558	if(plfo(jjjn).eq.3.)nsship=nsship+1
	1559	if(plfo(jjjn).eq.4.)nssynop=nssynop+1
	1560	if(plfo(jjjn).eq.5..and.lev_in_ob(jjjn).eq.1.) nstemp=nstemp+1
	1561	if(plfo(jjjn).eq.6..and.lev_in_ob(jjjn).eq.1.) nspilot=nspilot+1
	1562	if(plfo(jjjn).eq.7.)nssatwnds=nssatwnds+1
	1563	if(plfo(jjjn).eq.8.)nssams=nssams+1
	1564	if(plfo(jjjn).eq.9..and.lev_in_ob(jjjn).eq.1.) nsprofs=nsprofs+1
	1565	endif
	1566	enddo
	1567	write(91,912)nsthis
	1568	write(6,912)nsthis
	1569	912 FORMAT('total obs on this grid=',i8)
	1570	write(91,921)nsmetar
	1571	write(6,921)nsmetar
	1572	921 FORMAT('total metar obs on this grid=',i8)
	1573	write(91,922)nsspeci
	1574	write(6,922)nsspeci
	1575	922 FORMAT('total special obs on this grid=',i8)
	1576	write(91,923)nsship
	1577	write(6,923)nsship
	1578	923 FORMAT('total ship obs on this grid=',i8)
	1579	write(91,924)nssynop
	1580	write(6,924)nssynop
	1581	924 FORMAT('total synop obs on this grid=',i8)
	1582	write(91,925)nstemp
	1583	write(6,925)nstemp
	1584	925 FORMAT('total temp obs on this grid=',i8)
	1585	write(91,926)nspilot
	1586	write(6,926)nspilot
	1587	926 FORMAT('total pilot obs on this grid=',i8)
	1588	write(91,927)nssatwnds
	1589	write(6,927)nssatwnds
	1590	927 FORMAT('total sat-wind obs on this grid=',i8)
	1591	write(91,928)nssams
	1592	write(6,928)nssams
	1593	928 FORMAT('total sams obs on this grid=',i8)
	1594	write(91,929)nsprofs
	1595	write(6,929)nsprofs
	1596	929 FORMAT('total profiler obs on this grid=',i8)
	1597	close(91)
	1598	endif ! END OUTPUT OBS PER GRID EVERY HOUR
	1599
	1600
	1601	! INITIALIZE WEIGHTING ARRAYS TO ZERO
	1602	DO I=its,ite
	1603	DO K=1,kte
	1604	WT(I,K)=0.0
	1605	WT2ERR(I,K)=0.0
	1606	ENDDO
	1607	ENDDO
	1608
	1609	! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
	1610	! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
	1611	!
	1612	! COMPUTE P* AT OBS LOCATION (RA,RB). DO THIS AS SEPARATE VECTOR LOOP H
	1613	! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
	1614
	1615	! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
	1616	! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
	1617	! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
	1618	DO N=1,NSTAT
	1619	IF(IVAR.GE.3)THEN
	1620	PSURF(N)=ERRF(6,N)
	1621	ELSE
	1622	IF(IVAR.EQ.1)THEN
	1623	PSURF(N)=ERRF(7,N) ! U-points
	1624	ELSE
	1625	PSURF(N)=ERRF(8,N) ! V-points
	1626	ENDIF
	1627	ENDIF
	1628	ENDDO
	1629
	1630	! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
	1631	! J-STRIP
	1632
	1633	MAXJ=J+IFIX(RINDX*fdob%RINFMX+0.99) !ajb
	1634	MINJ=J-IFIX(RINDX*fdob%RINFMX+0.99) !ajb
	1635
	1636	! jc comment out this? want to use obs beyond the domain?
	1637	! MAXJ=MIN0(JL-IGRID,MAXJ) !yliu
	1638	! MINJ=MAX0(1,MINJ) !yliu
	1639
	1640	n=1
	1641
	1642	!***********************************************************************
	1643	DO nnn=1,NSTAT ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
	1644	!***********************************************************************
	1645	! Soundings are consecutive obs, but they need to be treated as a single
	1646	! entity. Thus change the looping to nnn, with n defined separately.
	1647
	1648
	1649	!yliu
	1650	! note for sfc data: nsndlev=1 and njcsnd=1
	1651	nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
	1652
	1653	! yliu start -- set together with the other parts
	1654	! test: do the sounding levels as individual obs
	1655	! nsndlev=1
	1656	! yliu end
	1657	njcsnd=nsndlev
	1658	! set pob here, to be used later
	1659	pob=varobs(5,n)
	1660
	1661	! print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
	1662	! print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
	1663	! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
	1664	! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP. THIS
	1665	! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
	1666	! ATION.
	1667
	1668	!yliu: Skip bad obs if it is sfc or single level sounding.
	1669	!yliu: Before this (020918), a snd will be skipped if its first
	1670	!yliu level has bad data- often true due to elevation
	1671
	1672	IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
	1673	! print *, " bad obs skipped"
	1674
	1675	ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
	1676	! print *, " skipped obs far away from this J-slice"
	1677
	1678	!----------------------------------------------------------------------
	1679	ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
	1680	!----------------------------------------------------------------------
	1681
	1682	! DETERMINE THE TIME-WEIGHT FACTOR FOR N
	1683	FDTIM=XTIME-DTMIN
	1684	! TWINDO IS IN MINUTES:
	1685	TW1=TWINDO/2.*60.
	1686	TW2=TWINDO*60.
	1687	TCONST=1./TW1
	1688	TIMEWT=0.0
	1689	TTIM=TIMEOB(N)*60.
	1690	!***********TTIM=TARGET TIME IN MINUTES
	1691	DIFT=ABS(FDTIM-TTIM)
	1692	IF(DIFT.LE.TW1)TIMEWT=1.0
	1693	IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
	1694	IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
	1695	IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
	1696	ENDIF
	1697
	1698	! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
	1699	! FOR THE VERTICAL WEIGHTING, WTSIG
	1700
	1701	! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
	1702	!ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
	1703	!ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
	1704
	1705	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
	1706	rko(n) = errf(9,n) !ajb 20021210
	1707	#endif
	1708	KOB=nint(RKO(N))
	1709	KOB=MIN0(kte,KOB)
	1710	KOB=MAX0(1,KOB)
	1711
	1712	! ASSIMILATE SURFACE LAYER DATA ON SIGMA
	1713	IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
	1714	DO K=1,kte
	1715	WTSIG(K)=0.0
	1716	ENDDO
	1717	! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
	1718	! WTSIG(1)=1.0
	1719	! WTSIG(2)=0.67
	1720	! WTSIG(3)=0.33
	1721	! KOMIN=3
	1722	! KOMAX=1
	1723	! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
	1724	! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
	1725	! fix this because kpblt at 1 and il is 0
	1726	MAXI=IFIX(RA(N)+0.99+RINDX)
	1727	MAXI=MIN0(ide-1,MAXI)
	1728	MINI=IFIX(RA(N)-RINDX-0.99)
	1729	MINI=MAX0(2,MINI)
	1730	!yliu start
	1731	! use also obs outside of this domain -- surface obs
	1732	! if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
	1733	! & RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
	1734	! print *, " skipped obs far away from this domain"
	1735	! currently can use obs within this domain or ones very close to (1/3
	1736	! influence of radius in the coarse domain) this
	1737	! domain. In later case, use BC column value to approximate the model value
	1738	! at obs point -- ERRF need model field in errob.F !!
	1739	if ( RA(N).GE.(0.-RINDX/3) &
	1740	.and. RA(N).LE.float(ide)+RINDX/3 &
	1741	.and. RB(N).GE.(0.-RINDX/3) &
	1742	.and. RB(N).LE.float(jde)+RINDX/3) then
	1743
	1744	! or use obs within this domain only
	1745	! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
	1746	! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
	1747	! print *, " skipped obs far outside of this domain"
	1748	! if(j.eq.3 .and. ivar.eq.3) then
	1749	! write(6,*) 'N = ',n,' exit 120 3'
	1750	! endif
	1751	!yliu end
	1752	!
	1753	! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
	1754	! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
	1755	! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
	1756	RJ=FLOAT(J)
	1757	RX=RJ-RB(N)
	1758	! WEIGHTS FOR THE 3-D VARIABLES
	1759	ERFIVR=ERRF(IVAR,N)
	1760	!
	1761	!JM I will be local, because it indexes into PDOC, WT, and others
	1762
	1763	! if((ivar.eq.1 .or. ivar.eq.3) .and. n.le.200) then
	1764	! write(6,'(a,i3,a,i3)')'SURF OBS NEAR: N = ',n,' nest = ',inest
	1765	! write(6,'(a,f10.3,a,f10.3,a,i3,a,i3,a,i3,a,i2)')
	1766	! $ ' RA =',RA(N),' RB =',RB(N),' J =',j,
	1767	! $ ' MINI =',MINI,' MAXI =',MAXI,' NEST =',inest
	1768	! endif
	1769
	1770	DO I=max0(its,MINI),min0(ite,MAXI)
	1771
	1772	RI=FLOAT(I)
	1773	RY=RI-RA(N)
	1774	RIS=RINDX*RINDX
	1775	RSQ=RXRX+RYRY
	1776	! DPRIM=SQRT(RSQ)
	1777	! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
	1778	! D=DPRIM+RINDXDCONABS(PSBO(N)-PDOC(I,J))
	1779	! DSQ=D*D
	1780	! WTIJ=(RIS-DSQ)/(RIS+DSQ)
	1781	wtij=(ris-rsq)/(ris+rsq)
	1782	scratch = (abs(psurf(n)-.001pbase(i,1))fdob%DCON)
	1783	pdfac=1.-AMIN1(1.0,scratch)
	1784	wtij=wtij*pdfac
	1785	WTIJ=AMAX1(0.0,WTIJ)
	1786
	1787	! try making sfc obs weighting go thru pbl
	1788	! jc kpbl is at dot or cross only - need to interpolate?
	1789	! wtsig(1)=1.
	1790	komax=max0(3,kpblt(i))
	1791
	1792	! jc arbitrary check here
	1793	IF (iprt) THEN
	1794	if (kpblt(i).gt.25 .and. ktau.ne.0) &
	1795	write(6,552)inest,i,j,kpblt(i)
	1796	552 FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
	1797	ENDIF
	1798
	1799	if(kpblt(i).gt.25) komax=3
	1800	komin=1
	1801	dk=float(komax)
	1802
	1803	do k=komin,komax
	1804
	1805	wtsig(k)=float(komax-k+1)/dk
	1806	WT(I,K)=WT(I,K)+TIMEWTWTSIG(K)WTIJ
	1807
	1808	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJWTIJWTSIG(K) &
	1809	WTSIG(K)ERFIVR
	1810
	1811	! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
	1812	!
	1813	! write(6,'(a,i2,a,f8.3,a,f8.3,a,f8.3,a,f8.3,a,f8.3)')
	1814	! 'Surface obs, after: k = ',k, &
	1815	! ' WT2ERR = ',wt2err(i,k), &
	1816	! ' TIMEWT = ',timewt, &
	1817	! ' WTIJ = ',wtij, &
	1818	! ' WSIG = ',wtsig(k), &
	1819	! ' ERFIVR = ',erfivr
	1820	! endif
	1821
	1822	enddo
	1823
	1824	ENDDO
	1825
	1826	! print *, " Surface "
	1827
	1828	endif ! end check for obs in domain
	1829	! END SURFACE-LAYER U OR V OBS NUDGING
	1830
	1831	ELSE
	1832	! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
	1833	!
	1834	! print *,'in upper air section'
	1835	! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
	1836	! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
	1837	! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
	1838	! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
	1839	!ajb SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
	1840
	1841	slope = (fdob%RINFMN-fdob%RINFMX)/(psurf(n)-fdob%PFREE)
	1842
	1843	RINFAC=SLOPEPOB+fdob%RINFMX-SLOPEfdob%pfree
	1844	RINFAC=AMAX1(RINFAC,fdob%RINFMN)
	1845	RINFAC=AMIN1(RINFAC,fdob%RINFMX)
	1846	!yliu: for multilevel upper-air data, take the maximum
	1847	! for the I loop.
	1848	if(nsndlev.gt.1) RINFAC = fdob%RINFMX
	1849	!yliu end
	1850
	1851	MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
	1852	MAXI=MIN0(ide-IGRID,MAXI)
	1853	MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
	1854	MINI=MAX0(1,MINI)
	1855
	1856	! yliu start
	1857	! use also obs outside of but close to this domain -- upr data
	1858	! if( RA(N).LT.(0.-RINFAC*RINDX)
	1859	! & .or. RA(N).GT.float(IL)+RINFAC*RINDX
	1860	! & .or. RB(N).LT.(0.-RINFAC*RINDX)
	1861	! & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then
	1862	! print *, " skipped obs far away from this I-range"
	1863	! currently can use obs within this domain or ones very close to (1/3
	1864	! influence of radius in the coarse domain) this
	1865	! domain. In later case, use BC column value to approximate the model value
	1866	! at obs point -- ERRF need model field in errob.F !!
	1867
	1868	!cc if (i.eq.39 .and. j.eq.34) then
	1869	!cc write(6,*) 'RA(N) = ',ra(n)
	1870	!cc write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
	1871	!cc endif
	1872	if( RA(N).GE.(0.-RINFAC*RINDX/3) &
	1873	.and. RA(N).LE.float(ide)+RINFAC*RINDX/3 &
	1874	.and. RB(N).GE.(0.-RINFAC*RINDX/3) &
	1875	.and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
	1876	! or use obs within this domain only
	1877	! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
	1878	! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
	1879	! print *, " skipped obs far outside of this domain"
	1880
	1881	! yliu end
	1882	! is this 2 needed here - kpbl not used?
	1883	! MINI=MAX0(2,MINI)
	1884
	1885	! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
	1886	! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
	1887	! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
	1888	RJ=FLOAT(J)
	1889	RX=RJ-RB(N)
	1890	! WEIGHTS FOR THE 3-D VARIABLES
	1891	!
	1892	ERFIVR=ERRF(IVAR,N)
	1893	! jc
	1894	nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
	1895	! yliu start
	1896	! test: do the sounding levels as individual obs
	1897	! nsndlev=1
	1898	! yliu end
	1899	njcsnd=nsndlev
	1900	!
	1901	DO I=max0(its,MINI),min0(ite,MAXI)
	1902	! jc
	1903	RI=FLOAT(I)
	1904	RY=RI-RA(N)
	1905	RIS=RINDXRINFACRINDX*RINFAC
	1906	RSQ=RXRX+RYRY
	1907	! yliu test: for upper-air data, keep D1 influence radii
	1908	! RIS=RIS /schnes /schnes
	1909	WTIJ=(RIS-RSQ)/(RIS+RSQ)
	1910	WTIJ=AMAX1(0.0,WTIJ)
	1911	! weight ob in vertical with +- 50 mb
	1912	! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
	1913	if(nsndlev.eq.1) then
	1914	rinprs=7.5
	1915	else
	1916	rinprs=3.0
	1917	endif
	1918	! yliu end
	1919	!
	1920	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	1921	! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
	1922	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	1923
	1924	if(nsndlev.eq.1)then
	1925	!----------------------------------------------------------------------
	1926	! --- HANDLE 1-LEVEL OBSERVATIONS ---
	1927	!----------------------------------------------------------------------
	1928
	1929	! if(I.eq.MINI) print *, " Single snd "
	1930	! ERFIVR is the residual (difference) between the ob and the model
	1931	! at that point. We can analyze that residual up and down.
	1932	! First find komin for ob.
	1933	!yliu start -- in the old code, komax and komin were reversed!
	1934	do k=kte,1,-1
	1935	pijk = .001*(pbase(i,k)+pp(i,k))
	1936	! print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
	1937	if(pijk.ge.(pob+rinprs)) then
	1938	komin=k
	1939	go to 325
	1940	endif
	1941	enddo
	1942	komin=1
	1943	325 continue
	1944	! now find komax for ob
	1945	do k=3,kte
	1946	pijk = .001*(pbase(i,k)+pp(i,k))
	1947	if(pijk.le.(pob-rinprs)) then
	1948	komax=k
	1949	go to 326
	1950	endif
	1951	enddo
	1952	komax=kte ! yliu 20050706
	1953	326 continue
	1954
	1955	! yliu: single-level upper-air data will act either above or below the PBL top
	1956	! ajb: Reset komin or komax. if kobs>kpblt, komin=kpblt+1, else komax=kpblt
	1957
	1958	if( (kpblt(i).le.komax) .and. (kpblt(i).ge.komin) ) then
	1959	kobs = 1
	1960	OBS_K: do k = komin, komax
	1961	if( pob .gt. .001*(pbase(i,k)+pp(i,k)) ) then
	1962	kobs = k
	1963	EXIT OBS_K
	1964	endif
	1965	enddo OBS_K
	1966
	1967	if(kobs.gt.kpblt(i)) then
	1968	komin=max0(kobs, komin) ! kobs here is kpblt(i)+1
	1969	else
	1970	komax=min0(kpblt(i), komax)
	1971	endif
	1972	endif
	1973	! yliu end
	1974	!
	1975	! print *,'1 level, komin,komax=',komin,komax
	1976	! if(i.eq.MINI) then
	1977	! print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
	1978	! ERFIVR=0
	1979	! endif
	1980	do k=1,kte
	1981	reserf(k)=0.0
	1982	wtsig(k)=0.0
	1983	enddo
	1984	!yliu end
	1985
	1986	!cc if (i.eq.39 .and. j.eq.34) then
	1987	!cc write(6,*) ' komin = ', komin,' komax = ',komax
	1988	!cc endif
	1989
	1990	do k=komin,komax
	1991	pijk = .001*(pbase(i,k)+pp(i,k))
	1992	reserf(k)=erfivr
	1993	wtsig(k)=1.-abs(pijk-pob)/rinprs
	1994	wtsig(k)=amax1(wtsig(k),0.0)
	1995	! print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
	1996	! Now calculate WT and WT2ERR for each i,j,k point cajb
	1997	WT(I,K)=WT(I,K)+TIMEWTWTIJwtsig(k)
	1998
	1999	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJWTIJ &
	2000	reserf(k)wtsig(k)wtsig(k)
	2001	enddo
	2002
	2003	else
	2004	!----------------------------------------------------------------------
	2005	! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
	2006	!----------------------------------------------------------------------
	2007	!yliu start
	2008	! if(I.eq.MINI) print *, " Multi-level snd "
	2009	! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
	2010	! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
	2011	if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
	2012	IF (iprt) THEN
	2013	print *, "n = ",n,"nsndlev = ",nsndlev
	2014	print *, "nlevs_ob,lev_in_ob", &
	2015	nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
	2016	print *, "in nudobs.F: sounding level messed up, stopping"
	2017	ENDIF
	2018	stop
	2019	endif
	2020	!yliu end
	2021	! This is for a multi-level observation
	2022	! The trick here is that the sounding is "one ob". You don't
	2023	! want multiple levels to each be treated like separate
	2024	! and independent observations.
	2025	! At each i,j want to interpolate sounding to the model levels at that
	2026	! particular point.
	2027	komin=1
	2028	komax=kte-2
	2029	! this loop goes to 1501
	2030	! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
	2031	!yliu start
	2032	do k=1,kte
	2033	reserf(k)=0.0
	2034	wtsig(k)=0.0
	2035	enddo
	2036	!yliu end
	2037
	2038	do k=komax,komin,-1
	2039
	2040	pijk = .001*(pbase(i,k)+pp(i,k))
	2041
	2042	! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
	2043	if(pijk.gt.varobs(5,n)) then
	2044	go to 1501
	2045	endif
	2046
	2047	! if sigma level pressure is .lt. than the highest ob level, don't interpolate
	2048	if(pijk.le.varobs(5,n+nsndlev-1)) then
	2049	go to 1501
	2050	endif
	2051
	2052	! now interpolate sounding to this k
	2053	! yliu start-- recalculate WTij for each k-level
	2054	!ajb SLOPE = (fdob%RINFMN-fdob%RINFMX)/(pdoc(i,j)+PTOP-fdob%PFREE)
	2055	slope = (fdob%RINFMN-fdob%RINFMX)/ (.001*pbase(i,1)-fdob%PFREE)
	2056	RINFAC=SLOPEpijk+fdob%RINFMX-SLOPEfdob%PFREE
	2057	RINFAC=AMAX1(RINFAC,fdob%RINFMN)
	2058	RINFAC=AMIN1(RINFAC,fdob%RINFMX)
	2059	RIS=RINDXRINFACRINDX*RINFAC
	2060	RSQ=RXRX+RYRY
	2061
	2062	! for upper-air data, keep D1 influence radii
	2063	! RIS=RIS /schnes /schnes
	2064	WTIJ=(RIS-RSQ)/(RIS+RSQ)
	2065	WTIJ=AMAX1(0.0,WTIJ)
	2066	! yliu end
	2067
	2068	! this loop goes to 1503
	2069	do nn=2,nsndlev
	2070	! only set pobhi if varobs(ivar) is ok
	2071	pobhi=-888888.
	2072
	2073	if(varobs(ivar,n+nn-1).gt.-800000. &
	2074	.and. varobs(5,n+nn-1).gt.-800000.) then
	2075	pobhi=varobs(5,n+nn-1)
	2076	nhi=n+nn-1
	2077	if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.20.) then
	2078	go to 1502 ! within 200mb of obs height
	2079	endif
	2080	endif
	2081
	2082	enddo
	2083
	2084	! did not find any ob above within 100 mb, so jump out
	2085	go to 1501
	2086	1502 continue
	2087
	2088	nlo=nhi-1
	2089	do nnjc=nhi-1,n,-1
	2090	if(varobs(ivar,nnjc).gt.-800000. &
	2091	.and. varobs(5,nnjc).gt.-800000.) then
	2092	poblo=varobs(5,nnjc)
	2093	nlo=nnjc
	2094	if(poblo.gt.pijk .and. abs(poblo-pijk).lt.20.) then
	2095	go to 1505 ! within 200mb of obs height
	2096	endif
	2097	endif
	2098	enddo
	2099	!yliu end --
	2100
	2101	! did not find any ob below within 200 mb, so jump out
	2102	go to 1501
	2103	1505 continue
	2104
	2105	! interpolate to model level
	2106	pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
	2107	reserf(k)=errf(ivar,nlo)+ &
	2108	(errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
	2109	wtsig(k)=1.
	2110
	2111	1501 continue
	2112
	2113	! now calculate WT and WT2ERR for each i,j,k point cajb
	2114	WT(I,K)=WT(I,K)+TIMEWTWTIJwtsig(k)
	2115
	2116	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJWTIJ &
	2117	reserf(k)wtsig(k)wtsig(k)
	2118
	2119	! if(ivar.eq.1 .and. i.eq.38 .and. j.eq.78) then
	2120	!
	2121	! if(wt(i,k) .ne. 0.0) then
	2122	! scratch = WT2ERR(I,K)/WT(I,K)
	2123	! else
	2124	! scratch = 999.
	2125	! endif
	2126	!
	2127	! write(6,'(a,i2,a,f8.3,a,f4.2,a,f7.4,a,f4.2,a,f5.3,a,f7.4)')
	2128	! $ 'Multi-level obs: k = ',k,
	2129	! $ ' WT2ERR = ',wt2err(i,k),
	2130	! $ ' WTIJ = ',wtij,
	2131	! $ ' RSF = ',reserf(k),
	2132	! $ ' WSIG = ',wtsig(k),
	2133	! $ ' WT = ',wt(i,k),
	2134	! $ ' W2EOWT = ',scratch
	2135	! endif
	2136
	2137
	2138	! end do k
	2139	enddo ! enddo k levels
	2140	! end multi-levels
	2141	endif ! end if(nsndlev.eq.1)
	2142	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	2143	! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
	2144	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
	2145	!
	2146	ENDDO ! END DO MINI,MAXI LOOP
	2147
	2148	endif ! check for obs in domain
	2149
	2150	! END OF NUDGING TO OBS ON PRESSURE LEVELS
	2151
	2152	ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
	2153
	2154	!----------------------------------------------------------------------
	2155	ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
	2156	!----------------------------------------------------------------------
	2157
	2158	! n=n+1
	2159	n=n+njcsnd
	2160
	2161	!yliu 1202 continue
	2162	if(n.gt.nstat)then
	2163	! print *,'n,nstat=',n,nstat,ivar,j
	2164	go to 1203
	2165	endif
	2166	! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
	2167
	2168	!***********************************************************************
	2169	ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
	2170	!***********************************************************************
	2171
	2172	1203 continue
	2173
	2174	! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
	2175	! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
	2176	! THE ATEN ARRAY
	2177	! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
	2178	! THEY ARE USED BELOW IN THE DENOMINATOR.
	2179	DO K=kts,kte
	2180	DO I=its,ite
	2181	IF(WT(I,K).EQ.0)THEN
	2182	WT2ERR(I,K)=0.0
	2183	ENDIF
	2184	IF(WT(I,K).EQ.0)THEN
	2185	WT(I,K)=1.0
	2186	ENDIF
	2187	ENDDO
	2188	ENDDO
	2189
	2190	126 CONTINUE
	2191
	2192	IF(IVAR.GE.3)GOTO 170
	2193	! this is for u,v
	2194	! 3-D DOT POINT TENDENCIES
	2195
	2196	! Calculate scales for converting nudge factor from u (v)
	2197	! to rho_u (or rho_v) units.
	2198
	2199	IF (IVAR == 1) THEN
	2200	call calc_rcouple_scales(mu,msfy,rscale,ims,ime,its,ite)
	2201	ELSE IF (IVAR == 2) THEN
	2202	call calc_rcouple_scales(mu,msfx,rscale,ims,ime,its,ite)
	2203	END IF
	2204
	2205	DO K=1,kte
	2206
	2207	DO I=i_s,i_e
	2208
	2209	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
	2210	W2EOWT=WT2ERR(I,K)/WT(I,K)
	2211	ELSE
	2212	W2EOWT=SAVWT(IPL,I,K)
	2213	ENDIF
	2214
	2215	! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
	2216	! scratch = GIVRSCALE(I)W2EOWTfdob%TFACIISWIND*GFACTOR
	2217	! write(6,*) 'ATEN calc: k = ',k
	2218	! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
	2219	! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
	2220	! $ ' W2EOWT = ',w2eowt
	2221	! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
	2222	! $ ' GFACTOR = ',gfactor
	2223	! endif
	2224	!
	2225	! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
	2226	! scratch = GIVRSCALE(I)W2EOWTfdob%TFACIISWIND*GFACTOR
	2227	! write(6,*) 'ATEN calc: k = ',k
	2228	! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
	2229	! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
	2230	! $ ' W2EOWT = ',w2eowt
	2231	! write(6,*) 'TFACI = ',fdob%tfaci,' ISWIND = ',iswind,
	2232	! $ ' GFACTOR = ',gfactor
	2233	! endif
	2234
	2235	ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
	2236	W2EOWTfdob%TFACI &
	2237	ISWIND GFACTOR !yliu *GFACTOR
	2238
	2239	! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
	2240	! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
	2241	! endif
	2242	! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
	2243	! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
	2244	! endif
	2245
	2246	ENDDO
	2247	ENDDO
	2248
	2249	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
	2250	DO K=1,kte
	2251	DO I=its,ite
	2252	SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
	2253	ENDDO
	2254	ENDDO
	2255	ENDIF
	2256
	2257	RETURN
	2258
	2259	170 CONTINUE
	2260
	2261	! 3-D CROSS-POINT TENDENCIES
	2262	! this is for t (ivar=3) and q (ivsr=4)
	2263	IF(3-IVAR.LT.0)THEN
	2264	GITQ=GIQ
	2265	ELSE
	2266	GITQ=GIT
	2267	ENDIF
	2268	IF(3-IVAR.LT.0)THEN
	2269	ISTQ=ISMOIS
	2270	ELSE
	2271	ISTQ=ISTEMP
	2272	ENDIF
	2273
	2274	DO K=1,kte
	2275	DO I=i_s,i_e
	2276	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
	2277	W2EOWT=WT2ERR(I,K)/WT(I,K)
	2278	ELSE
	2279	W2EOWT=SAVWT(IPL,I,K)
	2280	ENDIF
	2281
	2282	! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
	2283	! scratch = GITQMU(I)W2EOWTfdob%TFACIISTQ*GFACTOR
	2284	! write(6,*) 'ATEN calc: k = ',k
	2285	! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
	2286	! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
	2287	! $ ' W2EOWT = ',w2eowt
	2288	! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
	2289	! $ ' GFACTOR = ',gfactor
	2290	! endif
	2291	!
	2292	! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
	2293	! scratch = GITQMU(I)W2EOWTfdob%TFACIISTQ*GFACTOR
	2294	! write(6,*) 'ATEN calc: k = ',k
	2295	! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
	2296	! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
	2297	! $ ' W2EOWT = ',w2eowt
	2298	! write(6,*) ' TFACI = ',fdob%tfaci,' ISTQ = ',istq,
	2299	! $ ' GFACTOR = ',gfactor
	2300	! endif
	2301
	2302	ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
	2303	W2EOWTfdob%TFACIISTQ GFACTOR !yliu *GFACTOR
	2304
	2305	! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.29) then
	2306	! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
	2307	! endif
	2308	! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
	2309	! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
	2310	! endif
	2311
	2312	ENDDO
	2313	ENDDO
	2314
	2315	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
	2316	DO K=1,kte
	2317	DO I=its,ite
	2318	SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
	2319	ENDDO
	2320	ENDDO
	2321	ENDIF
	2322
	2323	RETURN
	2324	END SUBROUTINE nudob
	2325
	2326	SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
	2327	!-----------------------------------------------------------------------
	2328	IMPLICIT NONE
	2329	!-----------------------------------------------------------------------
	2330
	2331	INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
	2332	INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
	2333	REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
	2334	REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
	2335	REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
	2336
	2337	! Local variables
	2338	integer :: i
	2339
	2340	! Calculate scales to be used for producing rho-coupled nudging factors.
	2341	do i = its,ite
	2342	rscale(i) = a(i)/msf(i)
	2343	enddo
	2344
	2345	RETURN
	2346	END SUBROUTINE calc_rcouple_scales
	2347
	2348	!ajb: Not used
	2349	SUBROUTINE set_real_array(rscale, value, ims,ime, its,ite)
	2350	!-----------------------------------------------------------------------
	2351	IMPLICIT NONE
	2352	!-----------------------------------------------------------------------
	2353
	2354	INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
	2355	INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
	2356	REAL, INTENT(IN) :: value ! Constant array value
	2357	REAL, INTENT(OUT) :: rscale( ims:ime ) ! Output array
	2358
	2359	! Local variables
	2360	integer :: i
	2361
	2362	! Set array to constant value
	2363	do i = its,ite
	2364	rscale(i) = value
	2365	enddo
	2366
	2367	RETURN
	2368	END SUBROUTINE set_real_array
	2369
	2370	!ajb: Not used
	2371	SUBROUTINE calc_pottemp_scales(ivar, rcp, pb, p, tscale, &
	2372	ims,ime, its,ite, &
	2373	kms,kme, kts,kte)
	2374	!-----------------------------------------------------------------------
	2375	IMPLICIT NONE
	2376	!-----------------------------------------------------------------------
	2377
	2378	INTEGER, INTENT(IN) :: ims,ime, kms,kme ! Memory dimensions
	2379	INTEGER, INTENT(IN) :: its,ite, kts,kte ! Tile dimensions
	2380	INTEGER, INTENT(IN) :: ivar ! Variable identifier
	2381	REAL, INTENT(IN) :: rcp ! Constant (2./7.)
	2382	REAL, INTENT(IN) :: pb(ims:ime, kms:kme) ! Base pressure (Pa) array
	2383	REAL, INTENT(IN) :: p(ims:ime, kms:kme) ! Pressure pert. (Pa) array
	2384	REAL, INTENT(OUT) :: tscale(ims:ime, kms:kme) ! Scales for pot. temp.
	2385	! Local variables
	2386	integer :: i,k
	2387
	2388	if(ivar.eq.3) then
	2389
	2390	! Calculate scales to be used for producing potential temperature nudging factors.
	2391	do k = kts,kte
	2392	do i = its,ite
	2393	tscale(i,k) = ( 1000000. / ( pb(i,k)+p(i,k)) )**rcp
	2394	enddo
	2395	enddo
	2396	else
	2397	! Set to 1. for moisture scaling.
	2398	do k = kts,kte
	2399	do i = its,ite
	2400	tscale(i,k) = 1.0
	2401	enddo
	2402	enddo
	2403	endif
	2404
	2405	RETURN
	2406	END SUBROUTINE calc_pottemp_scales
	2407	#endif
	2408
	2409	END MODULE module_fddaobs_rtfdda
	2410

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format