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