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