Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_fddaobs_rtfdda.F @ 3988

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

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

Download in other formats:

Original Format