| blob | a242981e5c2b139f67edcd8a2a1e45816ce57c78 |
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE module_fddaobs_rtfdda
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.
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
38 CONTAINS
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 !-----------------------------------------------------------------------
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
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
168 #if ( EM_CORE == 1 )
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
175 call wrf_message("")
176 write(msg,fmt='(a,i2)') ' OBSERVATION NUDGING IS ACTIVATED FOR MESH ',inest
177 call wrf_message(msg)
179 ktau = itimestep
180 if(restart) then
181 fdob%ktaur = ktau
182 else
183 fdob%ktaur = 0
184 endif
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)
190 ! Set flag for nudging on pressure (not sigma) surfaces.
191 fdob%iwtsig = 0
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
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 !**************************************************************************
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
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
247 ! Store temporal and spatial scales
248 fdob%sfcfact = sfcfact
249 fdob%sfcfacr = sfcfacr
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("")
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
286 ! Initialize flags.
288 fdob%domain_tot=0
289 do nest=1,maxdom
290 fdob%domain_tot = fdob%domain_tot + nudge_opt(nest)
291 end do
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
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 )
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.
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
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("")
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
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
391 fdob%vif_fullmin = nudgezfullmin
392 fdob%vif_rampmin = nudgezrampmin
393 fdob%vif_max = nudgezmax
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.
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
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)
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
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
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)
457 ! Set parameters.
458 fdob%pfree = 50.0
459 fdob%rinfmn = 1.0
460 fdob%rinfmx = 2.0
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)
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
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
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
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
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
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 )
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
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 !------------------------------------------------------------------------------
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
571 !-----------------------------------------------------------------------
572 !
573 ! Description of input arguments.
574 !
575 !-----------------------------------------------------------------------
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.
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)
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.
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
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
656 real :: dummy = 99999.
658 real :: pbhi, pphi
659 real :: obs_pottemp ! Potential temperature at observation
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
671 real gridx,gridy,dxob,dyob,dzob,dxob_ms,dyob_ms
672 real pob
673 real hob
674 real uratiob,vratiob,tratiob,tratxob,fnpf
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
683 ! LOGICAL, EXTERNAL :: TILE_MASK
685 NSTA=NSTAT
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
691 ! ITYP=1 FOR U-POINTS, ITYP=2 for V-POINTS, and ITYP=3 FOR MASS POINTS
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
699 ERRF = 0.0 ! Zero out errf array
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 )
709 DO ITYP=1,3 ! Big loop: ityp=1 for U, ityp=2 for V, ityp=3 for T,Q,SP
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
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
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))
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)
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,
765 pob = 0.0
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)
771 IF ( MP_LOCAL_DUMMASK(N) ) THEN
772 #endif
774 ! Interpolate pressure to obs location column and convert from Pa to cb.
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
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)
795 if(abs(rko(n)+99).lt.1.)then
796 pob = varobs(5,n)
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
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
816 pphi = ppbo(kbot+1)
817 pbhi = pbbo(kbot+1)
819 rko(n) = real(kbot+1)- &
820 ( (pob-pbhi-pphi) / (pbbo(kbot)+ppbo(kbot)-pbhi-pphi) )
822 rko(n)=max(rko(n),1.0)
823 endif
824 endif
826 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
827 ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
828 #endif
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
835 RC(N)=RKO(N)
837 ENDDO ! END COARSE MESH LOOP OVER NSTA
839 ELSE ! FINE MESH CASE
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))
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)
872 ! ajb 20090423: BUGFIX TO OBS_IN_HEIGHT OPTION (see COARSE MESH comments)
873 pob = 0.0
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)
879 IF ( MP_LOCAL_DUMMASK(N) ) THEN
880 #endif
882 ! Interpolate pressure to obs location column and convert from Pa to cb.
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
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)
903 if(abs(rko(n)+99).lt.1.)then
904 pob = varobs(5,n)
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
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
924 pphi = ppbo(kbot+1)
925 pbhi = pbbo(kbot+1)
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
933 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
934 ENDIF !end IF( MP_LOCAL_DUMMASK(N) ) !ajb
935 #endif
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
942 RC(N)=RKO(N)
944 ENDDO ! END FINE MESH LOOP OVER NSTA
946 ENDIF ! end if(inest.eq.1)
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.
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
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.
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
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))
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)
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) )
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
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
1090 ! U COMPONENT WIND ERROR
1091 ERRF(1,N)=ERRF(1,N)+uratiob*VAROBS(1,N)-((1.-DZOB)* &
1092 ((1.-DyOB)*((1.- &
1093 DxOB)*UB(IOB,KOB,JOB)+DxOB*UB(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)+DxOB*UB(IOB+1,KOBP,JOB))+ &
1097 DyOB*((1.-DxOB)*UB(IOB,KOBP,JOB+1)+DxOB* &
1098 UB(IOB+1,KOBP,JOB+1))))
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
1121 ! Store model surface pressure (not the error!) at U point.
1122 ERRF(7,N)=.001*( (1.-DXOB)*PUG1 + DXOB*PUG2 )
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
1129 ENDIF ! end if(iswind.eq.1)
1131 ENDIF ! ITYP=1: PROCESS U
1133 !**********************************************************
1134 ! PROCESS V VARIABLE (IVAR=2)
1135 !**********************************************************
1136 IF(ITYP.EQ.2) THEN
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))
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)
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) )
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
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
1180 ! V COMPONENT WIND ERROR
1181 ERRF(2,N)=ERRF(2,N)+vratiob*VAROBS(2,N)-((1.-DZOB)* &
1182 ((1.-DyOB)*((1.- &
1183 DxOB)*VB(IOB,KOB,JOB)+DxOB*VB(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)+DxOB*VB(IOB+1,KOBP,JOB))+ &
1187 DyOB*((1.-DxOB)*VB(IOB,KOBP,JOB+1)+DxOB* &
1188 VB(IOB+1,KOBP,JOB+1))))
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
1203 ! Store model surface pressure (not the error!) at V point.
1204 ERRF(8,N)=.001*( (1.-DYOB)*PVG1 + DYOB*PVG2 )
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
1211 ENDIF ! end if(iswind.eq.1)
1213 ENDIF ! ITYP=1: PROCESS V
1215 !**********************************************************
1216 ! PROCESS MASS-POINT VARIABLES IVAR=3 (T) AND IVAR=4 (QV)
1217 !**********************************************************
1218 IF(ITYP.EQ.3) THEN
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)
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
1245 !yliu
1246 if(abs(tratxob) .lt. 1.0E-3) tratxob=1.
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
1254 ERRF(3,N)=ERRF(3,N)+tratxob*obs_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)))))
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
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))))
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) ))
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
1307 ENDIF !end if(istemp.eq.1 .or. ismois.eq.1)
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))
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)
1331 ENDDO ! END BIG LOOP
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)
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
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)
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
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
1377 ! For printout
1378 ! IF(MOD(KTAU,INPF).NE.0) THEN
1379 ! RETURN
1380 ! ENDIF
1382 RETURN
1383 END SUBROUTINE errob
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
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
1407 ! Local variables
1408 integer :: i, j
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
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
1429 RETURN
1430 END SUBROUTINE upoint
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
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
1454 ! Local variables
1455 integer :: i, j
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
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
1476 RETURN
1477 END SUBROUTINE vpoint
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
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
1497 ! Local variables
1498 LOGICAL :: retval
1500 TILE_MASK = (iloc .LE. ite .AND. iloc .GE. its .AND. &
1501 jloc .LE. jte .AND. jloc .GE. jts )
1503 RETURN
1504 END FUNCTION TILE_MASK
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
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 !-----------------------------------------------------------------------
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
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
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
1685 real :: scratch
1686 integer :: kk !ajb temp
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
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
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)
1725 IPL=IVAR + IPLO !yliu +IPLO
1727 ! DEFINE GRID-TYPE OFFSET FACTORS, IGRID AND GRID
1729 INPF=(IRATIO**LEVIDN(INEST))*NPFI
1730 INFR=(IRATIO**LEVIDN(INEST))*IONF
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
1749 ! TRANSFORM THE HORIZONTAL RADIUS OF INFLUENCE, RINXY, FROM
1750 ! KILOMETERS TO GRID LENGTHS, RINDX
1752 RINDX=RINXY*1000./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)
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
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
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
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
1809 ! DETERMINE THE LIMITS OF THE SEARCH REGION FOR THE CURRENT
1810 ! J-STRIP
1812 MAXJ=J+IFIX(RINDX*RINFMX+0.99) !ajb
1813 MINJ=J-IFIX(RINDX*RINFMX+0.99) !ajb
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
1819 n=1
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.
1828 !yliu
1829 ! note for sfc data: nsndlev=1 and njcsnd=1
1830 nsndlev=int(nlevs_ob(n)-lev_in_ob(n))+1
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)
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.
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
1851 IF( ABS(ERRF(IVAR,N)).GT.9.E4 .and. njcsnd.eq.1 ) THEN
1852 ! print *, " bad obs skipped"
1854 ELSEIF( RB(N).LT.FLOAT(MINJ) .OR. RB(N).GT.FLOAT(MAXJ) ) THEN
1855 ! print *, " skipped obs far away from this J-slice"
1857 !----------------------------------------------------------------------
1858 ELSE ! BEGIN SECTION FOR PROCESSING THE OBSERVATION
1859 !----------------------------------------------------------------------
1861 ! DETERMINE THE LIMITS OF APPLICATION OF THE OBS IN THE VERTICAL
1862 ! FOR THE VERTICAL WEIGHTING, WTSIG
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).
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)
1877 ! ASSIMILATE SURFACE LAYER DATA ON SIGMA
1878 IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5) THEN
1880 ! Compute temporal weight
1881 timewt = get_timewt(xtime,dtmin,twindo,sfcfact,timeob(n))
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
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)
1930 !ajb Compute and add weights to sum only if nudge_pbl switch is on.
1931 if(nudge_pbl) then
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=RINDX*RINDX*sfcfacr*sfcfacr
1939 RSQ=RX*RX+RY*RY
1940 ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
1941 DPRIM=SQRT(RSQ)
1942 D=DPRIM+RINDX*DCON*ABS(psurf(n)-.001*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=RINDX*RINDX*sfcfacr*sfcfacr
1952 RSQ=RX*RX+RY*RY
1953 ! THIS FUNCTION DECREASES WTIJ AS PSFC CHANGES WITHIN SEARCH RADIUS
1954 wtij(i)=(ris-rsq)/(ris+rsq)
1955 scratch = (abs (psurf(n)-.001*pbase(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
1962 ! Apply selected vertical spreading scheme.
1963 if(SFC_SCHEME_VERT.eq.1) then ! original simple scheme
1965 DO I=max0(its,MINI),min0(ite,MAXI)
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
1975 if(kpblt(i).gt.25) komax=3
1976 komin=1
1977 dk=float(komax)
1979 do k=komin,komax
1981 wtsig(k)=float(komax-k+1)/dk
1982 WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ(i)
1984 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*WTSIG(K) &
1985 *WTSIG(K)*ERFIVR
1986 enddo
1987 ENDDO
1989 else ! regime-based vif scheme
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
2017 DO I=max0(its,MINI),min0(ite,MAXI)
2019 dz_ramp = 1.0 / max( 1.0, vih2(i)-vih1(i) ) ! vih2 >= vih1 by construct
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))
2025 if(wtsig(k).le.0.0) EXIT LML
2026 WT(I,K)=WT(I,K)+TIMEWT*WTSIG(K)*WTIJ(i)
2027 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)*WTSIG(K) &
2028 *WTSIG(K)*ERFIVR
2029 enddo LML
2030 ENDDO
2031 endif
2033 endif ! end if nudge-pbl switch is on
2035 endif ! end check for obs in domain
2036 ! END SURFACE-LAYER OBS NUDGING
2038 ELSE
2040 ! Compute temporal weight
2041 timewt = get_timewt(xtime,dtmin,twindo,1.,timeob(n))
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)
2053 slope = (RINFMN-RINFMX)/(psurf(n)-PFREE)
2055 RINFAC=SLOPE*POB+RINFMX-SLOPE*pfree
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
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)
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 !!
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"
2093 ! yliu end
2094 ! is this 2 needed here - kpbl not used?
2095 ! MINI=MAX0(2,MINI)
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=RINDX*RINFAC*RINDX*RINFAC
2118 RSQ=RX*RX+RY*RY
2119 ! yliu test: for upper-air data, keep D1 influence radii
2120 WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
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
2133 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2134 ! --- HANDLE 1-LEVEL and MULTI-LEVEL OBSERVATIONS SEPARATELY ---
2135 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2137 if(nsndlev.eq.1)then
2138 !----------------------------------------------------------------------
2139 ! --- HANDLE 1-LEVEL OBSERVATIONS ---
2140 !----------------------------------------------------------------------
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
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
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
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
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)+TIMEWT*WTIJ(i)*wtsig(k)
2212 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)* &
2213 reserf(k)*wtsig(k)*wtsig(k)
2214 enddo
2215 endif
2217 else
2219 !----------------------------------------------------------------------
2220 ! --- HANDLE MULTI-LEVEL OBSERVATIONS ---
2221 !----------------------------------------------------------------------
2222 !yliu start
2223 ! if(I.eq.MINI) print *, " Multi-level snd "
2224 ! print *, " n=,nsndlev",n,nsndlev,nlevs_ob(n),lev_in_ob(n) &
2225 ! ,nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
2226 if(nlevs_ob(n+nsndlev-1).ne.lev_in_ob(n+nsndlev-1)) then
2227 IF (iprt) THEN
2228 write(msg,*) "n = ",n,"nsndlev = ",nsndlev
2229 call wrf_message(msg)
2230 write(msg,*) "nlevs_ob,lev_in_ob", &
2231 nlevs_ob(n+nsndlev-1),lev_in_ob(n+nsndlev-1)
2232 call wrf_message(msg)
2233 call wrf_message("in nudobs.F: sounding level messed up, stopping")
2234 ENDIF
2235 call wrf_error_fatal ( 'wrf_fddaobs_in: in4dob' )
2236 endif
2237 !yliu end
2238 ! This is for a multi-level observation
2239 ! The trick here is that the sounding is "one ob". You don't
2240 ! want multiple levels to each be treated like separate
2241 ! and independent observations.
2242 ! At each i,j want to interpolate sounding to the model levels at that
2243 ! particular point.
2244 komin=1
2245 komax=kte-2
2246 ! this loop goes to 1501
2247 ! do from kte-2 to 1 so don't adjust top of model. Arbitrary.
2248 !yliu start
2249 do k=1,kte
2250 reserf(k)=0.0
2251 wtsig(k)=0.0
2252 enddo
2253 !yliu end
2255 do k=komax,komin,-1
2257 pijk = .001*(pbase(i,k)+pp(i,k))
2259 ! if sigma level pressure is .gt. than the lowest ob level, don't interpolate
2260 if(pijk.gt.varobs(5,n)) then
2261 go to 1501
2262 endif
2264 ! if sigma level pressure is .lt. than the highest ob level, don't interpolate
2265 if(pijk.le.varobs(5,n+nsndlev-1)) then
2266 go to 1501
2267 endif
2269 ! now interpolate sounding to this k
2270 ! yliu start-- recalculate WTij for each k-level
2271 !ajb SLOPE = (RINFMN-RINFMX)/(pdoc(i,j)+PTOP-PFREE)
2272 slope = (RINFMN-RINFMX)/ (.001*pbase(i,1)-PFREE)
2273 RINFAC=SLOPE*pijk+RINFMX-SLOPE*PFREE
2274 RINFAC=AMAX1(RINFAC,RINFMN)
2275 RINFAC=AMIN1(RINFAC,RINFMX)
2276 RIS=RINDX*RINFAC*RINDX*RINFAC
2277 RSQ=RX*RX+RY*RY
2279 ! for upper-air data, keep D1 influence radii
2280 WTIJ(i)=(RIS-RSQ)/(RIS+RSQ)
2281 WTIJ(i)=AMAX1(0.0,WTIJ(i))
2282 ! yliu end
2284 ! this loop goes to 1503
2285 do nn=2,nsndlev
2287 ! only set pobhi if varobs(ivar) is ok
2288 pobhi=-888888.
2290 if(varobs(ivar,n+nn-1).gt.-800000. &
2291 .and. varobs(5,n+nn-1).gt.-800000.) then
2292 pobhi=varobs(5,n+nn-1)
2293 nhi=n+nn-1
2294 if(pobhi.lt.pijk .and. abs(pobhi-pijk).lt.0.5*maxsnd_gap) then
2295 go to 1502 ! within maxsnd_gap/2 of obs height
2296 endif
2297 endif
2299 enddo
2301 ! did not find any ob above within maxsnd_gap/2, so jump out
2302 go to 1501
2303 1502 continue
2305 nlo=nhi-1
2306 do nnjc=nhi-1,n,-1
2307 if(varobs(ivar,nnjc).gt.-800000. &
2308 .and. varobs(5,nnjc).gt.-800000.) then
2309 poblo=varobs(5,nnjc)
2310 nlo=nnjc
2311 if(poblo.gt.pijk .and. abs(poblo-pijk).lt.0.5*maxsnd_gap) then
2312 go to 1505 ! within maxsnd_gap/2 of obs height
2313 endif
2314 endif
2315 enddo
2316 !yliu end --
2318 ! did not find any ob below within maxsnd_gap/2, so jump out
2319 go to 1501
2320 1505 continue
2322 ! interpolate to model level
2323 pdiffj=alog(pijk/poblo)/alog(pobhi/poblo)
2324 reserf(k)=errf(ivar,nlo)+ &
2325 (errf(ivar,nhi)-errf(ivar,nlo))*pdiffj
2326 wtsig(k)=1.
2328 1501 continue
2330 ! now calculate WT and WT2ERR for each i,j,k point cajb
2331 !ajb Add weights to sum only if nudge_pbl switch is on OR k > kpblt.
2332 if(nudge_pbl .or. k.gt.kpblt(i)) then
2334 WT(I,K)=WT(I,K)+TIMEWT*WTIJ(i)*wtsig(k)
2336 WT2ERR(I,K)=WT2ERR(I,K)+TIMEWT*TIMEWT*WTIJ(i)*WTIJ(i)* &
2337 reserf(k)*wtsig(k)*wtsig(k)
2338 endif
2340 enddo ! enddo k levels
2342 ! end multi-levels
2343 endif ! end if(nsndlev.eq.1)
2344 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2345 ! END 1-LEVEL AND MULTI-LEVEL OBSERVATIONS
2346 !$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
2347 !
2348 ENDDO ! END DO MINI,MAXI LOOP
2350 endif ! check for obs in domain
2352 ! END OF NUDGING TO OBS ON PRESSURE LEVELS
2354 ENDIF !end IF(KOB.EQ.1.AND.IVAR.LE.4.and.nlevs_ob(n).lt.1.5)
2356 !----------------------------------------------------------------------
2357 ENDIF ! END SECTION FOR PROCESSING OF OBSERVATION
2358 !----------------------------------------------------------------------
2360 ! n=n+1
2361 n=n+njcsnd
2363 !yliu 1202 continue
2364 if(n.gt.nstat)then
2365 ! print *,'n,nstat=',n,nstat,ivar,j
2366 go to 1203
2367 endif
2368 ! print *, "e-- n=,nsndlev",n,njcsnd,nlevs_ob(n),lev_in_ob(n)
2370 !***********************************************************************
2371 ENDDO ! END OUTER LOOP FOR THE NSTAT OBSERVATIONS
2372 !***********************************************************************
2374 1203 continue
2376 ! WEIGHTS AND WEIGHTED DIFFERENCES HAVE BEEN SUMMED. NOW
2377 ! APPLY THE NUDGING FACTOR AND THE RESULTANT TENDENCY TO
2378 ! THE ATEN ARRAY
2379 ! ASSURE THAT WT(I,K) AND WTP(I,K) ARE NONZERO SINCE
2380 ! THEY ARE USED BELOW IN THE DENOMINATOR.
2381 DO K=kts,kte
2382 DO I=its,ite
2383 IF(WT(I,K).EQ.0)THEN
2384 WT2ERR(I,K)=0.0
2385 ENDIF
2386 IF(WT(I,K).EQ.0)THEN
2387 WT(I,K)=1.0
2388 ENDIF
2389 ENDDO
2390 ENDDO
2392 126 CONTINUE
2394 IF(IVAR.GE.3)GOTO 170
2395 ! this is for u,v
2396 ! 3-D DOT POINT TENDENCIES
2398 ! Calculate scales for converting nudge factor from u (v)
2399 ! to rho_u (or rho_v) units.
2401 IF (IVAR == 1) THEN
2402 call calc_rcouple_scales(mu,msfy,rscale,ims,ime,its,ite)
2403 ELSE IF (IVAR == 2) THEN
2404 call calc_rcouple_scales(mu,msfx,rscale,ims,ime,its,ite)
2405 END IF
2407 DO K=1,kte
2409 DO I=i_s,i_e
2411 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2412 W2EOWT=WT2ERR(I,K)/WT(I,K)
2413 ELSE
2414 W2EOWT=SAVWT(IPL,I,K)
2416 ! write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,W2EOWT = ',i,j,k,ipl,W2EOWT
2418 ENDIF
2420 ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2421 ! scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
2422 ! write(6,*) 'ATEN calc: k = ',k
2423 ! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
2424 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2425 ! $ ' W2EOWT = ',w2eowt
2426 ! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
2427 ! $ ' GFACTOR = ',gfactor
2428 ! endif
2429 !
2430 ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2431 ! scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
2432 ! write(6,*) 'ATEN calc: k = ',k
2433 ! write(6,*) 'V before: aten = ',aten(i,k),' scr = ',scratch
2434 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),
2435 ! $ ' W2EOWT = ',w2eowt
2436 ! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,
2437 ! $ ' GFACTOR = ',gfactor
2438 ! endif
2440 ! if(ivar .eq. 1 .and. i.eq.38 .and. (j.ge.36.and.j.le.62) .and. k.eq.7) then
2441 ! scratch = GIV*RSCALE(I)*W2EOWT*TFACI*ISWIND*GFACTOR
2442 ! write(6,*)
2443 ! write(6,*) 'ATEN calc: j = ',j,' k = ',k
2444 ! write(6,*) 'U before: aten = ',aten(i,k),' scr = ',scratch
2445 ! write(6,*) 'GIV = ',giv,' rscale = ',rscale(i),' W2EOWT = ',w2eowt
2446 ! write(6,*) 'TFACI = ',tfaci,' ISWIND = ',iswind,' GFACTOR = ',gfactor
2447 ! endif
2449 ATEN(i,k)=ATEN(i,k)+GIV*RSCALE(I) &
2450 *W2EOWT*TFACI &
2451 *ISWIND *GFACTOR !yliu *GFACTOR
2453 ! if(ivar .eq. 1 .and. i.eq.38 .and. j.eq.78 .and. k.eq.1) then
2454 ! write(6,*) 'U after: aten = ',aten(i,k),' scr = ',scratch
2455 ! endif
2456 ! if(ivar .eq. 2 .and. i.eq.39 .and. j.eq.29) then
2457 ! write(6,*) 'V after: aten = ',aten(i,k),' scr = ',scratch
2458 ! endif
2460 ENDDO
2461 ENDDO
2463 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2464 DO K=1,kte
2465 DO I=its,ite
2466 SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2468 ! write(6,'(a,4i4,f8.3)') 'i,j,k,ipl,savwt = ',i,j,k,ipl,savwt(ipl,i,k)
2469 ENDDO
2470 ENDDO
2471 ENDIF
2473 RETURN
2475 170 CONTINUE
2477 ! 3-D CROSS-POINT TENDENCIES
2478 ! this is for t (ivar=3) and q (ivsr=4)
2479 IF(3-IVAR.LT.0)THEN
2480 GITQ=GIQ
2481 ELSE
2482 GITQ=GIT
2483 ENDIF
2484 IF(3-IVAR.LT.0)THEN
2485 ISTQ=ISMOIS
2486 ELSE
2487 ISTQ=ISTEMP
2488 ENDIF
2490 DO K=1,kte
2491 DO I=i_s,i_e
2492 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR))THEN
2493 W2EOWT=WT2ERR(I,K)/WT(I,K)
2494 ELSE
2495 W2EOWT=SAVWT(IPL,I,K)
2496 ENDIF
2498 ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
2499 ! scratch = GITQ*MU(I)*W2EOWT*TFACI*ISTQ*GFACTOR
2500 ! write(6,*) 'ATEN calc: k = ',k
2501 ! write(6,*) 'T before: aten = ',aten(i,k),' scr = ',scratch
2502 ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i), &
2503 ! ' W2EOWT = ',w2eowt
2504 ! write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq, &
2505 ! ' GFACTOR = ',gfactor
2506 ! endif
2508 ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2509 ! scratch = GITQ*MU(I)*W2EOWT*TFACI*ISTQ*GFACTOR
2510 ! write(6,*) 'ATEN calc: k = ',k
2511 ! write(6,*) 'Q before: aten = ',aten(i,k),' scr = ',scratch
2512 ! write(6,*) 'GITQ = ',gitq,' MU = ',mu(i),
2513 ! $ ' W2EOWT = ',w2eowt
2514 ! write(6,*) ' TFACI = ',tfaci,' ISTQ = ',istq,
2515 ! $ ' GFACTOR = ',gfactor
2516 ! endif
2518 ATEN(i,k)=ATEN(i,k)+GITQ*MU(I) &
2519 *W2EOWT*TFACI*ISTQ *GFACTOR !yliu *GFACTOR
2521 ! if(ivar .eq. 3 .and. i.eq.39 .and. j.eq.49) then
2522 ! write(6,*) 'T after: aten = ',aten(i,k),' scr = ',scratch
2523 ! endif
2524 ! if(ivar .eq. 4 .and. i.eq.39 .and. j.eq.29) then
2525 ! write(6,*) 'Q after: aten = ',aten(i,k),' scr = ',scratch
2526 ! endif
2528 ENDDO
2529 ENDDO
2531 IF(MOD(KTAU,INFR).EQ.0.OR.(IFREST.AND.KTAU.EQ.KTAUR)) THEN
2532 DO K=1,kte
2533 DO I=its,ite
2534 SAVWT(IPL,I,K)=WT2ERR(I,K)/WT(I,K)
2535 ENDDO
2536 ENDDO
2537 ENDIF
2539 RETURN
2540 END SUBROUTINE nudob
2542 SUBROUTINE date_string(year, month, day, hour, minute, second, cdate)
2543 !-----------------------------------------------------------------------
2544 ! PURPOSE: Form a date string (YYYY-MM-DD_hh:mm:ss) from integer
2545 ! components.
2546 !-----------------------------------------------------------------------
2547 IMPLICIT NONE
2548 !-----------------------------------------------------------------------
2550 INTEGER, INTENT(IN) :: year
2551 INTEGER, INTENT(IN) :: month
2552 INTEGER, INTENT(IN) :: day
2553 INTEGER, INTENT(IN) :: hour
2554 INTEGER, INTENT(IN) :: minute
2555 INTEGER, INTENT(IN) :: second
2556 CHARACTER*19, INTENT(INOUT) :: cdate
2558 ! Local variables
2559 integer :: ic ! loop counter
2561 cdate(1:19) = "0000-00-00_00:00:00"
2562 write(cdate( 1: 4),'(i4)') year
2563 write(cdate( 6: 7),'(i2)') month
2564 write(cdate( 9:10),'(i2)') day
2565 write(cdate(12:13),'(i2)') hour
2566 write(cdate(15:16),'(i2)') minute
2567 write(cdate(18:19),'(i2)') second
2568 do ic = 1,19
2569 if(cdate(ic:ic) .eq. " ") cdate(ic:ic) = "0"
2570 enddo
2572 RETURN
2573 END SUBROUTINE date_string
2575 SUBROUTINE calc_rcouple_scales(a, msf, rscale, ims,ime, its,ite)
2576 !-----------------------------------------------------------------------
2577 IMPLICIT NONE
2578 !-----------------------------------------------------------------------
2580 INTEGER, INTENT(IN) :: ims,ime ! Memory dimensions
2581 INTEGER, INTENT(IN) :: its,ite ! Tile dimensions
2582 REAL, INTENT(IN) :: a( ims:ime ) ! Air mass array
2583 REAL, INTENT(IN) :: msf( ims:ime ) ! Map scale factor array
2584 REAL, INTENT(OUT) :: rscale( ims:ime ) ! Scales for rho-coupling
2586 ! Local variables
2587 integer :: i
2589 ! Calculate scales to be used for producing rho-coupled nudging factors.
2590 do i = its,ite
2591 rscale(i) = a(i)/msf(i)
2592 enddo
2594 RETURN
2595 END SUBROUTINE calc_rcouple_scales
2597 SUBROUTINE print_obs_info(iprt,inest,niobf,rio,rjo,rko, &
2598 prt_max,prt_freq,obs,stnid,lat,lon, &
2599 mlat,mlon,timeob,xtime)
2600 !*************************************************************************
2601 ! Purpose: Print obs information.
2602 !*************************************************************************
2604 IMPLICIT NONE
2606 LOGICAL, intent(in) :: iprt ! Print flag
2607 INTEGER, intent(in) :: inest ! Nest level
2608 INTEGER, intent(in) :: niobf ! Maximum number of observations
2609 REAL, intent(in) :: rio(niobf) ! West-east coord (non-stagger)
2610 REAL, intent(in) :: rjo(niobf) ! South-north coord (non-stagger)
2611 REAL, intent(in) :: rko(niobf) ! Bottom-top north coord (non-stagger)
2612 INTEGER, intent(in) :: prt_max ! Max no. of obs for diagnostic printout
2613 INTEGER, intent(in) :: prt_freq ! Frequency for diagnostic printout
2614 INTEGER, intent(in) :: obs(prt_max) ! Saved obs indices to print
2615 INTEGER, intent(in) :: stnid(40,prt_max) ! Saved station ids
2616 REAL, intent(in) :: lat(prt_max) ! Saved latitudes
2617 REAL, intent(in) :: lon(prt_max) ! Saved longitudes
2618 REAL, intent(in) :: mlat(prt_max) ! Saved model latitudes
2619 REAL, intent(in) :: mlon(prt_max) ! Saved longitudes
2620 REAL, intent(in) :: timeob(niobf) ! Times of each observation (hours)
2621 REAL, intent(in) :: xtime ! Model time in minutes
2623 ! Local variables
2624 integer :: i ! Loop counter over obs station chars
2625 integer :: n ! Loop counter over obs
2626 integer :: pnx ! Obs index for printout
2627 character(len=200) :: msg ! Argument to wrf_message
2628 character(len=20) :: station_id ! Station id of observation
2630 if(iprt) then
2631 if(prt_max.gt.0) then
2633 if(obs(1).ne.-999) then
2635 call wrf_message("")
2636 write(msg,fmt='(a,i4,a,f8.1,a)') 'REPORTING OBS MASS-PT LOCS FOR NEST ', &
2637 inest,' AT XTIME=',xtime,' MINUTES'
2638 call wrf_message(msg)
2640 write(msg,fmt='(a,i4,a,i5,a)') 'FREQ=',prt_freq,', MAX=',prt_max, &
2641 ' LOCS, NEWLY READ OBS ONLY, -999 => OBS OFF PROC'
2642 call wrf_message(msg)
2643 call wrf_message("")
2645 write(msg,fmt='(3a)') ' OBS# I J K OBS LAT', &
2646 ' OBS LON XLAT(I,J) XLONG(I,J) TIME(hrs)', &
2647 ' OBS STATION ID'
2648 call wrf_message(msg)
2650 endif
2651 endif
2653 ! Note: rio and rjo are referenced to non-staggered grid (not mass-point!)
2654 ! Hence subtract .5 from each to get mass-point coords.
2655 do n=1,prt_max
2656 pnx = obs(n)
2657 if(pnx.ne.-999) then
2658 ! Retrieve 15 chars of station id
2659 do i = 1,15
2660 station_id(i:i) = char(stnid(i,n))
2661 enddo
2662 write(msg,fmt='(2x,i7,3f8.3,2f9.3,2x,f9.3,2x,f9.3,3x,f6.2,7x,a15)') &
2663 pnx,rio(pnx)-.5,rjo(pnx)-.5,rko(pnx),lat(n),lon(n), &
2664 mlat(n),mlon(n),timeob(pnx),station_id
2665 call wrf_message(msg)
2666 endif
2667 enddo
2668 if(obs(1).ne.-999) call wrf_message("")
2669 endif
2670 END SUBROUTINE print_obs_info
2672 REAL FUNCTION ht_to_p( h, pbbc, ppbc, z, ic, jc, dx, dy, &
2673 k_start, k_end, kds,kde, ims,ime, jms,jme, kms,kme )
2675 !******************************************************************************
2676 ! Purpose: Interpolate pressure at a specified x (ic), y (jc), and height (h).
2677 ! The input pressure column pbbc+ppbc (base and perturbn) must already
2678 ! be horizontally interpolated to the x, y position. The subroutine
2679 ! get_height_column is called here to horizontally interpolated the
2680 ! 3D height field z to get a height column at (iob, job).
2681 !******************************************************************************
2683 IMPLICIT NONE
2685 REAL, INTENT(IN) :: h ! height value (m)
2686 INTEGER, INTENT(IN) :: k_start, k_end ! loop bounds
2687 INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
2688 INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
2689 REAL, INTENT(IN) :: pbbc(kds:kde) ! column base pressure (cb)
2690 REAL, INTENT(IN) :: ppbc(kds:kde) ! column pressure perturbation (cb)
2691 REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! ht (m) above sl on half-levels
2692 INTEGER, INTENT(IN) :: ic ! i-coord of desired p
2693 INTEGER, INTENT(IN) :: jc ! j-coord of desired p
2694 REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
2695 REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
2697 ! Local variables
2698 INTEGER :: k ! loop counter
2699 INTEGER :: klo ! lower k bound
2700 REAL :: zlo ! lower z bound for h
2701 REAL :: zhi ! upper z bound for h
2702 REAL :: p ! interpolated pressure value
2703 REAL :: ln_p ! log p
2704 REAL :: ln_plo ! log plo
2705 REAL :: ln_phi ! log phi
2706 REAL :: z_at_p( kms:kme ) ! height at p levels
2708 ! Get interpolated z column on pressure (half-) levels at (ic,jc)
2709 call get_height_column(z, ic, jc, dx, dy, z_at_p, &
2710 k_start, k_end, kds,kde, &
2711 ims,ime, jms,jme, kms,kme )
2713 ! Now we have pbbc, ppbc, z_at_p, so compute p at h. First, find
2714 ! bounding layers klo and khi so that z_at_p(klo) <= h <= z_at_p(khi)
2716 ZLEVS: do k = k_start+1, k_end
2717 klo = k-1
2718 if(h .le. z_at_p(k)) then
2719 EXIT ZLEVS
2720 endif
2721 enddo ZLEVS
2723 zlo = z_at_p(klo)
2724 zhi = z_at_p(klo+1)
2726 ! Interpolate natural log of pressure
2727 ln_plo = log( pbbc(klo+1) + ppbc(klo+1) )
2728 ln_phi = log( pbbc(klo) + ppbc(klo) )
2729 if(h.le.zlo) then
2730 ln_p = ln_phi ! set to k=1 pressure
2731 else if (h.ge.zhi) then
2732 ln_p = ln_plo ! set to k=k_end pressure
2733 else
2734 ln_p = ln_plo + (ln_phi-ln_plo)*((zhi-h)/(zhi-zlo))
2735 endif
2737 ! Return pressure
2738 p = exp(ln_p)
2739 ht_to_p = p
2740 RETURN
2741 END FUNCTION ht_to_p
2743 SUBROUTINE get_height_column( z, ic, jc, dx, dy, z_at_p, &
2744 k_start, k_end, kds,kde, &
2745 ims,ime, jms,jme, kms,kme )
2746 !*************************************************************************
2747 ! Purpose: Compute column height at ic, jc location on p levels
2748 !*************************************************************************
2750 IMPLICIT NONE
2752 INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
2753 INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
2754 INTEGER, INTENT(IN) :: ims,ime, jms,jme, kms,kme ! memory dims.
2755 REAL, INTENT(IN) :: z( ims:ime, kms:kme, jms:jme ) ! ht (m) on half-levels
2756 INTEGER, INTENT(IN) :: ic ! i-coord of desired p
2757 INTEGER, INTENT(IN) :: jc ! j-coord of desired p
2758 REAL, INTENT(IN) :: dx ! interp. fraction (x dir)
2759 REAL, INTENT(IN) :: dy ! interp. fraction (y dir)
2760 REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
2762 ! Local variables
2763 INTEGER :: k ! loop counter
2766 do k = kds, kde
2767 z_at_p(k) = &
2768 (1.-DY)*( (1.-DX)*z(IC,K,JC) + &
2769 DX *z(IC+1,K,JC) ) + &
2770 DY* ( (1.-DX)*z(IC,K,JC+1) + &
2771 DX *z(IC+1,K,JC+1) )
2772 enddo
2774 END SUBROUTINE get_height_column
2776 SUBROUTINE get_base_state_height_column( p_top, p00, t00, a, g, r_d, &
2777 znu, z_at_p, k_start, k_end, kds,kde, kms,kme )
2778 !********************************************************
2779 ! Purpose: Compute base-state column height on p levels.
2780 ! See, e.g., eqns 1.3-1.5 in MM5 User's Guide.
2781 ! Height is a function of pressure plus the constants:
2782 ! p00 - sea level pressure (Pa)
2783 ! t00 - sea level temperature (K)
2784 ! a - base state lapse rate (k/m)
2785 ! r_d - gas constant (J/Kg/K) (Joule=1kg/m**2/s**2)
2786 ! g - gravitational constant (m/s**2)
2787 !********************************************************
2789 IMPLICIT NONE
2791 REAL, INTENT(IN) :: p_top ! pressure at top of model
2792 REAL, INTENT(IN) :: p00 ! base state pressure
2793 REAL, INTENT(IN) :: t00 ! base state temperature
2794 REAL, INTENT(IN) :: a ! base state lapse rate
2795 REAL, INTENT(IN) :: g ! gravity constant
2796 REAL, INTENT(IN) :: r_d ! gas constant
2797 INTEGER, INTENT(IN) :: k_start, k_end ! Loop bounds
2798 INTEGER, INTENT(IN) :: kds,kde ! vertical dim.
2799 INTEGER, INTENT(IN) :: kms,kme ! vertical memory dim.
2800 REAL, INTENT(IN) :: znu( kms:kme ) ! eta values on half (mass) levels
2801 REAL, INTENT(OUT) :: z_at_p( kms:kme ) ! column height at p levels
2803 ! Local variables
2804 integer :: k ! loop counter
2805 real :: ps0 ! base state pressure at surface
2806 real :: pb(kms:kme) ! pressure on half eta levels
2807 real :: logterm ! temporary for Z calculation
2808 real :: ginv ! 1/g
2810 ginv = 1/g
2812 ! Compute base state pressure on half eta levels.
2813 do k = k_start, k_end
2814 pb(k) = znu(k)*(p00 - p_top) + p_top
2815 enddo
2817 ! Use hydrostatic relation to compute height at pressure levels.
2818 do k = k_start, k_end
2819 logterm = log(pb(k)/p00)
2820 z_at_p(k) = .5*r_d*a*ginv*logterm*logterm - r_d*t00*ginv*logterm
2821 enddo
2823 END SUBROUTINE get_base_state_height_column
2825 REAL FUNCTION get_timewt(xtime,dtmin,twindo,scalef,obtime)
2826 !*************************************************************************
2827 ! Purpose: Compute the temporal weight factor for an observation
2828 !*************************************************************************
2830 IMPLICIT NONE
2832 REAL, INTENT(IN) :: xtime ! model time (minutes)
2833 REAL, INTENT(IN) :: dtmin ! model timestep (minutes)
2834 REAL, INTENT(IN) :: twindo ! half window (hours)
2835 REAL, INTENT(IN) :: scalef ! window scale factor
2836 REAL, INTENT(IN) :: obtime ! observation time (hours)
2838 ! Local variables
2839 real :: fdtim ! reference time (minutes)
2840 real :: tw1 ! half of twindo, scaled, in minutes
2841 real :: tw2 ! twindo, scaled, in minutes
2842 real :: tconst ! reciprical of tw1
2843 real :: ttim ! obtime in minutes
2844 real :: dift ! | fdtim-ttim |
2845 real :: timewt ! returned weight
2847 ! DETERMINE THE TIME-WEIGHT FACTOR FOR N
2848 FDTIM=XTIME-DTMIN
2849 ! TWINDO IS IN MINUTES:
2850 TW1=TWINDO/2.*60.*scalef
2851 TW2=TWINDO*60.*scalef
2852 TCONST=1./TW1
2853 TIMEWT=0.0
2854 TTIM=obtime*60.
2855 !***********TTIM=TARGET TIME IN MINUTES
2856 DIFT=ABS(FDTIM-TTIM)
2857 IF(DIFT.LE.TW1)TIMEWT=1.0
2858 IF(DIFT.GT.TW1.AND.DIFT.LE.TW2) THEN
2859 IF(FDTIM.LT.TTIM)TIMEWT=(FDTIM-(TTIM-TW2))*TCONST
2860 IF(FDTIM.GT.TTIM)TIMEWT=((TTIM+TW2)-FDTIM)*TCONST
2861 ENDIF
2862 get_timewt = timewt
2863 END FUNCTION get_timewt
2865 SUBROUTINE print_vif_var(var, vif, nfullmin, nrampmin )
2866 !********************************************************
2867 ! Purpose: Print a description of the vertical influence
2868 ! function for a given variable.
2869 !********************************************************
2870 IMPLICIT NONE
2872 character(len=4), intent(in) :: var ! Variable (wind, temp, mois)
2873 real, intent(in) :: vif(6) ! Vertical influence function
2874 real, intent(in) :: nfullmin ! Vert infl fcn full nudge min
2875 real, intent(in) :: nrampmin ! Vert infl fcn ramp decr min
2877 ! Local variables
2878 character(len=200) :: msg1, msg2
2879 character(len=8) :: regime
2880 real :: nfullr1, nrampr1
2881 real :: nfullr2, nrampr2
2882 real :: nfullr4, nrampr4
2884 nfullr1 = vif(1)
2885 nrampr1 = vif(2)
2886 nfullr2 = vif(3)
2887 nrampr2 = vif(4)
2888 nfullr4 = vif(5)
2889 nrampr4 = vif(6)
2891 if(var.eq.'wind') then
2892 write(msg1,fmt='(a)') ' For winds:'
2893 elseif (var.eq.'temp') then
2894 write(msg1,fmt='(a)') ' For temperature:'
2895 elseif (var.eq.'mois') then
2896 write(msg1,fmt='(a)') ' For moisture:'
2897 else
2898 write(msg1,fmt='(a,a4)') 'Unknown variable type: ',var
2899 call wrf_message(msg1)
2900 call wrf_error_fatal ( 'print_vif_var: module_fddaobs_rtfdda STOP' )
2901 endif
2903 call wrf_message(msg1)
2905 ! For this variable, print a description of the vif for each regime
2906 call print_vif_regime(1, nfullr1, nrampr1, nfullmin, nrampmin)
2907 call print_vif_regime(2, nfullr2, nrampr2, nfullmin, nrampmin)
2908 call print_vif_regime(4, nfullr4, nrampr4, nfullmin, nrampmin)
2910 END SUBROUTINE print_vif_var
2912 SUBROUTINE print_vif_regime(reg, nfullr, nrampr, nfullmin, nrampmin )
2913 !********************************************************
2914 ! Purpose: Print a description of the vertical influence
2915 ! function for a given regime.
2916 !********************************************************
2917 IMPLICIT NONE
2919 integer, intent(in) :: reg ! Regime number (1, 2, 4)
2920 real, intent(in) :: nfullr ! Full nudge range for regime
2921 real, intent(in) :: nrampr ! Rampdown range for regime
2922 real, intent(in) :: nfullmin ! Vert infl fcn full nudge min
2923 real, intent(in) :: nrampmin ! Vert infl fcn ramp decr min
2925 ! Local variables
2926 character(len=200) :: msg1, msg2
2927 character(len=8) :: regime
2929 if(reg.eq.1) then
2930 write(regime,fmt='(a)') 'Regime 1'
2931 elseif (reg.eq.2) then
2932 write(regime,fmt='(a)') 'Regime 2'
2933 elseif (reg.eq.4) then
2934 write(regime,fmt='(a)') 'Regime 4'
2935 else
2936 write(msg1,fmt='(a,i3)') 'Unknown regime number: ',reg
2937 call wrf_message(msg1)
2938 call wrf_error_fatal ( 'print_vif_regime: module_fddaobs_rtfdda STOP' )
2939 endif
2941 !Set msg1 for description of full weighting range
2942 if(nfullr.lt.0) then
2943 if(nfullr.eq.-5000) then
2944 write(msg1,fmt='(2x,a8,a)') regime, ': Full weighting to the PBL top'
2945 elseif (nfullr.lt.-5000) then
2946 write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(-5000.-nfullr), &
2947 ' m above the PBL top'
2948 else
2949 write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting to ',int(nfullr+5000.), &
2950 ' m below the PBL top'
2951 endif
2952 else
2953 write(msg1,fmt='(2x,a8,a,i4,a)') regime, ': Full weighting through ', &
2954 int(max(nfullr,nfullmin)),' m'
2955 endif
2957 !Set msg2 for description of rampdown range
2958 if(nrampr.lt.0) then
2959 if(nrampr.eq.-5000) then
2960 write(msg2,fmt='(a)') ' and a vertical rampdown up to the PBL top.'
2961 elseif (nrampr.lt.-5000) then
2962 write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(-5000.-nrampr), &
2963 ' m above the PBL top.'
2964 else
2965 write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown to ',int(nrampr+5000.), &
2966 ' m below the PBL top.'
2967 endif
2968 else
2969 write(msg2,fmt='(a,i4,a)') ' and a vertical rampdown in the next ', &
2970 int(max(nrampr,nrampmin)),' m.'
2971 endif
2972 call wrf_message(TRIM(msg1)//msg2)
2974 END SUBROUTINE print_vif_regime
2975 #endif
2977 END MODULE module_fddaobs_rtfdda