Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_fddaobs_rtfdda.F @ 1544

Last change on this file since 1544 was 1, checked in by lfita, 10 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 117.1 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(nudge_opt, maxdom, inest, parid, &
42	idynin, dtramp, fdaend, restart, &
43	twindo_cg, twindo, itimestep, &
44	no_pbl_nudge_uv, &
45	no_pbl_nudge_t, &
46	no_pbl_nudge_q, &
47	sfc_scheme_horiz, sfc_scheme_vert, &
48	maxsnd_gap, &
49	sfcfact, sfcfacr, dpsmx, &
50	nudge_wind, nudge_temp, nudge_mois, &
51	nudgezfullr1_uv, nudgezrampr1_uv, &
52	nudgezfullr2_uv, nudgezrampr2_uv, &
53	nudgezfullr4_uv, nudgezrampr4_uv, &
54	nudgezfullr1_t, nudgezrampr1_t, &
55	nudgezfullr2_t, nudgezrampr2_t, &
56	nudgezfullr4_t, nudgezrampr4_t, &
57	nudgezfullr1_q, nudgezrampr1_q, &
58	nudgezfullr2_q, nudgezrampr2_q, &
59	nudgezfullr4_q, nudgezrampr4_q, &
60	nudgezfullmin, nudgezrampmin, nudgezmax, &
61	xlat, xlong, &
62	start_year, start_month, start_day, &
63	start_hour, start_minute, start_second, &
64	p00, t00, tlp, &
65	znu, p_top, &
66	#if ( EM_CORE == 1 )
67	fdob, &
68	#endif
69	iprt, &
70	ids,ide, jds,jde, kds,kde, &
71	ims,ime, jms,jme, kms,kme, &
72	its,ite, jts,jte, kts,kte)
73	!-----------------------------------------------------------------------
74	! This routine does initialization for real time fdda obs-nudging.
75	!
76	!-----------------------------------------------------------------------
77	USE module_model_constants, ONLY : g, r_d
78	USE module_domain
79	USE module_dm, ONLY : wrf_dm_min_real
80	!-----------------------------------------------------------------------
81	IMPLICIT NONE
82	!-----------------------------------------------------------------------
83
84	!=======================================================================
85	! Definitions
86	!-----------
87	INTEGER, intent(in) :: maxdom
88	INTEGER, intent(in) :: nudge_opt(maxdom)
89	INTEGER, intent(in) :: ids,ide, jds,jde, kds,kde, &
90	ims,ime, jms,jme, kms,kme, &
91	its,ite, jts,jte, kts,kte
92	INTEGER, intent(in) :: inest
93	INTEGER, intent(in) :: parid(maxdom)
94	INTEGER, intent(in) :: idynin ! flag for dynamic initialization
95	REAL, intent(in) :: dtramp ! time period for idynin ramping (min)
96	REAL, intent(in) :: fdaend(maxdom) ! nudging end time for domain (min)
97	LOGICAL, intent(in) :: restart
98	REAL, intent(in) :: twindo_cg ! time window on coarse grid
99	REAL, intent(in) :: twindo
100	INTEGER, intent(in) :: itimestep
101	INTEGER , INTENT(IN) :: no_pbl_nudge_uv(maxdom) ! flags for no wind nudging in pbl
102	INTEGER , INTENT(IN) :: no_pbl_nudge_t(maxdom) ! flags for no temperature nudging in pbl
103	INTEGER , INTENT(IN) :: no_pbl_nudge_q(maxdom) ! flags for no moisture nudging in pbl
104	INTEGER , INTENT(IN) :: sfc_scheme_horiz ! horizontal spreading scheme for surf obs (wrf or orig mm5)
105	INTEGER , INTENT(IN) :: sfc_scheme_vert ! vertical spreading scheme for surf obs (orig or regime vif)
106	REAL , INTENT(IN) :: maxsnd_gap ! max allowed pressure gap in soundings, for interp (centibars)
107	REAL, intent(in) :: sfcfact ! scale factor applied to time window for surface obs
108	REAL, intent(in) :: sfcfacr ! scale fac applied to horiz rad of infl for sfc obs
109	REAL, intent(in) :: dpsmx ! max press change allowed within hor rad of infl
110	INTEGER , INTENT(IN) :: nudge_wind(maxdom) ! wind-nudging flag
111	INTEGER , INTENT(IN) :: nudge_temp(maxdom) ! temperature-nudging flag
112	INTEGER , INTENT(IN) :: nudge_mois(maxdom) ! moisture-nudging flag
113	REAL, INTENT(IN) :: nudgezfullr1_uv ! vert infl fcn, regime=1 full-wt hght, winds
114	REAL, INTENT(IN) :: nudgezrampr1_uv ! vert infl fcn, regime=1 ramp down hght, winds
115	REAL, INTENT(IN) :: nudgezfullr2_uv ! vert infl fcn, regime=2 full-wt hght, winds
116	REAL, INTENT(IN) :: nudgezrampr2_uv ! vert infl fcn, regime=2 ramp down hght, winds
117	REAL, INTENT(IN) :: nudgezfullr4_uv ! vert infl fcn, regime=4 full-wt hght, winds
118	REAL, INTENT(IN) :: nudgezrampr4_uv ! vert infl fcn, regime=4 ramp down hght, winds
119	REAL, INTENT(IN) :: nudgezfullr1_t ! vert infl fcn, regime=1 full-wt hght, temp
120	REAL, INTENT(IN) :: nudgezrampr1_t ! vert infl fcn, regime=1 ramp down hght, temp
121	REAL, INTENT(IN) :: nudgezfullr2_t ! vert infl fcn, regime=2 full-wt hght, temp
122	REAL, INTENT(IN) :: nudgezrampr2_t ! vert infl fcn, regime=2 ramp down hght, temp
123	REAL, INTENT(IN) :: nudgezfullr4_t ! vert infl fcn, regime=4 full-wt hght, temp
124	REAL, INTENT(IN) :: nudgezrampr4_t ! vert infl fcn, regime=4 ramp down hght, temp
125	REAL, INTENT(IN) :: nudgezfullr1_q ! vert infl fcn, regime=1 full-wt hght, mois
126	REAL, INTENT(IN) :: nudgezrampr1_q ! vert infl fcn, regime=1 ramp down hght, mois
127	REAL, INTENT(IN) :: nudgezfullr2_q ! vert infl fcn, regime=2 full-wt hght, mois
128	REAL, INTENT(IN) :: nudgezrampr2_q ! vert infl fcn, regime=2 ramp down hght, mois
129	REAL, INTENT(IN) :: nudgezfullr4_q ! vert infl fcn, regime=4 full-wt hght, mois
130	REAL, INTENT(IN) :: nudgezrampr4_q ! vert infl fcn, regime=4 ramp down hght, mois
131	REAL, INTENT(IN) :: nudgezfullmin ! min dpth thru which vert infl fcn remains 1.0 (m)
132	REAL, INTENT(IN) :: nudgezrampmin ! min dpth thru which vif decreases 1.0 to 0.0 (m)
133	REAL, INTENT(IN) :: nudgezmax ! max dpth in which vif is nonzero (m)
134	REAL, INTENT(IN) :: xlat ( ims:ime, jms:jme ) ! latitudes on mass-point grid
135	REAL, INTENT(IN) :: xlong( ims:ime, jms:jme ) ! longitudes on mass-point grid
136	INTEGER, intent(in) :: start_year ! Model start year
137	INTEGER, intent(in) :: start_month ! Model start month
138	INTEGER, intent(in) :: start_day ! Model start day
139	INTEGER, intent(in) :: start_hour ! Model start hour
140	INTEGER, intent(in) :: start_minute ! Model start minute
141	INTEGER, intent(in) :: start_second ! Model start second
142	REAL, INTENT(IN) :: p00 ! base state pressure
143	REAL, INTENT(IN) :: t00 ! base state temperature
144	REAL, INTENT(IN) :: tlp ! base state lapse rate
145	REAL, INTENT(IN) :: p_top ! pressure at top of model
146	REAL, INTENT(IN) :: znu( kms:kme ) ! eta values on half (mass) levels
147	#if ( EM_CORE == 1 )
148	TYPE(fdob_type), intent(inout) :: fdob
149	#endif
150	LOGICAL, intent(in) :: iprt ! Flag enabling printing warning messages
151
152	! Local variables
153	logical :: nudge_flag ! nudging flag for this nest
154	integer :: ktau ! current timestep
155	integer :: nest ! loop counter
156	integer :: idom ! domain id
157	integer :: parent ! parent domain
158	real :: conv ! 180/pi
159	real :: tl1 ! Lambert standard parallel 1
160	real :: tl2 ! Lambert standard parallel 2
161	real :: xn1
162	real :: known_lat ! Latitude of domain point (i,j)=(1,1)
163	real :: known_lon ! Longitude of domain point (i,j)=(1,1)
164	character(len=200) :: msg ! Argument to wrf_message
165	real :: z_at_p( kms:kme ) ! height at p levels
166	integer :: i,j,k ! loop counters
167
168	#if ( EM_CORE == 1 )
169
170	! Check to see if the nudging flag has been set. If not,
171	! simply RETURN.
172	nudge_flag = (nudge_opt(inest) .eq. 1)
173	if (.not. nudge_flag) return
174
175	call wrf_message("")
176	write(msg,fmt='(a,i2)') ' OBSERVATION NUDGING IS ACTIVATED FOR MESH ',inest
177	call wrf_message(msg)
178
179	ktau = itimestep
180	if(restart) then
181	fdob%ktaur = ktau
182	else
183	fdob%ktaur = 0
184	endif
185
186	! Create character string containing model starting-date
187	CALL date_string(start_year, start_month, start_day, start_hour, &
188	start_minute, start_second, fdob%sdate)
189
190	! Set flag for nudging on pressure (not sigma) surfaces.
191	fdob%iwtsig = 0
192
193	!**************************************************************************
194	! * Initialize datend for dynamic initialization (ajb added 08132008) *
195	!**************************************************************************
196	! Set ending nudging date (used with dynamic ramp-down) to zero.
197	fdob%datend = 0.
198	fdob%ieodi = 0
199
200	! Check for dynamic initialization flag
201	if(idynin.eq.1)then
202	! Set datend to time in minutes after which data are assumed to have ended.
203	if(dtramp.gt.0.)then
204	fdob%datend = fdaend(inest) - dtramp
205	else
206	fdob%datend = fdaend(inest)
207	endif
208	if(iprt) then
209	call wrf_message("")
210	write(msg,fmt='(a,i3,a)') &
211	' * DYNAMIC-INITIALIZATION OPTION FOR INEST = ', inest, ' *'
212	call wrf_message(msg)
213	write(msg,*) ' FDAEND,DATEND,DTRAMP: ',fdaend(inest),fdob%datend,dtramp
214	call wrf_message(msg)
215	call wrf_message("")
216	endif
217	endif
218	! * end dynamic initialization section *
219	!**************************************************************************
220
221	! Store flags for surface obs spreading schemes
222	if(sfc_scheme_horiz.eq.1) then
223	call wrf_message('MM5 scheme selected for horizontal spreading of surface obs')
224	elseif (sfc_scheme_horiz.eq.0) then
225	call wrf_message('WRF scheme selected for horizontal spreading of surface obs')
226	else
227	write(msg,fmt='(a,i3)') 'Unknown h-spreading scheme for surface obs: ',sfc_scheme_horiz
228	call wrf_message(msg)
229	call wrf_message("Valid selections: 0=WRF scheme, 1=Original MM5 scheme")
230	call wrf_error_fatal ( 'fddaobs_init: module_fddaobs_rtfdda STOP' )
231	endif
232
233	if(sfc_scheme_vert.eq.1) then
234	call wrf_message('Original simple scheme selected for vertical spreading of surface obs')
235	elseif (sfc_scheme_vert.eq.0) then
236	call wrf_message("Regime-based VIF scheme selected for vertical spreading of surface obs")
237	else
238	write(msg,fmt='(a,i3)') 'Unknown v-spreading scheme for surface obs: ',sfc_scheme_vert
239	call wrf_message(msg)
240	call wrf_message("Valid selections: 0=Regime-based VIF scheme, 1=Original simple scheme")
241	call wrf_error_fatal ( 'fddaobs_init: module_fddaobs_rtfdda STOP' )
242	endif
243	fdob%sfc_scheme_horiz = sfc_scheme_horiz
244	fdob%sfc_scheme_vert = sfc_scheme_vert
245
246
247	! Store temporal and spatial scales
248	fdob%sfcfact = sfcfact
249	fdob%sfcfacr = sfcfacr
250
251	! Set time window.
252	fdob%window = twindo
253	call wrf_message("")
254	write(msg,fmt='(a,i3)') '*** TIME WINDOW SETTINGS FOR NEST ',inest
255	call wrf_message(msg)
256	write(msg,fmt='(a,f6.3,2(a,f5.3))') ' TWINDO (hrs) = ',twindo, &
257	' SFCFACT = ',sfcfact,' SFCFACR = ',sfcfacr
258	call wrf_message(msg)
259	call wrf_message("")
260
261	if(inest.eq.1) then
262	if(twindo .eq. 0.) then
263	if(iprt) then
264	call wrf_message("")
265	write(msg,) '** WARNING: TWINDO=0 on the coarse domain.'
266	call wrf_message(msg)
267	write(msg,) '** Did you forget to set twindo in the fdda namelist?'
268	call wrf_message(msg)
269	call wrf_message("")
270	endif
271	endif
272	else ! nest
273	if(twindo .eq. 0.) then
274	fdob%window = twindo_cg
275	if(iprt) then
276	call wrf_message("")
277	write(msg,fmt='(a,i2)') 'WARNING: TWINDO=0. for nest ',inest
278	call wrf_message(msg)
279	write(msg,fmt='(a,f12.5,a)') 'Default to coarse-grid value of ', twindo_cg,' hrs'
280	call wrf_message(msg)
281	call wrf_message("")
282	endif
283	endif
284	endif
285
286	! Initialize flags.
287
288	fdob%domain_tot=0
289	do nest=1,maxdom
290	fdob%domain_tot = fdob%domain_tot + nudge_opt(nest)
291	end do
292
293	! Calculate and store dcon from dpsmx
294	if(dpsmx.gt.0.) then
295	fdob%dpsmx = dpsmx
296	fdob%dcon = 1.0/fdob%dpsmx
297	else
298	call wrf_error_fatal('fddaobs_init: Namelist variable dpsmx must be greater than zero!')
299	endif
300
301	! Calculate and store base-state heights at half (mass) levels.
302	CALL get_base_state_height_column( p_top, p00, t00, tlp, g, r_d, znu, &
303	fdob%base_state, kts, kte, kds,kde, kms,kme )
304
305	! Initialize flags for nudging within PBL.
306	fdob%nudge_uv_pbl = .true.
307	fdob%nudge_t_pbl = .true.
308	fdob%nudge_q_pbl = .true.
309	if(no_pbl_nudge_uv(inest) .eq. 1) fdob%nudge_uv_pbl = .false.
310	if(no_pbl_nudge_t(inest) .eq. 1) fdob%nudge_t_pbl = .false.
311	if(no_pbl_nudge_q(inest) .eq. 1) fdob%nudge_q_pbl = .false.
312
313	if(no_pbl_nudge_uv(inest) .eq. 1) then
314	fdob%nudge_uv_pbl = .false.
315	write(msg,*) ' --> Obs nudging for U/V is turned off in PBL'
316	call wrf_message(msg)
317	endif
318	if(no_pbl_nudge_t(inest) .eq. 1) then
319	fdob%nudge_t_pbl = .false.
320	write(msg,*) ' --> Obs nudging for T is turned off in PBL'
321	call wrf_message(msg)
322	endif
323	if(no_pbl_nudge_q(inest) .eq. 1) then
324	fdob%nudge_q_pbl = .false.
325	write(msg,*) ' --> Obs nudging for Q is turned off in PBL'
326	call wrf_message(msg)
327	endif
328
329	! Store max allowed pressure gap for interpolating between soundings
330	fdob%max_sndng_gap = maxsnd_gap
331	write(msg,fmt='(a,f6.1)') &
332	'*** MAX PRESSURE GAP (cb) for interpolation between soundings = ',maxsnd_gap
333	call wrf_message(msg)
334	call wrf_message("")
335
336	! Initialize vertical influence fcn for LML obs
337	if(sfc_scheme_vert.eq.0) then
338	fdob%vif_uv(1) = nudgezfullr1_uv
339	fdob%vif_uv(2) = nudgezrampr1_uv
340	fdob%vif_uv(3) = nudgezfullr2_uv
341	fdob%vif_uv(4) = nudgezrampr2_uv
342	fdob%vif_uv(5) = nudgezfullr4_uv
343	fdob%vif_uv(6) = nudgezrampr4_uv
344	fdob%vif_t (1) = nudgezfullr1_t
345	fdob%vif_t (2) = nudgezrampr1_t
346	fdob%vif_t (3) = nudgezfullr2_t
347	fdob%vif_t (4) = nudgezrampr2_t
348	fdob%vif_t (5) = nudgezfullr4_t
349	fdob%vif_t (6) = nudgezrampr4_t
350	fdob%vif_q (1) = nudgezfullr1_q
351	fdob%vif_q (2) = nudgezrampr1_q
352	fdob%vif_q (3) = nudgezfullr2_q
353	fdob%vif_q (4) = nudgezrampr2_q
354	fdob%vif_q (5) = nudgezfullr4_q
355	fdob%vif_q (6) = nudgezrampr4_q
356
357	! Sanity checks
358	if(nudgezmax.le.0.) then
359	write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezmax MUST BE GREATER THAN ZERO.'
360	call wrf_message(msg)
361	write(msg,*) 'THE NAMELIST VALUE IS',nudgezmax
362	call wrf_message(msg)
363	call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgemax value' )
364	endif
365	if(nudgezfullmin.lt.0.) then
366	write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezfullmin MUST BE NONNEGATIVE.'
367	call wrf_message(msg)
368	write(msg,*) 'THE NAMELIST VALUE IS',nudgezfullmin
369	call wrf_message(msg)
370	call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgefullmin value' )
371	endif
372	if(nudgezrampmin.lt.0.) then
373	write(msg,*) 'STOP! OBS NAMELIST INPUT obs_nudgezrampmin MUST BE NONNEGATIVE.'
374	call wrf_message(msg)
375	write(msg,*) 'THE NAMELIST VALUE IS',nudgezrampmin
376	call wrf_message(msg)
377	call wrf_error_fatal ( 'fddaobs_init: STOP on bad obs_nudgerampmin value' )
378	endif
379	if(nudgezmax.lt.nudgezfullmin+nudgezrampmin) then
380	write(msg,*) 'STOP! INCONSISTENT OBS NAMELIST INPUTS.'
381	call wrf_message(msg)
382	write(msg,fmt='(3(a,f12.3))') 'obs_nudgezmax = ',nudgezmax, &
383	' obs_nudgezfullmin = ',nudgezfullmin, &
384	' obs_nudgezrampmin = ',nudgezrampmin
385	call wrf_message(msg)
386	write(msg,*) 'REQUIRE NUDGEZMAX >= NUDGEZFULLMIN + NUDGEZRAMPMIN'
387	call wrf_message(msg)
388	call wrf_error_fatal ( 'fddaobs_init: STOP on inconsistent namelist values' )
389	endif
390
391	fdob%vif_fullmin = nudgezfullmin
392	fdob%vif_rampmin = nudgezrampmin
393	fdob%vif_max = nudgezmax
394
395	! Check to make sure that if nudgzfullmin > 0, then it must be at least as large as the
396	! first model half-level will be anywhere in the domain at any time within the simulation.
397	! We use 1.1 times the base-state value fdob%base_state(1) for this purpose.
398
399	if(nudgezfullmin.gt.0.0) then
400	if(nudgezfullmin .lt. 1.1*fdob%base_state(1)) then
401	fdob%vif_fullmin = 1.1*fdob%base_state(1)
402	endif
403	endif
404
405	! Print vertical weight info only if wind, temperature, or moisture nudging is requested.
406	if( (nudge_wind(inest).eq.1) .or. (nudge_temp(inest).eq.1) &
407	.or. (nudge_mois(inest).eq.1) ) then
408	call wrf_message("")
409	write(msg,fmt='(a,i2,a)') ' *** SETUP DESCRIPTION FOR SURFACE OBS NUDGING ON MESH ',inest,' :'
410	call wrf_message(msg)
411
412	call wrf_message("")
413	write(msg,fmt='(a,i5,a)') ' NUDGEZMAX: The maximum height at which nudging will be'// &
414	' applied from surface obs is ', nint(nudgezmax),' m AGL.'
415	call wrf_message(msg)
416	call wrf_message("")
417	write(msg,fmt='(a,i3,a)') ' NUDGEZFULLMIN: The minimum height of full nudging weight'// &
418	' for surface obs is ', nint(fdob%vif_fullmin),' m.'
419	call wrf_message(msg)
420	if(nudgezfullmin.lt.fdob%vif_fullmin) then
421	write(msg,fmt='(a,i3,a)') ' *WARNING*: NUDGEZFULLMIN has been increased from'// &
422	' the user-input value of ',nint(nudgezfullmin),' m.'
423	call wrf_message(msg)
424	write(msg,fmt='(a,i3,a)') ' to ensure that at least the bottom model level is'// &
425	' included in full nudging.'
426	call wrf_message(msg)
427	endif
428	call wrf_message("")
429	write(msg,fmt='(a,i3,a)') ' NUDGEZRAMPMIN: The minimum height to ramp from full to no'// &
430	' nudging for surface obs is ', nint(nudgezrampmin),' m.'
431	call wrf_message(msg)
432	call wrf_message("")
433	endif !endif either wind, temperature, or moisture nudging is requested
434
435	! Print vif settings
436	if(nudge_wind(inest) .eq. 1) then
437	call print_vif_var('wind', fdob%vif_uv, nudgezfullmin, nudgezrampmin)
438	call wrf_message("")
439	endif
440	if(nudge_temp(inest) .eq. 1) then
441	call print_vif_var('temp', fdob%vif_t, nudgezfullmin, nudgezrampmin)
442	call wrf_message("")
443	endif
444	if(nudge_mois(inest) .eq. 1) then
445	call print_vif_var('mois', fdob%vif_q, nudgezfullmin, nudgezrampmin)
446	call wrf_message("")
447	endif
448
449	if( (nudge_wind(inest).eq.1) .or. (nudge_temp(inest).eq.1) &
450	.or. (nudge_mois(inest).eq.1) ) then
451	write(msg,fmt='(a,i2)') ' *** END SETUP DESCRIPTION FOR SURFACE OBS NUDGING ON MESH ',inest
452	call wrf_message(msg)
453	call wrf_message("")
454	endif
455	endif !endif(sfc_scheme_vert.EQ.0)
456
457	! Set parameters.
458	fdob%pfree = 50.0
459	fdob%rinfmn = 1.0
460	fdob%rinfmx = 2.0
461
462	! Get known lat and lon and store these on all processors in fdob structure, for
463	! later use in projection x-formation to map (lat,lon) to (x,y) for each obs.
464	IF (its .eq. 1 .AND. jts .eq. 1) then
465	known_lat = xlat(1,1)
466	known_lon = xlong(1,1)
467	ELSE
468	known_lat = 9999.
469	known_lon = 9999.
470	END IF
471	fdob%known_lat = wrf_dm_min_real(known_lat)
472	fdob%known_lon = wrf_dm_min_real(known_lon)
473
474	! Calculate the nest levels, levidn. Note that each nest
475	! must know the nest levels levidn(maxdom) of each domain.
476	do nest=1,maxdom
477
478	! Initialize nest level for each domain.
479	if (nest .eq. 1) then
480	fdob%levidn(nest) = 0 ! Mother domain has nest level 0
481	else
482	fdob%levidn(nest) = 1 ! All other domains start with 1
483	endif
484	idom = nest
485	100 parent = parid(idom) ! Go up the parent tree
486
487	if (parent .gt. 1) then ! If not up to mother domain
488	fdob%levidn(nest) = fdob%levidn(nest) + 1
489	idom = parent
490	goto 100
491	endif
492	enddo
493
494
495	! fdob%LML_OBS_HT1_LEV = kte
496	! HT1: do k = kte, kts, -1
497	! if( LML_HT1 .gt. z_at_p(k) ) then
498	! fdob%LML_OBS_HT1_LEV = k
499	! EXIT HT1
500	! endif
501	! enddo HT1
502
503	! fdob%LML_OBS_HT2_LEV = kte
504	! HT2: do k = kte, kts, -1
505	! if( LML_HT2 .gt. z_at_p(k) ) then
506	! fdob%LML_OBS_HT2_LEV = k
507	! EXIT HT2
508	! endif
509	! enddo HT2
510	RETURN
511	#endif
512	END SUBROUTINE fddaobs_init
513
514	#if ( EM_CORE == 1 )
515	!-----------------------------------------------------------------------
516	SUBROUTINE errob(inest, ub, vb, tb, t0, qvb, pbase, pp, rovcp, &
517	z, &
518	uratx, vratx, tratx, kpbl, &
519	nndgv, nerrf, niobf, maxdom, &
520	levidn, parid, nstat, nstaw, &
521	iswind, istemp, ismois, ispstr, &
522	timeob, rio, rjo, rko, &
523	varobs, errf, ktau, xtime, &
524	iratio, npfi, &
525	prt_max, prt_freq, iprt, &
526	obs_prt, stnid_prt, lat_prt, lon_prt, &
527	mlat_prt, mlon_prt, &
528	ids,ide, jds,jde, kds,kde, &
529	ims,ime, jms,jme, kms,kme, &
530	its,ite, jts,jte, kts,kte )
531
532	!-----------------------------------------------------------------------
533	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
534	USE module_dm, ONLY : get_full_obs_vector, wrf_dm_sum_real
535	#else
536	USE module_dm, ONLY : wrf_dm_sum_real
537	#endif
538	USE module_model_constants, ONLY : rcp
539
540	!-----------------------------------------------------------------------
541	IMPLICIT NONE
542	!-----------------------------------------------------------------------
543	!
544	! PURPOSE: THIS SUBROUTINE CALCULATES THE DIFFERENCE BETWEEN THE
545	! OBSERVED VALUES AND THE FORECAST VALUES AT THE OBSERVATION
546	! POINTS. THE OBSERVED VALUES CLOSEST TO THE CURRENT
547	! FORECAST TIME (XTIME) WERE DETERMINED IN SUBROUTINE
548	! IN4DOB AND STORED IN ARRAY VAROBS. THE DIFFERENCES
549	! CALCULATED BY SUBROUTINE ERROB WILL BE STORED IN ARRAY
550	! ERRF FOR THE NSTA OBSERVATION LOCATIONS. MISSING
551	! OBSERVATIONS ARE DENOTED BY THE DUMMY VALUE 99999.
552	!
553	! HISTORY: Original author: MM5 version???
554	! 02/04/2004 - Creation of WRF version. Al Bourgeois
555	! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
556	!------------------------------------------------------------------------------
557
558	! THE STORAGE ORDER IN ERRF IS AS FOLLOWS:
559	! IVAR VARIABLE TYPE(TAU-1)
560	! ---- --------------------
561	! 1 U error at obs loc
562	! 2 V error at obs loc
563	! 3 T error at obs loc
564	! 4 Q error at obs loc
565	! 5 Vertical layer index for PBL top at IOB, JOB
566	! 6 Model surface press at obs loc (T-points)
567	! 7 Model surface press at obs loc (U-points)
568	! 8 Model surface press at obs loc (V-points)
569	! 9 RKO at U-points
570
571	!-----------------------------------------------------------------------
572	!
573	! Description of input arguments.
574	!
575	!-----------------------------------------------------------------------
576
577	INTEGER, INTENT(IN) :: inest ! Domain index.
578	INTEGER, INTENT(IN) :: nndgv ! Number of nudge variables.
579	INTEGER, INTENT(IN) :: nerrf ! Number of error fields.
580	INTEGER, INTENT(IN) :: niobf ! Number of observations.
581	INTEGER, INTENT(IN) :: maxdom ! Maximum number of domains.
582	INTEGER, INTENT(IN) :: levidn(maxdom) ! Level of nest.
583	INTEGER, INTENT(IN) :: parid(maxdom) ! Id of parent grid.
584	INTEGER, INTENT(IN) :: ktau ! Model time step index
585	REAL, INTENT(IN) :: xtime ! Forecast time in minutes
586	INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
587	INTEGER, INTENT(IN) :: npfi ! Coarse-grid diagnostics freq.
588	INTEGER, INTENT(IN) :: prt_max ! Max number of obs to print.
589	INTEGER, INTENT(IN) :: prt_freq ! Frequency of obs to print.
590	LOGICAL, INTENT(IN) :: iprt ! Print flag
591	INTEGER, INTENT(IN) :: obs_prt(prt_max) ! Print obs indices
592	INTEGER, INTENT(IN) :: stnid_prt(40,prt_max) ! Print obs station ids
593	REAL, INTENT(IN) :: lat_prt(prt_max) ! Print obs latitude
594	REAL, INTENT(IN) :: lon_prt(prt_max) ! Print obs longitude
595	REAL, INTENT(IN) :: mlat_prt(prt_max) ! Print model lat at obs loc
596	REAL, INTENT(IN) :: mlon_prt(prt_max) ! Print model lon at obs loc
597	INTEGER, INTENT(IN) :: nstat ! # stations held for use
598	INTEGER, INTENT(IN) :: nstaw ! # stations in current window
599	INTEGER, intent(in) :: iswind
600	INTEGER, intent(in) :: istemp
601	INTEGER, intent(in) :: ismois
602	INTEGER, intent(in) :: ispstr
603	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
604	INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
605	INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
606
607	REAL, INTENT(IN) :: ub( ims:ime, kms:kme, jms:jme )
608	REAL, INTENT(IN) :: vb( ims:ime, kms:kme, jms:jme )
609	REAL, INTENT(IN) :: tb( ims:ime, kms:kme, jms:jme )
610	REAL, INTENT(IN) :: t0
611	REAL, INTENT(IN) :: qvb( ims:ime, kms:kme, jms:jme )
612	REAL, INTENT(IN) :: pbase( ims:ime, kms:kme, jms:jme )
613	REAL, INTENT(IN) :: pp( ims:ime, kms:kme, jms:jme ) ! Press. perturbation (Pa)
614	REAL, INTENT(IN) :: rovcp
615	REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! Ht above sl on half-levels
616	REAL, INTENT(IN) :: uratx( ims:ime, jms:jme ) ! U to U10 ratio on mass points.
617	REAL, INTENT(IN) :: vratx( ims:ime, jms:jme ) ! V to V10 ratio on mass points.
618	REAL, INTENT(IN) :: tratx( ims:ime, jms:jme ) ! T to TH2 ratio on mass points.
619	INTEGER,INTENT(IN) :: kpbl( ims:ime, jms:jme ) ! Vertical layer index for PBL top
620	REAL, INTENT(IN) :: timeob(niobf) ! Obs time (hrs)
621	REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
622	REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
623	REAL, INTENT(INOUT) :: rko(niobf) ! Obs bottom-top coordinate
624	REAL, INTENT(INOUT) :: varobs(nndgv, niobf)
625	REAL, INTENT(INOUT) :: errf(nerrf, niobf)
626
627	! Local variables
628	INTEGER :: iobmg(niobf) ! Obs i-coord on mass grid
629	INTEGER :: jobmg(niobf) ! Obs j-coord on mass grid
630	INTEGER :: ia(niobf)
631	INTEGER :: ib(niobf)
632	INTEGER :: ic(niobf)
633	REAL :: pbbo(kds:kde) ! column base pressure (cb) at obs loc.
634	REAL :: ppbo(kds:kde) ! column pressure perturbation (cb) at obs loc.
635
636	REAL :: ra(niobf)
637	REAL :: rb(niobf)
638	REAL :: rc(niobf)
639	REAL :: dxobmg(niobf) ! Interp. fraction (x dir) referenced to mass-grid
640	REAL :: dyobmg(niobf) ! Interp. fraction (y dir) referenced to mass-grid
641	INTEGER MM(MAXDOM)
642	INTEGER NNL
643	real :: uratio( ims:ime, jms:jme ) ! U to U10 ratio on momentum points.
644	real :: vratio( ims:ime, jms:jme ) ! V to V10 ratio on momentum points.
645	real :: pug1,pug2,pvg1,pvg2
646	character(len=200) :: msg ! Argument to wrf_message
647
648	! Define staggers for U, V, and T grids, referenced from non-staggered grid.
649	real, parameter :: gridx_t = 0.5 ! Mass-point x stagger
650	real, parameter :: gridy_t = 0.5 ! Mass-point y stagger
651	real, parameter :: gridx_u = 0.0 ! U-point x stagger
652	real, parameter :: gridy_u = 0.5 ! U-point y stagger
653	real, parameter :: gridx_v = 0.5 ! V-point x stagger
654	real, parameter :: gridy_v = 0.0 ! V-point y stagger
655
656	real :: dummy = 99999.
657
658	real :: pbhi, pphi
659	real :: obs_pottemp ! Potential temperature at observation
660
661	!*** DECLARATIONS FOR IMPLICIT NONE
662	integer nsta,ivar,n,ityp
663	integer iob,job,kob,iob_ms,job_ms
664	integer k,kbot,nml,nlb,nle
665	integer iobm,jobm,iobp,jobp,kobp,inpf,i,j
666	integer i_start,i_end,j_start,j_end ! loop ranges for uratio,vratio calc.
667	integer k_start,k_end
668	integer ips ! For printing obs information
669	integer pnx ! obs index for printout
670
671	real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
672	real pob
673	real hob
674	real uratiob,vratiob,tratiob,tratxob,fnpf
675
676	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
677	LOGICAL MP_LOCAL_DUMMASK(NIOBF) ! Mask for work to be done on this processor
678	LOGICAL MP_LOCAL_UOBMASK(NIOBF) ! Dot-point mask
679	LOGICAL MP_LOCAL_VOBMASK(NIOBF) ! Dot-point mask
680	LOGICAL MP_LOCAL_COBMASK(NIOBF) ! Cross-point mask
681	#endif
682
683	! LOGICAL, EXTERNAL :: TILE_MASK
684
685	NSTA=NSTAT
686
687	! FIRST, DETERMINE THE GRID TYPE CORRECTION FOR U-momentum, V-momentum,
688	! AND MASS POINTS, AND WHEN INEST=2, CONVERT THE STORED COARSE MESH INDICES
689	! TO THE FINE MESH INDEX EQUIVALENTS
690
691	! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
692
693	if (iprt) then
694	write(msg,fmt='(a,i5,a,i2,a,i5,a)') '++++++CALL ERROB AT KTAU = ', &
695	KTAU,' AND INEST = ',INEST,': NSTA = ',NSTAW,' ++++++'
696	call wrf_message(msg)
697	endif
698
699	ERRF = 0.0 ! Zero out errf array
700
701	! Set up loop bounds for this grid's boundary conditions
702	i_start = max( its-1,ids )
703	i_end = min( ite+1,ide-1 )
704	j_start = max( jts-1,jds )
705	j_end = min( jte+1,jde-1 )
706	k_start = kts
707	k_end = min( kte, kde-1 )
708
709	DO ITYP=1,3 ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP
710
711	! Set grid stagger
712	IF(ITYP.EQ.1) THEN ! U-POINT CASE
713	GRIDX = gridx_u
714	GRIDY = gridy_u
715	ELSE IF(ITYP.EQ.2) THEN ! V-POINT CASE
716	GRIDX = gridx_v
717	GRIDY = gridy_v
718	ELSE ! MASS-POINT CASE
719	GRIDX = gridx_t
720	GRIDY = gridy_t
721	ENDIF
722
723	! Compute URATIO and VRATIO fields on momentum (u,v) points.
724	IF(ityp.eq.1)THEN
725	call upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, uratx, uratio)
726	ELSE IF (ityp.eq.2) THEN
727	call vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, vratx, vratio)
728	ENDIF
729
730	IF(INEST.EQ.1) THEN ! COARSE MESH CASE...
731	DO N=1,NSTA
732	RA(N)=RIO(N)-GRIDX
733	RB(N)=RJO(N)-GRIDY
734	IA(N)=RA(N)
735	IB(N)=RB(N)
736	IOB=MAX0(1,IA(N))
737	IOB=MIN0(IOB,ide-1)
738	JOB=MAX0(1,IB(N))
739	JOB=MIN0(JOB,jde-1)
740	DXOB=RA(N)-FLOAT(IA(N))
741	DYOB=RB(N)-FLOAT(IB(N))
742
743	! Save mass-point arrays for computing rko for all var types
744	if(ityp.eq.1) then
745	iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
746	jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
747	dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
748	dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
749	endif
750	iob_ms = iobmg(n)
751	job_ms = jobmg(n)
752	dxob_ms = dxobmg(n)
753	dyob_ms = dyobmg(n)
754
755
756	! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION
757	! This is tricky: Initialize pob to zero in all procs. Only one proc actually
758	! calculates pob. If this is an obs to be converted from height-to-pressure, then
759	! by definition, varobs(5,n) will initially have the missing value -888888. After
760	! the pob calculation, pob will be zero in all procs except the one that calculated
761	! it, and so pob is dm_summed over all procs and stored into varobs(5,n). So on
762	! subsequent passes, the dm_sum will not occur because varobs(5,n) will not have a
763	! missing value. If this is not an obs-in-height,
764
765	pob = 0.0
766
767	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
768	! Set mask for obs to be handled by this processor
769	MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
770
771	IF ( MP_LOCAL_DUMMASK(N) ) THEN
772	#endif
773
774	! Interpolate pressure to obs location column and convert from Pa to cb.
775
776	do k = kds, kde
777	pbbo(k) = .001*( &
778	(1.-DYOB_MS)( (1.-DXOB_MS)pbase(IOB_MS,K,JOB_MS) + &
779	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
780	DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
781	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
782	ppbo(k) = .001*( &
783	(1.-DYOB_MS)( (1.-DXOB_MS)pp(IOB_MS,K,JOB_MS) + &
784	DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
785	DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
786	DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
787	enddo
788
789	!ajb 20040119: Note based on bugfix for dot/cross points split across processors,
790	!ajb which was actually a serial code fix: The ityp=2 (v-points) and
791	!ajb ityp=3 (mass-points) cases should not use the ityp=1 (u-points)
792	!ajb case rko! This is necessary for bit-for-bit reproducability
793	!ajb with the parallel run. (coarse mesh)
794
795	if(abs(rko(n)+99).lt.1.)then
796	pob = varobs(5,n)
797
798	if(pob .eq.-888888.) then
799	hob = varobs(6,n)
800	if(hob .gt. -800000. ) then
801	pob = ht_to_p( hob, ppbo, pbbo, z, iob_ms, job_ms, &
802	dxob_ms, dyob_ms, k_start, k_end, kds,kde, &
803	ims,ime, jms,jme, kms,kme )
804	endif
805	endif
806
807	if(pob .gt.-800000.)then
808	do k=k_end-1,1,-1
809	kbot = k
810	if(pob .le. pbbo(k)+ppbo(k)) then
811	goto 199
812	endif
813	enddo
814	199 continue
815
816	pphi = ppbo(kbot+1)
817	pbhi = pbbo(kbot+1)
818
819	rko(n) = real(kbot+1)- &
820	( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
821
822	rko(n)=max(rko(n),1.0)
823	endif
824	endif
825
826	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
827	ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
828	#endif
829
830	! ajb 20090423: If obs-in-height, varobs(5,n) is sum of pob (see note above).
831	if(varobs(5,n) .eq. -888888. .and. varobs(6,n) .gt. -800000.) then
832	varobs(5,n) = wrf_dm_sum_real ( pob )
833	endif
834
835	RC(N)=RKO(N)
836
837	ENDDO ! END COARSE MESH LOOP OVER NSTA
838
839	ELSE ! FINE MESH CASE
840
841	!**********************************************************************
842	!ajb_07012008: Conversion of obs coordinates to the fine mesh here
843	!ajb is no longer necessary, since implementation of the WRF map pro-
844	!ajb jections (to each nest directly) is implemented in sub in4dob.
845	!**********************************************************************
846	!ajb
847	!ajb GET (I,J,K) OF OBSERVATIONS ON FINE MESH VALUES.
848	DO N=1,NSTA
849	RA(N)=RIO(N)-GRIDX ! ajb added 07012008
850	RB(N)=RJO(N)-GRIDY ! ajb added 07012008
851	IA(N)=RA(N)
852	IB(N)=RB(N)
853	IOB=MAX0(1,IA(N))
854	IOB=MIN0(IOB,ide-1)
855	JOB=MAX0(1,IB(N))
856	JOB=MIN0(JOB,jde-1)
857	DXOB=RA(N)-FLOAT(IA(N))
858	DYOB=RB(N)-FLOAT(IB(N))
859
860	! Save mass-point arrays for computing rko for all var types
861	if(ityp.eq.1) then
862	iobmg(n) = MIN0(MAX0(1,int(RIO(n)-gridx_t)),ide-1)
863	jobmg(n) = MIN0(MAX0(1,int(RJO(n)-gridy_t)),jde-1)
864	dxobmg(n) = RIO(N)-gridx_t-FLOAT(int(RIO(N)-gridx_t))
865	dyobmg(n) = RJO(N)-gridy_t-FLOAT(int(RJO(N)-gridy_t))
866	endif
867	iob_ms = iobmg(n)
868	job_ms = jobmg(n)
869	dxob_ms = dxobmg(n)
870	dyob_ms = dyobmg(n)
871
872	! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION (see COARSE MESH comments)
873	pob = 0.0
874
875	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
876	! Set mask for obs to be handled by this processor
877	MP_LOCAL_DUMMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
878
879	IF ( MP_LOCAL_DUMMASK(N) ) THEN
880	#endif
881
882	! Interpolate pressure to obs location column and convert from Pa to cb.
883
884	do k = kds, kde
885	pbbo(k) = .001*( &
886	(1.-DYOB_MS)( (1.-DXOB_MS)pbase(IOB_MS,K,JOB_MS) + &
887	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS) ) + &
888	DYOB_MS* ( (1.-DXOB_MS)*pbase(IOB_MS,K,JOB_MS+1) + &
889	DXOB_MS *pbase(IOB_MS+1,K,JOB_MS+1) ) )
890	ppbo(k) = .001*( &
891	(1.-DYOB_MS)( (1.-DXOB_MS)pp(IOB_MS,K,JOB_MS) + &
892	DXOB_MS *pp(IOB_MS+1,K,JOB_MS) ) + &
893	DYOB_MS* ( (1.-DXOB_MS)*pp(IOB_MS,K,JOB_MS+1) + &
894	DXOB_MS *pp(IOB_MS+1,K,JOB_MS+1) ) )
895	enddo
896
897	!ajb 20040119: Note based on bugfix for dot/cross points split across processors,
898	!ajb which was actually a serial code fix: The ityp=2 (v-points) and
899	!ajb itype=3 (mass-points) cases should not use the ityp=1 (u-points)
900	!ajb case) rko! This is necessary for bit-for-bit reproducability
901	!ajb with parallel run. (fine mesh)
902
903	if(abs(rko(n)+99).lt.1.)then
904	pob = varobs(5,n)
905
906	if(pob .eq.-888888.) then
907	hob = varobs(6,n)
908	if(hob .gt. -800000. ) then
909	pob = ht_to_p( hob, ppbo, pbbo, z, iob_ms, job_ms, &
910	dxob_ms, dyob_ms, k_start, k_end, kds,kde, &
911	ims,ime, jms,jme, kms,kme )
912	endif
913	endif
914
915	if(pob .gt.-800000.)then
916	do k=k_end-1,1,-1
917	kbot = k
918	if(pob .le. pbbo(k)+ppbo(k)) then
919	goto 198
920	endif
921	enddo
922	198 continue
923
924	pphi = ppbo(kbot+1)
925	pbhi = pbbo(kbot+1)
926
927	rko(n) = real(kbot+1)- &
928	( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
929	rko(n)=max(rko(n),1.0)
930	endif
931	endif
932
933	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
934	ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
935	#endif
936
937	! ajb 20090423: If obs-in-height, varobs(5,n) is sum of pob (see note above).
938	if(varobs(5,n) .eq. -888888. .and. varobs(6,n) .gt. -800000.) then
939	varobs(5,n) = wrf_dm_sum_real ( pob )
940	endif
941
942	RC(N)=RKO(N)
943
944	ENDDO ! END FINE MESH LOOP OVER NSTA
945
946	ENDIF ! end if(inest.eq.1)
947
948	! INITIALIZE THE ARRAY OF DIFFERENCES BETWEEN THE OBSERVATIONS
949	! AND THE LOCAL FORECAST VALUES TO ZERO. FOR THE FINE MESH
950	! ONLY, SET THE DIFFERENCE TO A LARGE DUMMY VALUE IF THE
951	! OBSERVATION IS OUTSIDE THE FINE MESH DOMAIN.
952
953	! SET DIFFERENCE VALUE TO A DUMMY VALUE FOR OBS POINTS OUTSIDE
954	! CURRENT DOMAIN
955	IF(ITYP.EQ.1) THEN
956	NLB=1
957	NLE=1
958	ELSE IF(ITYP.EQ.2) THEN
959	NLB=2
960	NLE=2
961	ELSE
962	NLB=3
963	NLE=5
964	ENDIF
965	DO IVAR=NLB,NLE
966	DO N=1,NSTA
967	IF((RA(N)-1.).LT.0)THEN
968	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
969	ENDIF
970	IF((RB(N)-1.).LT.0)THEN
971	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
972	ENDIF
973	IF((FLOAT(ide)-2.0*GRIDX-RA(N)-1.E-10).LT.0)THEN
974	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
975	ENDIF
976	IF((FLOAT(jde)-2.0*GRIDY-RB(N)-1.E-10).LT.0)THEN
977	ERRF(IVAR,N)=ERRF(IVAR,N)+DUMMY
978	ENDIF
979	if(rc(n).lt.1.)errf(ivar,n)=errf(ivar,n)+dummy
980	ENDDO
981	ENDDO
982
983	! NOW FIND THE EXACT OFFSET OF EACH OBSERVATION FROM THE
984	! GRID POINT TOWARD THE LOWER LEFT
985	DO N=1,NSTA
986	IA(N)=RA(N)
987	IB(N)=RB(N)
988	IC(N)=RC(N)
989	ENDDO
990	DO N=1,NSTA
991	RA(N)=RA(N)-FLOAT(IA(N))
992	RB(N)=RB(N)-FLOAT(IB(N))
993	RC(N)=RC(N)-FLOAT(IC(N))
994	ENDDO
995	! PERFORM A TRILINEAR EIGHT-POINT (3-D) INTERPOLATION
996	! TO FIND THE FORECAST VALUE AT THE EXACT OBSERVATION
997	! POINTS FOR U, V, T, AND Q.
998
999	! Compute local masks for dot and cross points.
1000	if(ityp.eq.1) then
1001	DO N=1,NSTA
1002	IOB=MAX0(1,IA(N))
1003	IOB=MIN0(IOB,ide-1)
1004	JOB=MAX0(1,IB(N))
1005	JOB=MIN0(JOB,jde-1)
1006	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1007	! Set mask for U-momemtum points to be handled by this processor
1008	MP_LOCAL_UOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
1009	#endif
1010	ENDDO
1011	endif
1012	if(ityp.eq.2) then
1013	DO N=1,NSTA
1014	IOB=MAX0(1,IA(N))
1015	IOB=MIN0(IOB,ide-1)
1016	JOB=MAX0(1,IB(N))
1017	JOB=MIN0(JOB,jde-1)
1018	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1019	! Set mask for V-momentum points to be handled by this processor
1020	MP_LOCAL_VOBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
1021	#endif
1022	ENDDO
1023	endif
1024	if(ityp.eq.3) then
1025	DO N=1,NSTA
1026	IOB=MAX0(1,IA(N))
1027	IOB=MIN0(IOB,ide-1)
1028	JOB=MAX0(1,IB(N))
1029	JOB=MIN0(JOB,jde-1)
1030	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1031	! Set mask for cross (mass) points to be handled by this processor
1032	MP_LOCAL_COBMASK(N) = TILE_MASK(IOB, JOB, its, ite, jts, jte)
1033	#endif
1034	ENDDO
1035	endif
1036
1037	!**********************************************************
1038	! PROCESS U VARIABLE (IVAR=1)
1039	!**********************************************************
1040	IF(ITYP.EQ.1) THEN
1041	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1042	DO N=1,NSTA
1043	IF(MP_LOCAL_UOBMASK(N)) THEN
1044	ERRF(9,N)=rko(n) !RKO is needed by neighboring processors !ajb
1045	ENDIF
1046	ENDDO
1047	#endif
1048	IF(ISWIND.EQ.1) THEN
1049	DO N=1,NSTA
1050	IOB=MAX0(2,IA(N))
1051	IOB=MIN0(IOB,ide-1)
1052	IOBM=MAX0(1,IOB-1)
1053	IOBP=MIN0(ide-1,IOB+1)
1054	JOB=MAX0(2,IB(N))
1055	JOB=MIN0(JOB,jde-1)
1056	JOBM=MAX0(1,JOB-1)
1057	JOBP=MIN0(jde-1,JOB+1)
1058	KOB=MIN0(K_END,IC(N))
1059
1060	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1061	IF(MP_LOCAL_UOBMASK(N))THEN ! Do if obs on this processor
1062	#endif
1063	KOBP=MIN0(KOB+1,k_end)
1064	DXOB=RA(N)
1065	DYOB=RB(N)
1066	DZOB=RC(N)
1067
1068	! Compute surface pressure values at surrounding U and V points
1069	PUG1 = .5*( pbase(IOBM,1,JOB) + pbase(IOB,1,JOB) )
1070	PUG2 = .5*( pbase(IOB,1,JOB) + pbase(IOBP,1,JOB) )
1071
1072	! This is to correct obs to model sigma level using reverse similarity theory
1073	if(rko(n).eq.1.0)then
1074	uratiob=((1.-DYOB)((1.-DXOB)uratio(IOB,JOB)+ &
1075	DXOB*uratio(IOBP,JOB) &
1076	)+DYOB((1.-DXOB)uratio(IOB,JOBP)+ &
1077	DXOB*uratio(IOBP,JOBP)))
1078	else
1079	uratiob=1.
1080	endif
1081	!YLIU Some PBL scheme do not define the vratio/uratio
1082	if(abs(uratiob).lt.1.0e-3) then
1083	uratiob=1.
1084	endif
1085
1086	! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
1087	! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
1088	! OUTSIDE THE CURRENT DOMAIN
1089
1090	! U COMPONENT WIND ERROR
1091	ERRF(1,N)=ERRF(1,N)+uratiobVAROBS(1,N)-((1.-DZOB) &
1092	((1.-DyOB)*((1.- &
1093	DxOB)UB(IOB,KOB,JOB)+DxOBUB(IOB+1,KOB,JOB) &
1094	)+DyOB((1.-DxOB)UB(IOB,KOB,JOB+1)+DxOB* &
1095	UB(IOB+1,KOB,JOB+1)))+DZOB((1.-DyOB)((1.-DxOB) &
1096	UB(IOB,KOBP,JOB)+DxOBUB(IOB+1,KOBP,JOB))+ &
1097	DyOB((1.-DxOB)UB(IOB,KOBP,JOB+1)+DxOB* &
1098	UB(IOB+1,KOBP,JOB+1))))
1099
1100	! if(n.le.10) then
1101	! write(6,*)
1102	! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF1 at ',iob,job,kob, &
1103	! ' N = ',n,' inest = ',inest
1104	! write(6,*) 'VAROBS(1,N) = ',varobs(1,n)
1105	! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
1106	! write(6,*) 'UB(IOB,KOB,JOB) = ',UB(IOB,KOB,JOB)
1107	! write(6,*) 'UB(IOB+1,KOB,JOB) = ',UB(IOB+1,KOB,JOB)
1108	! write(6,*) 'UB(IOB,KOB,JOB+1) = ',UB(IOB,KOB,JOB+1)
1109	! write(6,*) 'UB(IOB+1,KOB,JOB+1) = ',UB(IOB+1,KOB,JOB+1)
1110	! write(6,*) 'UB(IOB,KOBP,JOB) = ',UB(IOB,KOBP,JOB)
1111	! write(6,*) 'UB(IOB+1,KOBP,JOB) = ',UB(IOB+1,KOBP,JOB)
1112	! write(6,*) 'UB(IOB,KOBP,JOB+1) = ',UB(IOB,KOBP,JOB+1)
1113	! write(6,*) 'UB(IOB+1,KOBP,JOB+1) = ',UB(IOB+1,KOBP,JOB+1)
1114	! write(6,*) 'uratiob = ',uratiob
1115	! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
1116	! write(6,*) 'ERRF(1,N) = ',errf(1,n)
1117	! write(6,*)
1118	! endif
1119
1120
1121	! Store model surface pressure (not the error!) at U point.
1122	ERRF(7,N)=.001( (1.-DXOB)PUG1 + DXOB*PUG2 )
1123
1124	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1125	ENDIF ! end IF( MP_LOCAL_UOBMASK(N) )
1126	#endif
1127	ENDDO ! END U-POINT LOOP OVER OBS
1128
1129	ENDIF ! end if(iswind.eq.1)
1130
1131	ENDIF ! ITYP=1: PROCESS U
1132
1133	!**********************************************************
1134	! PROCESS V VARIABLE (IVAR=2)
1135	!**********************************************************
1136	IF(ITYP.EQ.2) THEN
1137
1138	IF(ISWIND.EQ.1) THEN
1139	DO N=1,NSTA
1140	IOB=MAX0(2,IA(N))
1141	IOB=MIN0(IOB,ide-1)
1142	IOBM=MAX0(1,IOB-1)
1143	IOBP=MIN0(ide-1,IOB+1)
1144	JOB=MAX0(2,IB(N))
1145	JOB=MIN0(JOB,jde-1)
1146	JOBM=MAX0(1,JOB-1)
1147	JOBP=MIN0(jde-1,JOB+1)
1148	KOB=MIN0(K_END,IC(N))
1149
1150	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1151	IF(MP_LOCAL_VOBMASK(N))THEN ! Do if obs on this processor
1152	#endif
1153	KOBP=MIN0(KOB+1,k_end)
1154	DXOB=RA(N)
1155	DYOB=RB(N)
1156	DZOB=RC(N)
1157
1158	! Compute surface pressure values at surrounding U and V points
1159	PVG1 = .5*( pbase(IOB,1,JOBM) + pbase(IOB,1,JOB) )
1160	PVG2 = .5*( pbase(IOB,1,JOB) + pbase(IOB,1,JOBP) )
1161
1162	! This is to correct obs to model sigma level using reverse similarity theory
1163	if(rko(n).eq.1.0)then
1164	vratiob=((1.-DYOB)((1.-DXOB)vratio(IOB,JOB)+ &
1165	DXOB*vratio(IOBP,JOB) &
1166	)+DYOB((1.-DXOB)vratio(IOB,JOBP)+ &
1167	DXOB*vratio(IOBP,JOBP)))
1168	else
1169	vratiob=1.
1170	endif
1171	!YLIU Some PBL scheme do not define the vratio/uratio
1172	if(abs(vratiob).lt.1.0e-3) then
1173	vratiob=1.
1174	endif
1175
1176	! INITIAL ERRF(IVAR,N) IS ZERO FOR OBSERVATIONS POINTS
1177	! WITHIN THE DOMAIN, AND A LARGE DUMMY VALUE FOR POINTS
1178	! OUTSIDE THE CURRENT DOMAIN
1179
1180	! V COMPONENT WIND ERROR
1181	ERRF(2,N)=ERRF(2,N)+vratiobVAROBS(2,N)-((1.-DZOB) &
1182	((1.-DyOB)*((1.- &
1183	DxOB)VB(IOB,KOB,JOB)+DxOBVB(IOB+1,KOB,JOB) &
1184	)+DyOB((1.-DxOB)VB(IOB,KOB,JOB+1)+DxOB* &
1185	VB(IOB+1,KOB,JOB+1)))+DZOB((1.-DyOB)((1.-DxOB) &
1186	VB(IOB,KOBP,JOB)+DxOBVB(IOB+1,KOBP,JOB))+ &
1187	DyOB((1.-DxOB)VB(IOB,KOBP,JOB+1)+DxOB* &
1188	VB(IOB+1,KOBP,JOB+1))))
1189
1190	! if(n.le.10) then
1191	! write(6,*)
1192	! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF2 at ',iob,job,kob, &
1193	! ' N = ',n,' inest = ',inest
1194	! write(6,*) 'VAROBS(2,N) = ',varobs(2,n)
1195	! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
1196	! write(6,*) 'VB(IOB,KOB,JOB) = ',VB(IOB,KOB,JOB)
1197	! write(6,*) 'ERRF(2,N) = ',errf(2,n)
1198	! write(6,*) 'vratiob = ',vratiob
1199	! write(6,*)
1200	! endif
1201
1202
1203	! Store model surface pressure (not the error!) at V point.
1204	ERRF(8,N)=.001( (1.-DYOB)PVG1 + DYOB*PVG2 )
1205
1206	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1207	ENDIF ! end IF( MP_LOCAL_VOBMASK(N) )
1208	#endif
1209	ENDDO ! END V-POINT LOOP OVER OBS
1210
1211	ENDIF ! end if(iswind.eq.1)
1212
1213	ENDIF ! ITYP=1: PROCESS V
1214
1215	!**********************************************************
1216	! PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
1217	!**********************************************************
1218	IF(ITYP.EQ.3) THEN
1219
1220	IF(ISTEMP.EQ.1 .OR. ISMOIS.EQ.1) THEN
1221	DO N=1,NSTA
1222	IOB=MAX0(1,IA(N))
1223	IOB=MIN0(IOB,ide-1)
1224	JOB=MAX0(1,IB(N))
1225	JOB=MIN0(JOB,jde-1)
1226	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1227	IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
1228	#endif
1229	KOB=MIN0(k_end,IC(N))
1230	KOBP=MIN0(KOB+1,K_END)
1231	DXOB=RA(N)
1232	DYOB=RB(N)
1233	DZOB=RC(N)
1234
1235	! This is to correct obs to model sigma level using reverse similarity theory
1236	if(rko(n).eq.1.0)then
1237	tratxob=((1.-DYOB)((1.-DXOB)tratx(IOB,JOB)+ &
1238	DXOB*tratx(IOB+1,JOB) &
1239	)+DYOB((1.-DXOB)tratx(IOB,JOB+1)+ &
1240	DXOB*tratx(IOB+1,JOB+1)))
1241	else
1242	tratxob=1.
1243	endif
1244
1245	!yliu
1246	if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
1247
1248	!ajb We must convert temperature to potential temperature
1249	obs_pottemp = -888888.
1250	if(varobs(3,n).gt.-800000. .and. varobs(5,n).gt.-800000) then
1251	obs_pottemp = varobs(3,n)(100./varobs(5,n))*RCP - t0
1252	endif
1253
1254	ERRF(3,N)=ERRF(3,N)+tratxobobs_pottemp-((1.-DZOB) &
1255	((1.-DyOB)*((1.- &
1256	DxOB)(TB(IOB,KOB,JOB))+DxOB &
1257	(TB(IOB+1,KOB,JOB)))+DyOB((1.-DxOB) &
1258	(TB(IOB,KOB,JOB+1))+DxOB* &
1259	(TB(IOB+1,KOB,JOB+1))))+DZOB*((1.- &
1260	DyOB)((1.-DxOB)(TB(IOB,KOBP,JOB))+DxOB* &
1261	(TB(IOB+1,KOBP,JOB)))+DyOB((1.-DxOB) &
1262	(TB(IOB,KOBP,JOB+1))+DxOB* &
1263	(TB(IOB+1,KOBP,JOB+1)))))
1264
1265	! if(n.le.10) then
1266	! write(6,*)
1267	! write(6,'(a,i3,i3,i3,a,i3,a,i2)') 'ERRF3 at ',iob,job,kob, &
1268	! ' N = ',n,' inest = ',inest
1269	! write(6,*) 'VAROBS(3,N) = ',varobs(3,n)
1270	! write(6,*) 'VAROBS(5,N) = ',varobs(5,n)
1271	! write(6,*) 'TB(IOB,KOB,JOB) = ',TB(IOB,KOB,JOB)
1272	! write(6,*) 'TB(IOB+1,KOB,JOB) = ',TB(IOB+1,KOB,JOB)
1273	! write(6,*) 'TB(IOB,KOB,JOB+1) = ',TB(IOB,KOB,JOB+1)
1274	! write(6,*) 'TB(IOB+1,KOB,JOB+1) = ',TB(IOB+1,KOB,JOB+1)
1275	! write(6,*) 'TB(IOB,KOBP,JOB) = ',TB(IOB,KOBP,JOB)
1276	! write(6,*) 'TB(IOB+1,KOBP,JOB) = ',TB(IOB+1,KOBP,JOB)
1277	! write(6,*) 'TB(IOB,KOBP,JOB+1) = ',TB(IOB,KOBP,JOB+1)
1278	! write(6,*) 'TB(IOB+1,KOBP,JOB+1) = ',TB(IOB+1,KOBP,JOB+1)
1279	! write(6,*) 'tratxob = ',tratxob
1280	! write(6,*) 'DXOB,DYOB,DZOB = ',DXOB,DYOB,DZOB
1281	! write(6,*) 'ERRF(3,N) = ',errf(3,n)
1282	! write(6,*)
1283	! endif
1284
1285	! MOISTURE ERROR
1286	ERRF(4,N)=ERRF(4,N)+VAROBS(4,N)-((1.-DZOB)((1.-DyOB)((1.- &
1287	DxOB)QVB(IOB,KOB,JOB)+DxOB &
1288	QVB(IOB+1,KOB,JOB))+DyOB((1.-DxOB) &
1289	QVB(IOB,KOB,JOB+1)+DxOB* &
1290	QVB(IOB+1,KOB,JOB+1)))+DZOB*((1.- &
1291	DyOB)((1.-DxOB)QVB(IOB,KOBP,JOB)+DxOB &
1292	QVB(IOB+1,KOBP,JOB))+DyOB((1.-DxOB &
1293	)QVB(IOB,KOBP,JOB+1)+DxOB &
1294	QVB(IOB+1,KOBP,JOB+1))))
1295
1296	! Store model surface pressure (not the error!) at T-point
1297	ERRF(6,N)= .001* &
1298	((1.-DyOB)((1.-DxOB)pbase(IOB,1,JOB)+DxOB* &
1299	pbase(IOB+1,1,JOB))+DyOB((1.-DxOB) &
1300	pbase(IOB,1,JOB+1)+DxOB*pbase(IOB+1,1,JOB+1) ))
1301
1302	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1303	ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
1304	#endif
1305	ENDDO ! END T and QV LOOP OVER OBS
1306
1307	ENDIF !end if(istemp.eq.1 .or. ismois.eq.1)
1308
1309	!*******************************************************
1310	! USE ERRF(5,N) TO STORE KPBL AT (I,J) NEAREST THE OBS
1311	!*******************************************************
1312	DO N=1,NSTA
1313	IOB=MAX0(1,IA(N))
1314	IOB=MIN0(IOB,ide-1)
1315	JOB=MAX0(1,IB(N))
1316	JOB=MIN0(JOB,jde-1)
1317	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1318	IF(MP_LOCAL_COBMASK(N)) THEN ! Do if obs on this processor
1319	#endif
1320	DXOB=RA(N)
1321	DYOB=RB(N)
1322	ERRF(5,N) = kpbl(iob+nint(dxob),job+nint(dyob))
1323
1324	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1325	ENDIF ! end IF( MP_LOCAL_COBMASK(N) )
1326	#endif
1327	ENDDO
1328	!!**********************************************************
1329	ENDIF ! end if(ityp.eq.3)
1330
1331	ENDDO ! END BIG LOOP
1332
1333	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1334	! Gather the errf values calculated by the processors with the obs.
1335	CALL get_full_obs_vector(nsta, nerrf, niobf, mp_local_uobmask, &
1336	mp_local_vobmask, mp_local_cobmask, errf)
1337
1338	! 02252010: Go ahead and assign rko for "obs-off" processors here, to
1339	! fix the problem where duplicate obs can be dropped from
1340	! the "obs-on" processor, but not from the others, due to
1341	! rko being -99 on the "obs-off" processors.
1342	do n=1,nsta
1343	rko(n) = errf(9,n)
1344	enddo
1345	! 02252010: End bugfix.
1346	#endif
1347
1348	! Print obs information.
1349	CALL print_obs_info(iprt,inest,niobf,rio,rjo,rko, &
1350	prt_max,prt_freq,obs_prt,stnid_prt, &
1351	lat_prt,lon_prt,mlat_prt,mlon_prt, &
1352	timeob,xtime)
1353
1354	! DIFFERENCE BETWEEN OBS AND FCST IS COMPLETED
1355	IF(INEST.EQ.1)THEN
1356	INPF=NPFI
1357	ELSE
1358	FNPF=IRATIO**LEVIDN(INEST)
1359	INPF=FNPF*NPFI
1360	ENDIF
1361	! Gross error check for temperature. Set all vars bad.
1362	do n=1,nsta
1363	if((abs(errf(3,n)).gt.20.).and. &
1364	(errf(3,n).gt.-800000.))then
1365
1366	errf(1,n)=-888888.
1367	errf(2,n)=-888888.
1368	errf(3,n)=-888888.
1369	errf(4,n)=-888888.
1370	varobs(1,n)=-888888.
1371	varobs(2,n)=-888888.
1372	varobs(3,n)=-888888.
1373	varobs(4,n)=-888888.
1374	endif
1375	enddo
1376
1377	! For printout
1378	! IF(MOD(KTAU,INPF).NE.0) THEN
1379	! RETURN
1380	! ENDIF
1381
1382	RETURN
1383	END SUBROUTINE errob
1384
1385	SUBROUTINE upoint(i_start,i_end, j_start,j_end, ids,ide, ims,ime, jms,jme, &
1386	arrin, arrout)
1387	!------------------------------------------------------------------------------
1388	! PURPOSE: This subroutine interpolates a real 2D array defined over mass
1389	! coordinate points, to U (momentum) points.
1390	!
1391	!------------------------------------------------------------------------------
1392	IMPLICIT NONE
1393
1394	INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
1395	INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
1396	INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
1397	INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
1398	INTEGER, INTENT(IN) :: ids ! Starting i index for entire model domain
1399	INTEGER, INTENT(IN) :: ide ! Ending i index for entire model domain
1400	INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
1401	INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
1402	INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
1403	INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
1404	REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
1405	REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on U points
1406
1407	! Local variables
1408	integer :: i, j
1409
1410	! Do domain interior first
1411	do j = j_start, j_end
1412	do i = max(2,i_start), i_end
1413	arrout(i,j) = 0.5*(arrin(i,j)+arrin(i-1,j))
1414	enddo
1415	enddo
1416
1417	! Do west-east boundaries
1418	if(i_start .eq. ids) then
1419	do j = j_start, j_end
1420	arrout(i_start,j) = arrin(i_start,j)
1421	enddo
1422	endif
1423	if(i_end .eq. ide-1) then
1424	do j = j_start, j_end
1425	arrout(i_end+1,j) = arrin(i_end,j)
1426	enddo
1427	endif
1428
1429	RETURN
1430	END SUBROUTINE upoint
1431
1432	SUBROUTINE vpoint(i_start,i_end, j_start,j_end, jds,jde, ims,ime, jms,jme, &
1433	arrin, arrout)
1434	!------------------------------------------------------------------------------
1435	! PURPOSE: This subroutine interpolates a real 2D array defined over mass
1436	! coordinate points, to V (momentum) points.
1437	!
1438	!------------------------------------------------------------------------------
1439	IMPLICIT NONE
1440
1441	INTEGER, INTENT(IN) :: i_start ! Starting i index for this model tile
1442	INTEGER, INTENT(IN) :: i_end ! Ending i index for this model tile
1443	INTEGER, INTENT(IN) :: j_start ! Starting j index for this model tile
1444	INTEGER, INTENT(IN) :: j_end ! Ending j index for this model tile
1445	INTEGER, INTENT(IN) :: jds ! Starting j index for entire model domain
1446	INTEGER, INTENT(IN) :: jde ! Ending j index for entire model domain
1447	INTEGER, INTENT(IN) :: ims ! Starting i index for model patch
1448	INTEGER, INTENT(IN) :: ime ! Ending i index for model patch
1449	INTEGER, INTENT(IN) :: jms ! Starting j index for model patch
1450	INTEGER, INTENT(IN) :: jme ! Ending j index for model patch
1451	REAL, INTENT(IN) :: arrin ( ims:ime, jms:jme ) ! input array on mass points
1452	REAL, INTENT(OUT) :: arrout( ims:ime, jms:jme ) ! output array on V points
1453
1454	! Local variables
1455	integer :: i, j
1456
1457	! Do domain interior first
1458	do j = max(2,j_start), j_end
1459	do i = i_start, i_end
1460	arrout(i,j) = 0.5*(arrin(i,j)+arrin(i,j-1))
1461	enddo
1462	enddo
1463
1464	! Do south-north boundaries
1465	if(j_start .eq. jds) then
1466	do i = i_start, i_end
1467	arrout(i,j_start) = arrin(i,j_start)
1468	enddo
1469	endif
1470	if(j_end .eq. jde-1) then
1471	do i = i_start, i_end
1472	arrout(i,j_end+1) = arrin(i,j_end)
1473	enddo
1474	endif
1475
1476	RETURN
1477	END SUBROUTINE vpoint
1478
1479	LOGICAL FUNCTION TILE_MASK(iloc, jloc, its, ite, jts, jte)
1480	!------------------------------------------------------------------------------
1481	! PURPOSE: Check to see if an i, j grid coordinate is in the tile index range.
1482	!
1483	! Returns: TRUE if the grid coordinate (ILOC,JLOC) is in the tile defined by
1484	! tile-range indices (its,jts) and (ite,jte)
1485	! FALSE otherwise.
1486	!
1487	!------------------------------------------------------------------------------
1488	IMPLICIT NONE
1489
1490	INTEGER, INTENT(IN) :: iloc
1491	INTEGER, INTENT(IN) :: jloc
1492	INTEGER, INTENT(IN) :: its
1493	INTEGER, INTENT(IN) :: ite
1494	INTEGER, INTENT(IN) :: jts
1495	INTEGER, INTENT(IN) :: jte
1496
1497	! Local variables
1498	LOGICAL :: retval
1499
1500	TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND. &
1501	jloc .LE. jte .AND. jloc .GE. jts )
1502
1503	RETURN
1504	END FUNCTION TILE_MASK
1505
1506	!-----------------------------------------------------------------------
1507	SUBROUTINE nudob(j, ivar, aten, inest, ifrest, ktau, ktaur, &
1508	xtime, mu, msfx, msfy, nndgv, nerrf, niobf, maxdom, &
1509	npfi, ionf, rinxy, twindo, &
1510	nudge_pbl, &
1511	sfcfact, sfcfacr, &
1512	levidn, &
1513	parid, nstat, &
1514	rinfmn, rinfmx, pfree, dcon, tfaci, &
1515	sfc_scheme_horiz, sfc_scheme_vert, maxsnd_gap, &
1516	lev_in_ob, plfo, nlevs_ob, &
1517	iratio, dx, dtmin, rio, rjo, rko, &
1518	timeob, varobs, errf, pbase, ptop, pp, &
1519	iswind, istemp, ismois, giv, git, giq, &
1520	savwt, kpblt, nscan, &
1521	vih1, vih2, terrh, zslab, &
1522	iprt, &
1523	ids,ide, jds,jde, kds,kde, & ! domain dims
1524	ims,ime, jms,jme, kms,kme, & ! memory dims
1525	its,ite, jts,jte, kts,kte ) ! tile dims
1526
1527	!-----------------------------------------------------------------------
1528	USE module_model_constants
1529	USE module_domain
1530	!-----------------------------------------------------------------------
1531	IMPLICIT NONE
1532	!-----------------------------------------------------------------------
1533	!
1534	! PURPOSE: THIS SUBROUTINE GENERATES NUDGING TENDENCIES FOR THE J-TH
1535	! VERTICAL SLICE (I-K PLANE) FOR FOUR-DIMENSIONAL DATA
1536	! ASSIMILATION FROM INDIVIDUAL OBSERVATIONS. THE NUDGING
1537	! TENDENCIES ARE FOUND FROM A ONE-PASS CALCULATION OF
1538	! WEIGHTING FACTORS SIMILAR TO THE BENJAMIN-SEAMAN OBJECTIVE
1539	! ANALYSIS. THIS SUBROUTINE IS DESIGNED FOR RAPID EXECUTION
1540	! AND MINIMAL STORAGE REQUIREMENTS. ALGORITHMS SHOULD BE
1541	! VECTORIZED WHEREVER POSSIBLE.
1542	!
1543	! HISTORY: Original author: MM5 version???
1544	! 02/04/2004 - Creation of WRF version. Al Bourgeois
1545	! 08/28/2006 - Conversion from F77 to F90 Al Bourgeois
1546	!------------------------------------------------------------------------------
1547	!
1548	! NOTE: This routine was originally designed for MM5, which uses
1549	! a nonstandard (I,J) coordinate system. For WRF, I is the
1550	! east-west running coordinate, and J is the south-north
1551	! running coordinate. So a "J-slab" here is west-east in
1552	! extent, not south-north as for MM5. -ajb 06/10/2004
1553	!
1554	! NET WEIGHTING (WT) OF THE DIFFERENCE BETWEEN THE OBSERVATIONS
1555	! AND LOCAL FORECAST VALUES IS BASED ON THE MULTIPLE OF THREE
1556	! TYPES OF FACTORS:
1557	! 1) TIME WEIGHTING - ONLY OBSERVATIONS WITHIN A SELECTED
1558	! TIME WINDOW (TWINDO) CENTERED AT THE CURRENT FORECAST
1559	! TIME (XTIME) ARE USED. OBSERVATIONS CLOSEST TO
1560	! XTIME ARE TIME-WEIGHTED MOST HEAVILY (TIMEWT)
1561	! 2) VERTICAL WEIGHTING - NON-ZERO WEIGHTS (WTSIG) ARE
1562	! CALCULATED WITHIN A VERTICAL REGION OF INFLUENCE
1563	! (RINSIG).
1564	! 3) HORIZONTAL WEIGHTING - NON-ZERO WEIGHTS (WTIJ) ARE
1565	! CALCULATED WITHIN A RADIUS OF INFLUENCE (RINXY). THE
1566	! VALUE OF RIN IS DEFINED IN KILOMETERS, AND CONVERTED
1567	! TO GRID LENGTHS FOR THE APPROPRIATE MESH SIZE.
1568	!
1569	! THE FIVE FORECAST VARIABLES ARE PROCESSED BY CHANGING THE
1570	! VALUE OF IVAR AS FOLLOWS:
1571	! IVAR VARIABLE(TAU-1)
1572	! ---- ---------------
1573	! 1 U
1574	! 2 V
1575	! 3 T
1576	! 4 QV
1577	! 5 PSB(CROSS) REMOVED IN V3
1578	! (6) PSB(DOT)
1579	!
1580	!-----------------------------------------------------------------------
1581	!
1582	! Description of input arguments.
1583	!
1584	!-----------------------------------------------------------------------
1585
1586	INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde ! domain dims.
1587	INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
1588	INTEGER, INTENT(IN) :: its,ite, jts,jte, kts,kte ! tile dims.
1589	INTEGER, INTENT(IN) :: j ! south-north running coordinate.
1590	INTEGER, INTENT(IN) :: ivar
1591	INTEGER, INTENT(IN) :: inest ! domain index
1592	LOGICAL, INTENT(IN) :: ifrest
1593	INTEGER, INTENT(IN) :: ktau
1594	INTEGER, INTENT(IN) :: ktaur
1595	REAL, INTENT(IN) :: xtime ! forecast time in minutes
1596	INTEGER, INTENT(IN) :: nndgv ! number of nudge variables
1597	INTEGER, INTENT(IN) :: nerrf ! number of error fields
1598	INTEGER, INTENT(IN) :: niobf ! number of observations
1599	INTEGER, INTENT(IN) :: maxdom ! maximum number of domains
1600	INTEGER, INTENT(IN) :: npfi
1601	INTEGER, INTENT(IN) :: ionf
1602	REAL, INTENT(IN) :: rinxy
1603	REAL, INTENT(IN) :: twindo
1604	REAL, intent(in) :: sfcfact ! scale for time window (surface obs)
1605	REAL, intent(in) :: sfcfacr ! scale for hor rad inf (surface obs)
1606	LOGICAL, intent(in) :: nudge_pbl ! flag for nudging in pbl
1607	INTEGER, INTENT(IN) :: levidn(maxdom) ! level of nest
1608	INTEGER, INTENT(IN) :: parid(maxdom) ! parent domain id
1609	INTEGER, INTENT(IN) :: nstat ! number of obs stations
1610	REAL, INTENT(IN) :: rinfmn ! minimum radius of influence
1611	REAL, INTENT(IN) :: rinfmx ! maximum radius of influence
1612	REAL, INTENT(IN) :: pfree ! pressure level (cb) where terrain effect becomes small
1613	REAL, INTENT(IN) :: dcon ! 1/DPSMX
1614	REAL, INTENT(IN) :: tfaci ! scale factor used for ramp-down in dynamic initialization
1615	INTEGER , INTENT(IN) :: sfc_scheme_horiz ! horizontal spreading scheme for surf obs (wrf or orig mm5)
1616	INTEGER , INTENT(IN) :: sfc_scheme_vert ! vertical spreading scheme for surf obs (orig or regime vif)
1617	REAL , INTENT(IN) :: maxsnd_gap ! max allowed pressure gap in soundings, for interp (centibars)
1618	REAL, INTENT(IN) :: lev_in_ob(niobf) ! Level in sounding-type obs.
1619	REAL, intent(IN) :: plfo(niobf) ! Index for type of obs platform
1620	REAL, INTENT(IN) :: nlevs_ob(niobf) ! Number of levels in sounding.
1621	INTEGER, INTENT(IN) :: iratio ! Nest to parent gridsize ratio.
1622	REAL, INTENT(IN) :: dx ! This domain grid cell-size (m)
1623	REAL, INTENT(IN) :: dtmin ! Model time step in minutes
1624	REAL, INTENT(IN) :: rio(niobf) ! Obs west-east coordinate (non-stag grid).
1625	REAL, INTENT(IN) :: rjo(niobf) ! Obs south-north coordinate (non-stag grid).
1626	REAL, INTENT(INOUT) :: rko(niobf) ! Obs vertical coordinate.
1627	REAL, INTENT(IN) :: timeob(niobf)
1628	REAL, INTENT(IN) :: varobs(nndgv,niobf)
1629	REAL, INTENT(IN) :: errf(nerrf, niobf)
1630	REAL, INTENT(IN) :: pbase( ims:ime, kms:kme ) ! Base pressure.
1631	REAL, INTENT(IN) :: ptop
1632	REAL, INTENT(IN) :: pp( ims:ime, kms:kme ) ! Pressure perturbation (Pa)
1633	REAL, INTENT(IN) :: mu(ims:ime) ! Air mass on u, v, or mass-grid
1634	REAL, INTENT(IN) :: msfx(ims:ime) ! Map scale (only used for vars u & v)
1635	REAL, INTENT(IN) :: msfy(ims:ime) ! Map scale (only used for vars u & v)
1636	INTEGER, intent(in) :: iswind ! Nudge flag for wind
1637	INTEGER, intent(in) :: istemp ! Nudge flag for temperature
1638	INTEGER, intent(in) :: ismois ! Nudge flag for moisture
1639	REAL, intent(in) :: giv ! Coefficient for wind
1640	REAL, intent(in) :: git ! Coefficient for temperature
1641	REAL, intent(in) :: giq ! Coefficient for moisture
1642	REAL, INTENT(INOUT) :: aten( ims:ime, kms:kme)
1643	REAL, INTENT(INOUT) :: savwt( nndgv, ims:ime, kms:kme )
1644	INTEGER, INTENT(IN) :: kpblt(ims:ime)
1645	INTEGER, INTENT(IN) :: nscan ! number of scans
1646	REAL, INTENT(IN) :: vih1(its:ite) ! Vert infl ht (m) abv grd for full wts
1647	REAL, INTENT(IN) :: vih2(its:ite) ! Vert infl ht (m) abv grd for ramp
1648	REAL, INTENT(IN) :: terrh(ims:ime) ! Terrain height (m)
1649	! INTEGER, INTENT(IN) :: vik1(its:ite) ! Vertical infl k-level for full wts
1650	! INTEGER, INTENT(IN) :: vik2(its:ite) ! Vertical infl k-level for ramp
1651	REAL, INTENT(IN) :: zslab(ims:ime, kms:kme) ! model ht above ground (m)
1652	LOGICAL, INTENT(IN) :: iprt ! print flag
1653
1654	! Local variables
1655	integer :: mm(maxdom)
1656	integer :: kobs ! k-lev below obs (for obs straddling pblt)
1657	integer :: kpbl_obs(nstat) ! kpbl at grid point (IOB,JOB) (ajb 20090519)
1658	real :: ra(niobf)
1659	real :: rb(niobf)
1660	real :: psurf(niobf)
1661	real :: wtsig(kms:kme),wt(ims:ime,kms:kme),wt2err(ims:ime,kms:kme)
1662	real :: rscale(ims:ime) ! For converting to rho-coupled units.
1663	real :: wtij(ims:ime) ! For holding weights in i-loop
1664	real :: reserf(100)
1665	character*40 name
1666	character*3 chr_hr
1667	character(len=200) :: msg ! Argument to wrf_message
1668
1669	!*** DECLARATIONS FOR IMPLICIT NONE
1670	integer :: i,k,iplo,icut,ipl,inpf,infr,jjjn
1671	integer :: igrid,n,nml,nnl,nsthis,nsmetar,nsspeci,nsship
1672	integer :: nssynop,nstemp,nspilot,nssatwnds,nssams,nsprofs
1673	integer :: maxi,mini,maxj,minj,nnn,nsndlev,njcsnd,kob
1674	integer :: komin,komax,nn,nhi,nlo,nnjc
1675	integer :: i_s,i_e
1676	integer :: istq
1677	real :: gfactor,rfactor,gridx,gridy,rindx,ris
1678	real :: grfacx,grfacy
1679	real :: timewt,pob
1680	real :: ri,rj,rx,ry,rsq,pdfac,erfivr,dk,slope,rinfac
1681	real :: dprim,dsq,d ! vars for mm5 surface-ob weighting
1682	real :: rinprs,pijk,pobhi,poblo,pdiffj,w2eowt,gitq
1683	real :: dz_ramp ! For ramping weights for surface obs
1684
1685	real :: scratch
1686	integer :: kk !ajb temp
1687
1688	! print *,'start nudob, nstat,j,ivar=',nstat,j,ivar
1689	! if(ivar.ne.4)return
1690	!yliu start -- for multi-scans: NSCAN=0: original
1691	! NSCAN=1: added a scan with a larger Ri and smaller G
1692	! if(NSCAN.ne.0 .and. NSCAN.ne.1) stop
1693	! ajb note: Will need to increase memory for SAVWT if NSCAN=1:
1694	if(NSCAN.ne.0) then
1695	IF (iprt) then
1696	write(msg,*) 'SAVWT must be resized for NSCAN=1'
1697	call wrf_message(msg)
1698	ENDIF
1699	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
1700	endif
1701	IPLO=0 + NSCAN*4
1702	GFACTOR=1. + NSCAN*(-1. + 0.33333)
1703	RFACTOR=1. + NSCAN*(-1. + 3.0)
1704	!yliu end
1705	! jc
1706
1707	! return if too close to j boundary
1708	if(inest.eq.1.and.ivar.lt.3.and.(j.le.2.or.j.ge.jde-1)) then
1709	! write(6,*) '1 RETURN: IVAR = ',ivar,' J = ',j,
1710	! $ ' too close to boundary.'
1711	return
1712	endif
1713	if(inest.eq.1.and.ivar.ge.3.and.(j.le.2.or.j.ge.jde-2)) then
1714	! write(6,*) '2 RETURN: IVAR = ',ivar,' J = ',j,
1715	! $ ' too close to boundary.'
1716	return
1717	endif
1718
1719	! COMPUTE IPL WHICH REPRESENTS IVAR FOR EACH MESH IN SAVWT MODS
1720	ICUT=0
1721	IF(INEST.GT.1)ICUT=1
1722	i_s = max0(2+icut,its)
1723	i_e = min0(ide-1-icut,ite)
1724
1725	IPL=IVAR + IPLO !yliu +IPLO
1726
1727	! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
1728
1729	INPF=(IRATIO*LEVIDN(INEST))NPFI
1730	INFR=(IRATIO*LEVIDN(INEST))IONF
1731
1732	GRIDX=0.0
1733	GRIDY=0.0
1734	IGRID=0
1735	IF(IVAR.GE.3)THEN
1736	GRIDX=0.5
1737	GRIDY=0.5
1738	IGRID=1
1739	ELSEIF(IVAR.eq.1) THEN
1740	GRIDY=0.5
1741	GRIDX=0.0
1742	IGRID=1
1743	ELSEIF(IVAR.eq.2) THEN
1744	GRIDX=0.5
1745	GRIDY=0.0
1746	IGRID=1
1747	ENDIF
1748
1749	! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
1750	! KILOMETERS TO GRID LENGTHS, RINDX
1751
1752	RINDX=RINXY1000./DX RFACTOR !yliu *RFACTOR
1753	RIS=RINDX*RINDX
1754	IF(IFREST.AND.KTAU.EQ.KTAUR)GOTO 5
1755	IF(MOD(KTAU,INFR).NE.0)GOTO 126
1756	5 CONTINUE
1757	IF (iprt) THEN
1758	IF(J.EQ.10) then
1759	write(msg,6) INEST,J,KTAU,XTIME,IVAR,IPL,rindx
1760	call wrf_message(msg)
1761	ENDIF
1762	ENDIF
1763	6 FORMAT(1X,'OBS NUDGING FOR IN,J,KTAU,XTIME,', &
1764	'IVAR,IPL: ',I2,1X,I2,1X,I5,1X,F8.2,1X,I2,1X,I2, &
1765	' rindx=',f4.1)
1766
1767	!********************************************************************
1768	!ajb_07012008 Setting ra and rb for the fine mesh her is now simple.
1769	! Values are no longer calculated here based on the
1770	! coarse mesh, since direct use of WRF map projections
1771	! on each nest was implemented in subroutine in4dob.
1772	!********************************************************************
1773	! SET RA AND RB
1774	DO N=1,NSTAT
1775	RA(N)=RIO(N)-GRIDX
1776	RB(N)=RJO(N)-GRIDY
1777	ENDDO
1778
1779	! INITIALIZE WEIGHTING ARRAYS TO ZERO
1780	DO I=its,ite
1781	DO K=1,kte
1782	WT(I,K)=0.0
1783	WT2ERR(I,K)=0.0
1784	ENDDO
1785	ENDDO
1786
1787	! DO P* COMPUTATIONS ON DOT POINTS FOR IVAR.LT.3 (U AND V)
1788	! AND CROSS POINTS FOR IVAR.GE.3 (T,Q,P*).
1789	!
1790	! COMPUTE P* AT OBS LOCATION (RA,RB). DO THIS AS SEPARATE VECTOR LOOP H
1791	! SO IT IS ALREADY AVAILABLE IN NSTAT LOOP 120 BELOW
1792
1793	! PSURF IS NOT AVAILABLE GLOBALLY, THEREFORE, THE BILINEAR INTERPOLATION
1794	! AROUND THE OBS POINT IS DONE IN ERROB() AND STORED IN ERRF([678],N) FOR
1795	! THE POINT (6=PRESS, 7=U-MOM, 8=V-MOM).
1796	! ajb 05052009: psurf is actually pbase(k=1) interpolated to obs (i,j).
1797	DO N=1,NSTAT
1798	IF(IVAR.GE.3)THEN
1799	PSURF(N)=ERRF(6,N)
1800	ELSE
1801	IF(IVAR.EQ.1)THEN
1802	PSURF(N)=ERRF(7,N) ! U-points
1803	ELSE
1804	PSURF(N)=ERRF(8,N) ! V-points
1805	ENDIF
1806	ENDIF
1807	ENDDO
1808
1809	! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
1810	! J-STRIP
1811
1812	MAXJ=J+IFIX(RINDX*RINFMX+0.99) !ajb
1813	MINJ=J-IFIX(RINDX*RINFMX+0.99) !ajb
1814
1815	! jc comment out this? want to use obs beyond the domain?
1816	! MAXJ=MIN0(JL-IGRID,MAXJ) !yliu
1817	! MINJ=MAX0(1,MINJ) !yliu
1818
1819	n=1
1820
1821	!***********************************************************************
1822	DO nnn=1,NSTAT ! BEGIN OUTER LOOP FOR THE NSTAT OBSERVATIONS
1823	!***********************************************************************
1824	! Soundings are consecutive obs, but they need to be treated as a single
1825	! entity. Thus change the looping to nnn, with n defined separately.
1826
1827
1828	!yliu
1829	! note for sfc data: nsndlev=1 and njcsnd=1
1830	nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
1831
1832	! yliu start -- set together with the other parts
1833	! test: do the sounding levels as individual obs
1834	! nsndlev=1
1835	! yliu end
1836	njcsnd=nsndlev
1837	! set pob here, to be used later
1838	pob=varobs(5,n)
1839
1840	! print *, "s-- n=,nsndlev",n,njcsnd,J, ipl
1841	! print *, "s--",ivar,(errf(ivar,i),i=n,n+njcsnd)
1842	! CHECK TO SEE OF STATION N HAS DATA FOR VARIABLE IVAR
1843	! AND IF IT IS SUFFICIENTLY CLOSE TO THE J STRIP. THIS
1844	! SHOULD ELIMINATE MOST STATIONS FROM FURTHER CONSIDER-
1845	! ATION.
1846
1847	!yliu: Skip bad obs if it is sfc or single level sounding.
1848	!yliu: Before this (020918), a snd will be skipped if its first
1849	!yliu level has bad data- often true due to elevation
1850
1851	IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
1852	! print *, " bad obs skipped"
1853
1854	ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
1855	! print *, " skipped obs far away from this J-slice"
1856
1857	!----------------------------------------------------------------------
1858	ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
1859	!----------------------------------------------------------------------
1860
1861	! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
1862	! FOR THE VERTICAL WEIGHTING, WTSIG
1863
1864	! ASSIMILATE OBSERVATIONS ON PRESSURE LEVELS, EXCEPT FOR SURFACE
1865	!ajb 20021210: (Bugfix) RKO is not available globally. It is computed in
1866	!ajb ERROB() by the processor handling the obs point, and stored in ERRF(9,N).
1867
1868	#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
1869	rko(n) = errf(9,n) !ajb 20021210
1870	kpbl_obs(n) = errf(5,n) !ajb 20090519
1871	#endif
1872	!ajb 20090427 added .45 to rko so KOB is equal to 1 only if RKO > 1.05
1873	KOB=nint(RKO(N)+0.45)
1874	KOB=MIN0(kte,KOB)
1875	KOB=MAX0(1,KOB)
1876
1877	! ASSIMILATE SURFACE LAYER DATA ON SIGMA
1878	IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
1879
1880	! Compute temporal weight
1881	timewt = get_timewt(xtime,dtmin,twindo,sfcfact,timeob(n))
1882
1883	DO K=1,kte
1884	WTSIG(K)=0.0
1885	ENDDO
1886	! DEFINE WTSIG: (FOR SRP: SPREAD SURFACE DATA THROUGH LOWEST 200 M)
1887	! WTSIG(1)=1.0
1888	! WTSIG(2)=0.67
1889	! WTSIG(3)=0.33
1890	! KOMIN=3
1891	! KOMAX=1
1892	! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
1893	! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX (IN GRID LENGTHS).
1894	! fix this because kpblt at 1 and il is 0
1895	MAXI=IFIX(RA(N)+0.99+RINDX*sfcfacr)
1896	MAXI=MIN0(ide-1,MAXI)
1897	MINI=IFIX(RA(N)-RINDX*sfcfacr-0.99)
1898	MINI=MAX0(2,MINI)
1899	!yliu start
1900	! use also obs outside of this domain -- surface obs
1901	! if(RA(N).LT.0.-RINDX .or. RA(N).GT.float(IL+RINDX) .or.
1902	! & RB(N).LT.0.-RINDX .or. RB(N).GT.float(JL+RINDX)) then
1903	! print *, " skipped obs far away from this domain"
1904	! currently can use obs within this domain or ones very close to (1/3
1905	! influence of radius in the coarse domain) this
1906	! domain. In later case, use BC column value to approximate the model value
1907	! at obs point -- ERRF need model field in errob.F !!
1908	if ( RA(N).GE.(0.-RINDX*sfcfacr/3) &
1909	.and. RA(N).LE.float(ide)+RINDX*sfcfacr/3 &
1910	.and. RB(N).GE.(0.-RINDX*sfcfacr/3) &
1911	.and. RB(N).LE.float(jde)+RINDX*sfcfacr/3) then
1912
1913	! or use obs within this domain only
1914	! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
1915	! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
1916	! print *, " skipped obs far outside of this domain"
1917	! if(j.eq.3 .and. ivar.eq.3) then
1918	! write(6,*) 'N = ',n,' exit 120 3'
1919	! endif
1920	!yliu end
1921	!
1922	! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
1923	! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
1924	! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
1925	RJ=FLOAT(J)
1926	RX=RJ-RB(N)
1927	! WEIGHTS FOR THE 3-D VARIABLES
1928	ERFIVR=ERRF(IVAR,N)
1929
1930	!ajb Compute and add weights to sum only if nudge_pbl switch is on.
1931	if(nudge_pbl) then
1932
1933	! Apply selected horizontal spreading scheme.
1934	if(SFC_SCHEME_HORIZ.eq.1) then ! original mm5 scheme
1935	DO I=max0(its,MINI),min0(ite,MAXI)
1936	RI=FLOAT(I)
1937	RY=RI-RA(N)
1938	RIS=RINDXRINDXsfcfacr*sfcfacr
1939	RSQ=RXRX+RYRY
1940	! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
1941	DPRIM=SQRT(RSQ)
1942	D=DPRIM+RINDXDCONABS(psurf(n)-pbase(i,1))
1943	DSQ=D*D
1944	WTIJ(i)=(RIS-DSQ)/(RIS+DSQ)
1945	WTIJ(i)=AMAX1(0.0,WTIJ(i))
1946	ENDDO
1947	else ! wrf scheme
1948	DO I=max0(its,MINI),min0(ite,MAXI)
1949	RI=FLOAT(I)
1950	RY=RI-RA(N)
1951	RIS=RINDXRINDXsfcfacr*sfcfacr
1952	RSQ=RXRX+RYRY
1953	! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
1954	wtij(i)=(ris-rsq)/(ris+rsq)
1955	scratch = (abs (psurf(n)-.001pbase(i,1))DCON)
1956	pdfac=1.-AMIN1(1.0,scratch)
1957	wtij(i)=wtij(i)*pdfac
1958	WTIJ(i)=AMAX1(0.0,WTIJ(i))
1959	ENDDO
1960	endif
1961
1962	! Apply selected vertical spreading scheme.
1963	if(SFC_SCHEME_VERT.eq.1) then ! original simple scheme
1964
1965	DO I=max0(its,MINI),min0(ite,MAXI)
1966
1967	! try making sfc obs weighting go thru pbl
1968	komax=max0(3,kpblt(i))
1969	IF (iprt) THEN
1970	if (kpblt(i).gt.25 .and. ktau.ne.0) &
1971	write(6,552)inest,i,j,kpblt(i)
1972	552 FORMAT('kpblt is gt 25, inest,i,j,kpblt=',4i4)
1973	ENDIF
1974
1975	if(kpblt(i).gt.25) komax=3
1976	komin=1
1977	dk=float(komax)
1978
1979	do k=komin,komax
1980
1981	wtsig(k)=float(komax-k+1)/dk
1982	WT(I,K)=WT(I,K)+TIMEWTWTSIG(K)WTIJ(i)
1983
1984	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJ(i)WTIJ(i)WTSIG(K) &
1985	WTSIG(K)ERFIVR
1986	enddo
1987	ENDDO
1988
1989	else ! regime-based vif scheme
1990
1991	! Here we calculate weights in the vertical coordinate, based on vih1 and vih2.
1992	! In the equation for wtsig(k), Z=zslab(i,k)-terrh(i) contains the model Z-values
1993	! (height above ground in meters) on a J-slab. The equation produces wtsig = 1.0 at
1994	! levels 1 to K, where z(K) < vih1 < z(K+1). For the example below, in the equation
1995	! for wtsig(k), the expression vih1(i)-Z(i,k) is strictly positive for k=1,2,3 since
1996	! levels 1, 2, and 3 are below vih1. So xtsig(k)=min(1.0, 1.0-x) where x > 0 ==>
1997	! wtsig(k)=1 for k=1,2,3.
1998	!
1999	! For levels K+1 and up, wtsig will decrease linearly with height, with values
2000	! along the ramp that has value 1.0 at vih1 and 0.0 at vih2. In the example:
2001	!
2002	! dz_ramp = 1/(200-150) = 1/50
2003	! xtsig(4) = 1 + (150-175)/50 = 1 - 1/2 = 1/2
2004	!
2005	! WTSIG
2006	! 1 -\|* * * * * *
2007	! \|
2008	! \| *
2009	! \|
2010	! \| *
2011	! \|
2012	! \| *
2013	! 0 -\|--\|-------\|-----------\|------\|----\|----\|---------\|----> Z = HT ABOVE
2014	! 15 55 115 150 175 200 250 GROUND
2015	! k=1 k=2 k=3 vih1 k=4 vih2 k=5
2016
2017	DO I=max0(its,MINI),min0(ite,MAXI)
2018
2019	dz_ramp = 1.0 / max( 1.0, vih2(i)-vih1(i) ) ! vih2 >= vih1 by construct
2020
2021	LML: do k = kts, kte
2022	wtsig(k) = min( 1.0, 1.0 + ( vih1(i)-zslab(i,k)+terrh(i) ) * dz_ramp )
2023	wtsig(k) = max( 0.0, wtsig(k))
2024
2025	if(wtsig(k).le.0.0) EXIT LML
2026	WT(I,K)=WT(I,K)+TIMEWTWTSIG(K)WTIJ(i)
2027	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJ(i)WTIJ(i)WTSIG(K) &
2028	WTSIG(K)ERFIVR
2029	enddo LML
2030	ENDDO
2031	endif
2032
2033	endif ! end if nudge-pbl switch is on
2034
2035	endif ! end check for obs in domain
2036	! END SURFACE-LAYER OBS NUDGING
2037
2038	ELSE
2039
2040	! Compute temporal weight
2041	timewt = get_timewt(xtime,dtmin,twindo,1.,timeob(n))
2042
2043
2044	! BEGIN CALCULATIONS TO SPREAD OBS INFLUENCE ALONG PRESSURE LEVELS
2045	!
2046	! print *,'in upper air section'
2047	! DEFINE THE MAX AND MIN I VALUES FOR POSSIBLE NONZERO
2048	! WEIGHTS, BASED ON THE RADIUS OF INFLUENCE, RINDX, AND RINFAC.
2049	! RINFAC VARIES AS A LINEAR FUNCTION FROM FROM RINFMN AT P*+PTOP
2050	! TO RINFMX AT PFREE AND "ABOVE" (LOWER PRESSURE).
2051	!ajb SLOPE=(RINFMN-RINFMX)/(PSBO(N)+PTOP-PFREE)
2052
2053	slope = (RINFMN-RINFMX)/(psurf(n)-PFREE)
2054
2055	RINFAC=SLOPEPOB+RINFMX-SLOPEpfree
2056	RINFAC=AMAX1(RINFAC,RINFMN)
2057	RINFAC=AMIN1(RINFAC,RINFMX)
2058	!yliu: for multilevel upper-air data, take the maximum
2059	! for the I loop.
2060	if(nsndlev.gt.1) RINFAC = RINFMX
2061	!yliu end
2062
2063	MAXI=IFIX(RA(N)+0.99+RINDX*RINFAC)
2064	MAXI=MIN0(ide-IGRID,MAXI)
2065	MINI=IFIX(RA(N)-RINDX*RINFAC-0.99)
2066	MINI=MAX0(1,MINI)
2067
2068	! yliu start
2069	! use also obs outside of but close to this domain -- upr data
2070	! if( RA(N).LT.(0.-RINFAC*RINDX)
2071	! & .or. RA(N).GT.float(IL)+RINFAC*RINDX
2072	! & .or. RB(N).LT.(0.-RINFAC*RINDX)
2073	! & .or. RB(N).GT.float(JL)+RINFAC*RINDX)then
2074	! print *, " skipped obs far away from this I-range"
2075	! currently can use obs within this domain or ones very close to (1/3
2076	! influence of radius in the coarse domain) this
2077	! domain. In later case, use BC column value to approximate the model value
2078	! at obs point -- ERRF need model field in errob.F !!
2079
2080	!cc if (i.eq.39 .and. j.eq.34) then
2081	!cc write(6,*) 'RA(N) = ',ra(n)
2082	!cc write(6,*) 'rinfac = ',rinfac,' rindx = ',rindx
2083	!cc endif
2084	if( RA(N).GE.(0.-RINFAC*RINDX/3) &
2085	.and. RA(N).LE.float(ide)+RINFAC*RINDX/3 &
2086	.and. RB(N).GE.(0.-RINFAC*RINDX/3) &
2087	.and. RB(N).LE.float(jde)+RINFAC*RINDX/3) then
2088	! or use obs within this domain only
2089	! if(RA(N).LT.1 .or. RA(N).GT.float(IL) .or.
2090	! & RB(N).LT.1 .or. RB(N).GT.float(JL)) then
2091	! print *, " skipped obs far outside of this domain"
2092
2093	! yliu end
2094	! is this 2 needed here - kpbl not used?
2095	! MINI=MAX0(2,MINI)
2096
2097	! LOOP THROUGH THE NECESSARY GRID POINTS SURROUNDING
2098	! OBSERVATION N. COMPUTE THE HORIZONTAL DISTANCE TO
2099	! THE OBS AND FIND THE WEIGHTING SUM OVER ALL OBS
2100	RJ=FLOAT(J)
2101	RX=RJ-RB(N)
2102	! WEIGHTS FOR THE 3-D VARIABLES
2103	!
2104	ERFIVR=ERRF(IVAR,N)
2105	! jc
2106	nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
2107	! yliu start
2108	! test: do the sounding levels as individual obs
2109	! nsndlev=1
2110	! yliu end
2111	njcsnd=nsndlev
2112	!
2113	DO I=max0(its,MINI),min0(ite,MAXI)
2114	! jc
2115	RI=FLOAT(I)
2116	RY=RI-RA(N)
2117	RIS=RINDXRINFACRINDX*RINFAC
2118	RSQ=RXRX+RYRY
2119	! yliu test: for upper-air data, keep D1 influence radii
2120	WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
2121
2122	WTIJ=AMAX1(0.0,WTIJ(i))
2123	! weight ob in vertical with +- 50 mb
2124	! yliu: 75 hba for single upper-air, 30hba for multi-level soundings
2125	if(nsndlev.eq.1) then
2126	rinprs=7.5
2127	!
2128	else
2129	rinprs=3.0
2130	endif
2131	! yliu end
2132
2133	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2134	! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
2135	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2136
2137	if(nsndlev.eq.1)then
2138	!----------------------------------------------------------------------
2139	! --- HANDLE 1-LEVEL OBSERVATIONS ---
2140	!----------------------------------------------------------------------
2141
2142	! if(I.eq.MINI) print *, " Single snd "
2143	! ERFIVR is the residual (difference) between the ob and the model
2144	! at that point. We can analyze that residual up and down.
2145	! First find komin for ob.
2146	!yliu start -- in the old code, komax and komin were reversed!
2147	do k=kte,1,-1
2148	pijk = .001*(pbase(i,k)+pp(i,k))
2149	! print *,'k,pijk,pob,rinprs=',k,pijk,pob,rinprs
2150	if(pijk.ge.(pob+rinprs)) then
2151	komin=k
2152	go to 325
2153	endif
2154	enddo
2155	komin=1
2156	325 continue
2157	! now find komax for ob
2158	do k=3,kte
2159	pijk = .001*(pbase(i,k)+pp(i,k))
2160	if(pijk.le.(pob-rinprs)) then
2161	komax=k
2162	go to 326
2163	endif
2164	enddo
2165	komax=kte ! yliu 20050706
2166	326 continue
2167
2168	! yliu: single-level upper-air data will act either above or below the PBL top
2169	! ajb: Reset komin or komax. if kobs>kpbl_obs, komin=kpbl_obs+1, else komax=kpbl_obs
2170
2171	if( (kpblt(i).le.komax) .and. (kpblt(i).ge.komin) ) then
2172	kobs = 1
2173	OBS_K: do k = komin, komax
2174	if( pob .gt. .001*(pbase(i,k)+pp(i,k)) ) then
2175	kobs = k
2176	EXIT OBS_K
2177	endif
2178	enddo OBS_K
2179
2180	if(kobs.gt.kpbl_obs(n)) then
2181	! Obs will act only above the PBL top
2182	komin=max0(kobs, komin) ! kobs here is kpblt(i)+1
2183	else ! Obs acts below PBL top
2184	! Obs will act only below the PBL top
2185	komax=min0(kpblt(i), komax)
2186	endif
2187	endif
2188	! yliu end
2189	!
2190	! print *,'1 level, komin,komax=',komin,komax
2191	! if(i.eq.MINI) then
2192	! print *, "yyyyyyyyyyS IPL erfivr=", IPL, ERFIVR,J,pob
2193	! ERFIVR=0
2194	! endif
2195	do k=1,kte
2196	reserf(k)=0.0
2197	wtsig(k)=0.0
2198	enddo
2199	!yliu end
2200
2201	!ajb Add weights to sum only if nudge_pbl switch is on OR obs is above pbl top.
2202	if(nudge_pbl .or. komin.ge.kpblt(i)) then
2203	do k=komin,komax
2204	pijk = .001*(pbase(i,k)+pp(i,k))
2205	reserf(k)=erfivr
2206	wtsig(k)=1.-abs(pijk-pob)/rinprs
2207	wtsig(k)=amax1(wtsig(k),0.0)
2208	! print *,'k,pijk,pob,rinprs,wtsig=',k,pijk,pob,rinprs,wtsig(k)
2209	! Now calculate WT and WT2ERR for each i,j,k point cajb
2210	WT(I,K)=WT(I,K)+TIMEWTWTIJ(i)wtsig(k)
2211
2212	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJ(i)WTIJ(i) &
2213	reserf(k)wtsig(k)wtsig(k)
2214	enddo
2215	endif
2216
2217	else
2218	!----------------------------------------------------------------------
2219	! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
2220	!----------------------------------------------------------------------
2221	!yliu start
2222	! if(I.eq.MINI) print *, " Multi-level snd "
2223	! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
2224	! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
2225	if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
2226	IF (iprt) THEN
2227	write(msg,*) "n = ",n,"nsndlev = ",nsndlev
2228	call wrf_message(msg)
2229	write(msg,*) "nlevs_ob,lev_in_ob", &
2230	nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
2231	call wrf_message(msg)
2232	call wrf_message("in nudobs.F: sounding level messed up, stopping")
2233	ENDIF
2234	call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
2235	endif
2236	!yliu end
2237	! This is for a multi-level observation
2238	! The trick here is that the sounding is "one ob". You don't
2239	! want multiple levels to each be treated like separate
2240	! and independent observations.
2241	! At each i,j want to interpolate sounding to the model levels at that
2242	! particular point.
2243	komin=1
2244	komax=kte-2
2245	! this loop goes to 1501
2246	! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
2247	!yliu start
2248	do k=1,kte
2249	reserf(k)=0.0
2250	wtsig(k)=0.0
2251	enddo
2252	!yliu end
2253
2254	do k=komax,komin,-1
2255
2256	pijk = .001*(pbase(i,k)+pp(i,k))
2257
2258	! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
2259	if(pijk.gt.varobs(5,n)) then
2260	go to 1501
2261	endif
2262
2263	! if sigma level pressure is .lt. than the highest ob level, don't interpolate
2264	if(pijk.le.varobs(5,n+nsndlev-1)) then
2265	go to 1501
2266	endif
2267
2268	! now interpolate sounding to this k
2269	! yliu start-- recalculate WTij for each k-level
2270	!ajb SLOPE = (RINFMN-RINFMX)/(pdoc(i,j)+PTOP-PFREE)
2271	slope = (RINFMN-RINFMX)/ (.001*pbase(i,1)-PFREE)
2272	RINFAC=SLOPEpijk+RINFMX-SLOPEPFREE
2273	RINFAC=AMAX1(RINFAC,RINFMN)
2274	RINFAC=AMIN1(RINFAC,RINFMX)
2275	RIS=RINDXRINFACRINDX*RINFAC
2276	RSQ=RXRX+RYRY
2277
2278	! for upper-air data, keep D1 influence radii
2279	WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
2280	WTIJ(i)=AMAX1(0.0,WTIJ(i))
2281	! yliu end
2282
2283	! this loop goes to 1503
2284	do nn=2,nsndlev
2285
2286	! only set pobhi if varobs(ivar) is ok
2287	pobhi=-888888.
2288
2289	if(varobs(ivar,n+nn-1).gt.-800000. &
2290	.and. varobs(5,n+nn-1).gt.-800000.) then
2291	pobhi=varobs(5,n+nn-1)
2292	nhi=n+nn-1
2293	if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.20.) then
2294	go to 1502 ! within 200mb of obs height
2295	endif
2296	endif
2297
2298	enddo
2299
2300	! did not find any ob above within 200 mb, so jump out
2301	go to 1501
2302	1502 continue
2303
2304	nlo=nhi-1
2305	do nnjc=nhi-1,n,-1
2306	if(varobs(ivar,nnjc).gt.-800000. &
2307	.and. varobs(5,nnjc).gt.-800000.) then
2308	poblo=varobs(5,nnjc)
2309	nlo=nnjc
2310	if(poblo.gt.pijk .and. abs(poblo-pijk).lt.20.) then
2311	go to 1505 ! within 200mb of obs height
2312	endif
2313	endif
2314	enddo
2315	!yliu end --
2316
2317	! did not find any ob below within 200 mb, so jump out
2318	go to 1501
2319	1505 continue
2320
2321	! interpolate to model level
2322	pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
2323	reserf(k)=errf(ivar,nlo)+ &
2324	(errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
2325	wtsig(k)=1.
2326
2327	1501 continue
2328
2329	! now calculate WT and WT2ERR for each i,j,k point cajb
2330	!ajb Add weights to sum only if nudge_pbl switch is on OR k > kpblt.
2331	if(nudge_pbl .or. k.gt.kpblt(i)) then
2332
2333	WT(I,K)=WT(I,K)+TIMEWTWTIJ(i)wtsig(k)
2334
2335	WT2ERR(I,K)=WT2ERR(I,K)+TIMEWTTIMEWTWTIJ(i)WTIJ(i) &
2336	reserf(k)wtsig(k)wtsig(k)
2337	endif
2338
2339	enddo ! enddo k levels
2340
2341	! end multi-levels
2342	endif ! end if(nsndlev.eq.1)
2343	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2344	! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
2345	!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2346	!
2347	ENDDO ! END DO MINI,MAXI LOOP
2348
2349	endif ! check for obs in domain
2350
2351	! END OF NUDGING TO OBS ON PRESSURE LEVELS
2352
2353	ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
2354
2355	!----------------------------------------------------------------------
2356	ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
2357	!----------------------------------------------------------------------
2358
2359	! n=n+1
2360	n=n+njcsnd
2361
2362	!yliu 1202 continue
2363	if(n.gt.nstat)then
2364	! print *,'n,nstat=',n,nstat,ivar,j
2365	go to 1203
2366	endif
2367	! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
2368
2369	!***********************************************************************
2370	ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
2371	!***********************************************************************
2372
2373	1203 continue
2374
2375	! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
2376	! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
2377	! THE ATEN ARRAY
2378	! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
2379	! THEY ARE USED BELOW IN THE DENOMINATOR.
2380	DO K=kts,kte
2381	DO I=its,ite
2382	IF(WT(I,K).EQ.0)THEN
2383	WT2ERR(I,K)=0.0
2384	ENDIF
2385	IF(WT(I,K).EQ.0)THEN
2386	WT(I,K)=1.0
2387	ENDIF
2388	ENDDO
2389	ENDDO
2390
2391	126 CONTINUE
2392
2393	IF(IVAR.GE.3)GOTO 170
2394	! this is for u,v
2395	! 3-D DOT POINT TENDENCIES
2396
2397	! Calculate scales for converting nudge factor from u (v)
2398	! to rho_u (or rho_v) units.
2399
2400	IF (IVAR == 1) THEN
2401	call calc_rcouple_scales(mu,msfy,rscale,ims,ime,its,ite)
2402	ELSE IF (IVAR == 2) THEN
2403	call calc_rcouple_scales(mu,msfx,rscale,ims,ime,its,ite)
2404	END IF
2405
2406	DO K=1,kte
2407
2408	DO I=i_s,i_e
2409
2410	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2411	W2EOWT=WT2ERR(I,K)/WT(I,K)
2412	ELSE
2413	W2EOWT=SAVWT(IPL,I,K)
2414
2415	! write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,W2EOWT = ',i,j,k,ipl,W2EOWT
2416
2417	ENDIF
2418
2419	! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2420	! scratch = GIVRSCALE(I)W2EOWTTFACIISWIND*GFACTOR
2421	! write(6,*) 'ATEN calc: k = ',k
2422	! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
2423	! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2424	! $ ' W2EOWT = ',w2eowt
2425	! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
2426	! $ ' GFACTOR = ',gfactor
2427	! endif
2428	!
2429	! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2430	! scratch = GIVRSCALE(I)W2EOWTTFACIISWIND*GFACTOR
2431	! write(6,*) 'ATEN calc: k = ',k
2432	! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
2433	! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2434	! $ ' W2EOWT = ',w2eowt
2435	! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
2436	! $ ' GFACTOR = ',gfactor
2437	! endif
2438
2439	! if(ivar .eq. 1 .and. i.eq.38 .and. (j.ge.36.and.j.le.62) .and. k.eq.7) then
2440	! scratch = GIVRSCALE(I)W2EOWTTFACIISWIND*GFACTOR
2441	! write(6,*)
2442	! write(6,*) 'ATEN calc: j = ',j,' k = ',k
2443	! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
2444	! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),' W2EOWT = ',w2eowt
2445	! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,' GFACTOR = ',gfactor
2446	! endif
2447
2448	ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
2449	W2EOWTTFACI &
2450	ISWIND GFACTOR !yliu *GFACTOR
2451
2452	! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2453	! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
2454	! endif
2455	! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2456	! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
2457	! endif
2458
2459	ENDDO
2460	ENDDO
2461
2462	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2463	DO K=1,kte
2464	DO I=its,ite
2465	SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2466
2467	! write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,savwt = ',i,j,k,ipl,savwt(ipl,i,k)
2468	ENDDO
2469	ENDDO
2470	ENDIF
2471
2472	RETURN
2473
2474	170 CONTINUE
2475
2476	! 3-D CROSS-POINT TENDENCIES
2477	! this is for t (ivar=3) and q (ivsr=4)
2478	IF(3-IVAR.LT.0)THEN
2479	GITQ=GIQ
2480	ELSE
2481	GITQ=GIT
2482	ENDIF
2483	IF(3-IVAR.LT.0)THEN
2484	ISTQ=ISMOIS
2485	ELSE
2486	ISTQ=ISTEMP
2487	ENDIF
2488
2489	DO K=1,kte
2490	DO I=i_s,i_e
2491	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2492	W2EOWT=WT2ERR(I,K)/WT(I,K)
2493	ELSE
2494	W2EOWT=SAVWT(IPL,I,K)
2495	ENDIF
2496
2497	! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
2498	! scratch = GITQMU(I)W2EOWTTFACIISTQ*GFACTOR
2499	! write(6,*) 'ATEN calc: k = ',k
2500	! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
2501	! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i), &
2502	! ' W2EOWT = ',w2eowt
2503	! write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq, &
2504	! ' GFACTOR = ',gfactor
2505	! endif
2506
2507	! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2508	! scratch = GITQMU(I)W2EOWTTFACIISTQ*GFACTOR
2509	! write(6,*) 'ATEN calc: k = ',k
2510	! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
2511	! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
2512	! $ ' W2EOWT = ',w2eowt
2513	! write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq,
2514	! $ ' GFACTOR = ',gfactor
2515	! endif
2516
2517	ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
2518	W2EOWTTFACIISTQ GFACTOR !yliu *GFACTOR
2519
2520	! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
2521	! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
2522	! endif
2523	! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2524	! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
2525	! endif
2526
2527	ENDDO
2528	ENDDO
2529
2530	IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
2531	DO K=1,kte
2532	DO I=its,ite
2533	SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2534	ENDDO
2535	ENDDO
2536	ENDIF
2537
2538	RETURN
2539	END SUBROUTINE nudob
2540
2541	SUBROUTINE date_string(year, month, day, hour, minute, second, cdate)
2542	!-----------------------------------------------------------------------
2543	! PURPOSE: Form a date string (YYYY-MM-DD_hh:mm:ss) from integer
2544	! components.
2545	!-----------------------------------------------------------------------
2546	IMPLICIT NONE
2547	!-----------------------------------------------------------------------
2548
2549	INTEGER, INTENT(IN) :: year
2550	INTEGER, INTENT(IN) :: month
2551	INTEGER, INTENT(IN) :: day
2552	INTEGER, INTENT(IN) :: hour
2553	INTEGER, INTENT(IN) :: minute
2554	INTEGER, INTENT(IN) :: second
2555	CHARACTER*19, INTENT(INOUT) :: cdate
2556
2557	! Local variables
2558	integer :: ic ! loop counter
2559
2560	cdate(1:19) = "0000-00-00_00:00:00"
2561	write(cdate( 1: 4),'(i4)') year
2562	write(cdate( 6: 7),'(i2)') month
2563	write(cdate( 9:10),'(i2)') day
2564	write(cdate(12:13),'(i2)') hour
2565	write(cdate(15:16),'(i2)') minute
2566	write(cdate(18:19),'(i2)') second
2567	do ic = 1,19
2568	if(cdate(ic:ic) .eq. " ") cdate(ic:ic) = "0"
2569	enddo
2570
2571	RETURN
2572	END SUBROUTINE date_string
2573
2574	SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
2575	!-----------------------------------------------------------------------
2576	IMPLICIT NONE
2577	!-----------------------------------------------------------------------
2578
2579	INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
2580	INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
2581	REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
2582	REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
2583	REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
2584
2585	! Local variables
2586	integer :: i
2587
2588	! Calculate scales to be used for producing rho-coupled nudging factors.
2589	do i = its,ite
2590	rscale(i) = a(i)/msf(i)
2591	enddo
2592
2593	RETURN
2594	END SUBROUTINE calc_rcouple_scales
2595
2596	SUBROUTINE print_obs_info(iprt,inest,niobf,rio,rjo,rko, &
2597	prt_max,prt_freq,obs,stnid,lat,lon, &
2598	mlat,mlon,timeob,xtime)
2599	!*************************************************************************
2600	! Purpose: Print obs information.
2601	!*************************************************************************
2602
2603	IMPLICIT NONE
2604
2605	LOGICAL, intent(in) :: iprt ! Print flag
2606	INTEGER, intent(in) :: inest ! Nest level
2607	INTEGER, intent(in) :: niobf ! Maximum number of observations
2608	REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
2609	REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
2610	REAL, intent(in) :: rko(niobf) ! Bottom-top north coord (non-stagger)
2611	INTEGER, intent(in) :: prt_max ! Max no. of obs for diagnostic printout
2612	INTEGER, intent(in) :: prt_freq ! Frequency for diagnostic printout
2613	INTEGER, intent(in) :: obs(prt_max) ! Saved obs indices to print
2614	INTEGER, intent(in) :: stnid(40,prt_max) ! Saved station ids
2615	REAL, intent(in) :: lat(prt_max) ! Saved latitudes
2616	REAL, intent(in) :: lon(prt_max) ! Saved longitudes
2617	REAL, intent(in) :: mlat(prt_max) ! Saved model latitudes
2618	REAL, intent(in) :: mlon(prt_max) ! Saved longitudes
2619	REAL, intent(in) :: timeob(niobf) ! Times of each observation (hours)
2620	REAL, intent(in) :: xtime ! Model time in minutes
2621
2622	! Local variables
2623	integer :: i ! Loop counter over obs station chars
2624	integer :: n ! Loop counter over obs
2625	integer :: pnx ! Obs index for printout
2626	character(len=200) :: msg ! Argument to wrf_message
2627	character(len=20) :: station_id ! Station id of observation
2628
2629	if(iprt) then
2630	if(prt_max.gt.0) then
2631
2632	if(obs(1).ne.-999) then
2633
2634	call wrf_message("")
2635	write(msg,fmt='(a,i4,a,f8.1,a)') 'REPORTING OBS MASS-PT LOCS FOR NEST ', &
2636	inest,' AT XTIME=',xtime,' MINUTES'
2637	call wrf_message(msg)
2638
2639	write(msg,fmt='(a,i4,a,i5,a)') 'FREQ=',prt_freq,', MAX=',prt_max, &
2640	' LOCS, NEWLY READ OBS ONLY, -999 => OBS OFF PROC'
2641	call wrf_message(msg)
2642	call wrf_message("")
2643
2644	write(msg,fmt='(3a)') ' OBS# I J K OBS LAT', &
2645	' OBS LON XLAT(I,J) XLONG(I,J) TIME(hrs)', &
2646	' OBS STATION ID'
2647	call wrf_message(msg)
2648
2649	endif
2650	endif
2651
2652	! Note: rio and rjo are referenced to non-staggered grid (not mass-point!)
2653	! Hence subtract .5 from each to get mass-point coords.
2654	do n=1,prt_max
2655	pnx = obs(n)
2656	if(pnx.ne.-999) then
2657	! Retrieve 15 chars of station id
2658	do i = 1,15
2659	station_id(i:i) = char(stnid(i,n))
2660	enddo
2661	write(msg,fmt='(2x,i7,3f8.3,2f9.3,2x,f9.3,2x,f9.3,3x,f6.2,7x,a15)') &
2662	pnx,rio(pnx)-.5,rjo(pnx)-.5,rko(pnx),lat(n),lon(n), &
2663	mlat(n),mlon(n),timeob(pnx),station_id
2664	call wrf_message(msg)
2665	endif
2666	enddo
2667	if(obs(1).ne.-999) call wrf_message("")
2668	endif
2669	END SUBROUTINE print_obs_info
2670
2671	REAL FUNCTION ht_to_p( h, pbbc, ppbc, z, ic, jc, dx, dy, &
2672	k_start, k_end, kds,kde, ims,ime, jms,jme, kms,kme )
2673
2674	!******************************************************************************
2675	! Purpose: Interpolate pressure at a specified x (ic), y (jc), and height (h).
2676	! The input pressure column pbbc+ppbc (base and perturbn) must already
2677	! be horizontally interpolated to the x, y position. The subroutine
2678	! get_height_column is called here to horizontally interpolated the
2679	! 3D height field z to get a height column at (iob, job).
2680	!******************************************************************************
2681
2682	IMPLICIT NONE
2683
2684	REAL, INTENT(IN) :: h ! height value (m)
2685	INTEGER, INTENT(IN) :: k_start, k_end ! loop bounds
2686	INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
2687	INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
2688	REAL, INTENT(IN) :: pbbc(kds:kde) ! column base pressure (cb)
2689	REAL, INTENT(IN) :: ppbc(kds:kde) ! column pressure perturbation (cb)
2690	REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! ht (m) above sl on half-levels
2691	INTEGER, INTENT(IN) :: ic ! i-coord of desired p
2692	INTEGER, INTENT(IN) :: jc ! j-coord of desired p
2693	REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
2694	REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
2695
2696	! Local variables
2697	INTEGER :: k ! loop counter
2698	INTEGER :: klo ! lower k bound
2699	REAL :: zlo ! lower z bound for h
2700	REAL :: zhi ! upper z bound for h
2701	REAL :: p ! interpolated pressure value
2702	REAL :: ln_p ! log p
2703	REAL :: ln_plo ! log plo
2704	REAL :: ln_phi ! log phi
2705	REAL :: z_at_p( kms:kme ) ! height at p levels
2706
2707	! Get interpolated z column on pressure (half-) levels at (ic,jc)
2708	call get_height_column(z, ic, jc, dx, dy, z_at_p, &
2709	k_start, k_end, kds,kde, &
2710	ims,ime, jms,jme, kms,kme )
2711
2712	! Now we have pbbc, ppbc, z_at_p, so compute p at h. First, find
2713	! bounding layers klo and khi so that z_at_p(klo) <= h <= z_at_p(khi)
2714
2715	ZLEVS: do k = k_start+1, k_end
2716	klo = k-1
2717	if(h .le. z_at_p(k)) then
2718	EXIT ZLEVS
2719	endif
2720	enddo ZLEVS
2721
2722	zlo = z_at_p(klo)
2723	zhi = z_at_p(klo+1)
2724
2725	! Interpolate natural log of pressure
2726	ln_plo = log( pbbc(klo+1) + ppbc(klo+1) )
2727	ln_phi = log( pbbc(klo) + ppbc(klo) )
2728	if(h.le.zlo) then
2729	ln_p = ln_phi ! set to k=1 pressure
2730	else if (h.ge.zhi) then
2731	ln_p = ln_plo ! set to k=k_end pressure
2732	else
2733	ln_p = ln_plo + (ln_phi-ln_plo)*((zhi-h)/(zhi-zlo))
2734	endif
2735
2736	! Return pressure
2737	p = exp(ln_p)
2738	ht_to_p = p
2739	RETURN
2740	END FUNCTION ht_to_p
2741
2742	SUBROUTINE get_height_column( z, ic, jc, dx, dy, z_at_p, &
2743	k_start, k_end, kds,kde, &
2744	ims,ime, jms,jme, kms,kme )
2745	!*************************************************************************
2746	! Purpose: Compute column height at ic, jc location on p levels
2747	!*************************************************************************
2748
2749	IMPLICIT NONE
2750
2751	INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
2752	INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
2753	INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
2754	REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! ht (m) on half-levels
2755	INTEGER, INTENT(IN) :: ic ! i-coord of desired p
2756	INTEGER, INTENT(IN) :: jc ! j-coord of desired p
2757	REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
2758	REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
2759	REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
2760
2761	! Local variables
2762	INTEGER :: k ! loop counter
2763
2764
2765	do k = kds, kde
2766	z_at_p(k) = &
2767	(1.-DY)( (1.-DX)z(IC,K,JC) + &
2768	DX *z(IC+1,K,JC) ) + &
2769	DY* ( (1.-DX)*z(IC,K,JC+1) + &
2770	DX *z(IC+1,K,JC+1) )
2771	enddo
2772
2773	END SUBROUTINE get_height_column
2774
2775	SUBROUTINE get_base_state_height_column( p_top, p00, t00, a, g, r_d, &
2776	znu, z_at_p, k_start, k_end, kds,kde, kms,kme )
2777	!********************************************************
2778	! Purpose: Compute base-state column height on p levels.
2779	! See, e.g., eqns 1.3-1.5 in MM5 User's Guide.
2780	! Height is a function of pressure plus the constants:
2781	! p00 - sea level pressure (Pa)
2782	! t00 - sea level temperature (K)
2783	! a - base state lapse rate (k/m)
2784	! r_d - gas constant (J/Kg/K) (Joule=1kg/m2/s2)
2785	! g - gravitational constant (m/s**2)
2786	!********************************************************
2787
2788	IMPLICIT NONE
2789
2790	REAL, INTENT(IN) :: p_top ! pressure at top of model
2791	REAL, INTENT(IN) :: p00 ! base state pressure
2792	REAL, INTENT(IN) :: t00 ! base state temperature
2793	REAL, INTENT(IN) :: a ! base state lapse rate
2794	REAL, INTENT(IN) :: g ! gravity constant
2795	REAL, INTENT(IN) :: r_d ! gas constant
2796	INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
2797	INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
2798	INTEGER, INTENT(IN) :: kms,kme ! vertical memory dim.
2799	REAL, INTENT(IN) :: znu( kms:kme ) ! eta values on half (mass) levels
2800	REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
2801
2802	! Local variables
2803	integer :: k ! loop counter
2804	real :: ps0 ! base state pressure at surface
2805	real :: pb(kms:kme) ! pressure on half eta levels
2806	real :: logterm ! temporary for Z calculation
2807	real :: ginv ! 1/g
2808
2809	ginv = 1/g
2810
2811	! Compute base state pressure on half eta levels.
2812	do k = k_start, k_end
2813	pb(k) = znu(k)*(p00 - p_top) + p_top
2814	enddo
2815
2816	! Use hydrostatic relation to compute height at pressure levels.
2817	do k = k_start, k_end
2818	logterm = log(pb(k)/p00)
2819	z_at_p(k) = .5r_daginvlogtermlogterm - r_dt00ginvlogterm
2820	enddo
2821
2822	END SUBROUTINE get_base_state_height_column
2823
2824	REAL FUNCTION get_timewt(xtime,dtmin,twindo,scalef,obtime)
2825	!*************************************************************************
2826	! Purpose: Compute the temporal weight factor for an observation
2827	!*************************************************************************
2828
2829	IMPLICIT NONE
2830
2831	REAL, INTENT(IN) :: xtime ! model time (minutes)
2832	REAL, INTENT(IN) :: dtmin ! model timestep (minutes)
2833	REAL, INTENT(IN) :: twindo ! half window (hours)
2834	REAL, INTENT(IN) :: scalef ! window scale factor
2835	REAL, INTENT(IN) :: obtime ! observation time (hours)
2836
2837	! Local variables
2838	real :: fdtim ! reference time (minutes)
2839	real :: tw1 ! half of twindo, scaled, in minutes
2840	real :: tw2 ! twindo, scaled, in minutes
2841	real :: tconst ! reciprical of tw1
2842	real :: ttim ! obtime in minutes
2843	real :: dift ! \| fdtim-ttim \|
2844	real :: timewt ! returned weight
2845
2846	! DETERMINE THE TIME-WEIGHT FACTOR FOR N
2847	FDTIM=XTIME-DTMIN
2848	! TWINDO IS IN MINUTES:
2849	TW1=TWINDO/2.60.scalef
2850	TW2=TWINDO60.scalef
2851	TCONST=1./TW1
2852	TIMEWT=0.0
2853	TTIM=obtime*60.
2854	!***********TTIM=TARGET TIME IN MINUTES
2855	DIFT=ABS(FDTIM-TTIM)
2856	IF(DIFT.LE.TW1)TIMEWT=1.0
2857	IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
2858	IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
2859	IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
2860	ENDIF
2861	get_timewt = timewt
2862	END FUNCTION get_timewt
2863
2864	SUBROUTINE print_vif_var(var, vif, nfullmin, nrampmin )
2865	!********************************************************
2866	! Purpose: Print a description of the vertical influence
2867	! function for a given variable.
2868	!********************************************************
2869	IMPLICIT NONE
2870
2871	character(len=4), intent(in) :: var ! Variable (wind, temp, mois)
2872	real, intent(in) :: vif(6) ! Vertical influence function
2873	real, intent(in) :: nfullmin ! Vert infl fcn full nudge min
2874	real, intent(in) :: nrampmin ! Vert infl fcn ramp decr min
2875
2876	! Local variables
2877	character(len=200) :: msg1, msg2
2878	character(len=8) :: regime
2879	real :: nfullr1, nrampr1
2880	real :: nfullr2, nrampr2
2881	real :: nfullr4, nrampr4
2882
2883	nfullr1 = vif(1)
2884	nrampr1 = vif(2)
2885	nfullr2 = vif(3)
2886	nrampr2 = vif(4)
2887	nfullr4 = vif(5)
2888	nrampr4 = vif(6)
2889
2890	if(var.eq.'wind') then
2891	write(msg1,fmt='(a)') ' For winds:'
2892	elseif (var.eq.'temp') then
2893	write(msg1,fmt='(a)') ' For temperature:'
2894	elseif (var.eq.'mois') then
2895	write(msg1,fmt='(a)') ' For moisture:'
2896	else
2897	write(msg1,fmt='(a,a4)') 'Unknown variable type: ',var
2898	call wrf_message(msg1)
2899	call wrf_error_fatal ( 'print_vif_var: module_fddaobs_rtfdda STOP' )
2900	endif
2901
2902	call wrf_message(msg1)
2903
2904	! For this variable, print a description of the vif for each regime
2905	call print_vif_regime(1, nfullr1, nrampr1, nfullmin, nrampmin)
2906	call print_vif_regime(2, nfullr2, nrampr2, nfullmin, nrampmin)
2907	call print_vif_regime(4, nfullr4, nrampr4, nfullmin, nrampmin)
2908
2909	END SUBROUTINE print_vif_var
2910
2911	SUBROUTINE print_vif_regime(reg, nfullr, nrampr, nfullmin, nrampmin )
2912	!********************************************************
2913	! Purpose: Print a description of the vertical influence
2914	! function for a given regime.
2915	!********************************************************
2916	IMPLICIT NONE
2917
2918	integer, intent(in) :: reg ! Regime number (1, 2, 4)
2919	real, intent(in) :: nfullr ! Full nudge range for regime
2920	real, intent(in) :: nrampr ! Rampdown range for regime
2921	real, intent(in) :: nfullmin ! Vert infl fcn full nudge min
2922	real, intent(in) :: nrampmin ! Vert infl fcn ramp decr min
2923
2924	! Local variables
2925	character(len=200) :: msg1, msg2
2926	character(len=8) :: regime
2927
2928	if(reg.eq.1) then
2929	write(regime,fmt='(a)') 'Regime 1'
2930	elseif (reg.eq.2) then
2931	write(regime,fmt='(a)') 'Regime 2'
2932	elseif (reg.eq.4) then
2933	write(regime,fmt='(a)') 'Regime 4'
2934	else
2935	write(msg1,fmt='(a,i3)') 'Unknown regime number: ',reg
2936	call wrf_message(msg1)
2937	call wrf_error_fatal ( 'print_vif_regime: module_fddaobs_rtfdda STOP' )
2938	endif
2939
2940	!Set msg1 for description of full weighting range
2941	if(nfullr.lt.0) then
2942	if(nfullr.eq.-5000) then
2943	write(msg1,fmt='(2x,a8,a)') regime, ': Full weighting to the PBL top'
2944	elseif (nfullr.lt.-5000) then
2945	write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(-5000.-nfullr), &
2946	' m above the PBL top'
2947	else
2948	write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(nfullr+5000.), &
2949	' m below the PBL top'
2950	endif
2951	else
2952	write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting through ', &
2953	int(max(nfullr,nfullmin)),' m'
2954	endif
2955
2956	!Set msg2 for description of rampdown range
2957	if(nrampr.lt.0) then
2958	if(nrampr.eq.-5000) then
2959	write(msg2,fmt='(a)') ' and a vertical rampdown up to the PBL top.'
2960	elseif (nrampr.lt.-5000) then
2961	write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(-5000.-nrampr), &
2962	' m above the PBL top.'
2963	else
2964	write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(nrampr+5000.), &
2965	' m below the PBL top.'
2966	endif
2967	else
2968	write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown in the next ', &
2969	int(max(nrampr,nrampmin)),' m.'
2970	endif
2971	call wrf_message(TRIM(msg1)//msg2)
2972
2973	END SUBROUTINE print_vif_regime
2974	#endif
2975
2976	END MODULE module_fddaobs_rtfdda
2977

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

Download in other formats:

Original Format