Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
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

Line
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