| blob | 6a64a62ae2665d193348c1f10960d9a277ee9f15 |
1 !wrf:model_layer:physics
2 !
3 !ckay=KIRAN ALAPATY @ US EPA -- November 01, 2015
4 !
5 ! Tim Glotfelty@CNSU; AJ Deng@PSU
6 !modified for use with FASDAS
7 !Flux Adjusting Surface Data Assimilation System to assimilate
8 !surface layer and soil layers temperature and moisture using
9 ! surfance reanalsys
10 !Reference: Alapaty et al., 2008: Development of the flux-adjusting surface
11 ! data assimilation system for mesoscale models. JAMC, 47, 2331-2350
13 !
14 !
15 MODULE module_fdda_psufddagd
17 CONTAINS
18 !
19 !-------------------------------------------------------------------
20 !
21 SUBROUTINE fddagd(itimestep,dx,dt,xtime, &
22 id,analysis_interval, end_fdda_hour, &
23 if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
24 if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
25 guv, gt, gq, if_ramping, dtramp_min, &
27 grid_sfdda, &
28 analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
29 rinblw, &
31 u3d,v3d,th3d,t3d, &
32 qv3d, &
33 p3d,pi3d, &
34 u_ndg_old,v_ndg_old,t_ndg_old,q_ndg_old,mu_ndg_old, &
35 u_ndg_new,v_ndg_new,t_ndg_new,q_ndg_new,mu_ndg_new, &
36 u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
37 rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
38 u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
39 rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
40 RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
41 fasdas, SDA_HFX, SDA_QFX, & ! fasdas
42 HFX_FDDA,dz8w, & ! fasdas
43 pblh, ht, regime, znt, z, z_at_w, &
44 ids,ide, jds,jde, kds,kde, &
45 ims,ime, jms,jme, kms,kme, &
46 its,ite, jts,jte, kts,kte )
48 !-------------------------------------------------------------------
49 implicit none
50 !-------------------------------------------------------------------
51 !
52 ! This code was implemented by Aijun Deng (Penn State). The 3-D analysis nudging was comleted
53 ! and released in December 2006. The surface analysis nudging capability was added and
54 ! released in March 2009 with WRFV3.1.
55 !
56 !-- u3d 3d u-velocity staggered on u points
57 !-- v3d 3d v-velocity staggered on v points
58 !-- th3d 3d potential temperature (k)
59 !-- t3d temperature (k)
60 !-- qv3d 3d water vapor mixing ratio (kg/kg)
61 !-- p3d 3d pressure (pa)
62 !-- pi3d 3d exner function (dimensionless)
63 !-- rundgdten staggered u tendency due to
64 ! fdda grid nudging (m/s/s)
65 !-- rvndgdten staggered v tendency due to
66 ! fdda grid nudging (m/s/s)
67 !-- rthndgdten theta tendency due to
68 ! fdda grid nudging (K/s)
69 !-- rqvndgdten qv tendency due to
70 ! fdda grid nudging (kg/kg/s)
71 !-- rmundgdten mu tendency due to
72 ! fdda grid nudging (Pa/s)
73 !-- ids start index for i in domain
74 !-- ide end index for i in domain
75 !-- jds start index for j in domain
76 !-- jde end index for j in domain
77 !-- kds start index for k in domain
78 !-- kde end index for k in domain
79 !-- ims start index for i in memory
80 !-- ime end index for i in memory
81 !-- jms start index for j in memory
82 !-- jme end index for j in memory
83 !-- kms start index for k in memory
84 !-- kme end index for k in memory
85 !-- its start index for i in tile
86 !-- ite end index for i in tile
87 !-- jts start index for j in tile
88 !-- jte end index for j in tile
89 !-- kts start index for k in tile
90 !-- kte end index for k in tile
91 !-------------------------------------------------------------------
92 !
93 INTEGER, INTENT(IN) :: itimestep, analysis_interval, end_fdda_hour
94 INTEGER, INTENT(IN) :: analysis_interval_sfc, end_fdda_hour_sfc
95 INTEGER, INTENT(IN) :: grid_sfdda
97 INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
98 if_no_pbl_nudging_q
99 INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_q
100 INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_q
101 INTEGER, INTENT(IN) :: if_ramping
103 INTEGER , INTENT(IN) :: id
104 REAL, INTENT(IN) :: DT, dx, xtime, dtramp_min
106 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
107 ims,ime, jms,jme, kms,kme, &
108 its,ite, jts,jte, kts,kte
110 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
111 INTENT(IN) :: qv3d, &
112 p3d, &
113 pi3d, &
114 th3d, &
115 t3d, &
116 z, &
117 z_at_w
119 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
120 INTENT(INOUT) :: rundgdten, &
121 rvndgdten, &
122 rthndgdten, &
123 rqvndgdten
124 !
125 ! FASDAS
126 !
127 INTEGER, INTENT(IN) :: fasdas
128 REAL, DIMENSION( ims:ime, jms:jme ), &
129 INTENT(INOUT) :: SDA_HFX, &
130 SDA_QFX
131 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
132 INTENT(INOUT) :: HFX_FDDA
133 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
134 INTENT(IN) :: dz8w
135 !
136 ! END FASDAS
137 !
139 REAL, DIMENSION( ims:ime, jms:jme ), &
140 INTENT(INOUT) :: rmundgdten
142 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
143 INTENT(IN) :: u_ndg_old, &
144 v_ndg_old, &
145 t_ndg_old, &
146 q_ndg_old, &
147 u_ndg_new, &
148 v_ndg_new, &
149 t_ndg_new, &
150 q_ndg_new
152 REAL, DIMENSION( ims:ime, jms:jme ), &
153 INTENT(IN) :: u10_ndg_old, &
154 v10_ndg_old, &
155 t2_ndg_old, &
156 th2_ndg_old, &
157 q2_ndg_old, &
158 rh_ndg_old, &
159 psl_ndg_old, &
160 ps_ndg_old, &
161 u10_ndg_new, &
162 v10_ndg_new, &
163 t2_ndg_new, &
164 th2_ndg_new, &
165 q2_ndg_new, &
166 rh_ndg_new, &
167 psl_ndg_new, &
168 ps_ndg_new
170 REAL, DIMENSION( ims:ime, jms:jme ), &
171 INTENT(IN) :: tob_ndg_old, &
172 tob_ndg_new
174 REAL, DIMENSION( ims:ime, jms:jme ), &
175 INTENT(INOUT) :: mu_ndg_old, &
176 mu_ndg_new
178 REAL, DIMENSION( ims:ime, jms:jme ), &
179 INTENT(IN) :: odis_ndg_old, odis_ndg_new
181 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
182 INTENT(IN) :: u3d, &
183 v3d
185 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
186 ht, &
187 regime, &
188 znt
190 REAL, INTENT(IN) :: guv, gt, gq
191 REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc, rinblw
193 INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, k0, j0
194 REAL :: xtime_old, xtime_new, coef, val_analysis
195 INTEGER :: kpbl, dbg_level
197 REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac, actual_end_fdda_min
198 !BPR BEGIN
199 REAL :: tfac_sfc, actual_end_fdda_min_sfc
200 !BPR END
201 REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ) :: wpbl ! 1: u, 2: v, 3: t, 4: q
202 REAL, DIMENSION( kts:kte, 4 ) :: wzfac ! 1: u, 2: v, 3: t, 4: q
204 ! TWG 2015 Pseudo Radiative Flux
205 REAL, DIMENSION( its:ite, kts:kte, jts:jte) :: rho_air
206 ! END TWG 2015
208 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
210 CHARACTER (LEN=256) :: message
211 INTEGER :: int4
214 int4 = 1 ! 1: temporal interpolation. else: target nudging toward *_ndg_new values
216 actual_end_fdda_min = end_fdda_hour*60.0
217 !BPR BEGIN
218 actual_end_fdda_min_sfc = end_fdda_hour*60.0
219 !BPR END
220 IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
221 actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
222 !BPR BEGIN
223 actual_end_fdda_min_sfc = end_fdda_hour_sfc*60.0 + ABS(dtramp_min)
224 !BPR END
225 !BPR BEGIN
226 !IF( xtime > actual_end_fdda_min ) THEN
227 IF( ( xtime > actual_end_fdda_min ) .AND. ( xtime > actual_end_fdda_min_sfc ) ) THEN
228 !BPR END
229 ! If xtime is greater than the end time, no need to calculate tendencies. Just set the tendencies
230 ! to zero to turn off nudging and return.
231 DO j = jts, jte
232 DO k = kts, kte
233 DO i = its, ite
234 RUNDGDTEN(i,k,j) = 0.0
235 RVNDGDTEN(i,k,j) = 0.0
236 RTHNDGDTEN(i,k,j) = 0.0
237 RQVNDGDTEN(i,k,j) = 0.0
238 IF( k .EQ. kts ) RMUNDGDTEN(i,j) = 0.
239 ENDDO
240 ENDDO
241 ENDDO
242 RETURN
243 ENDIF
245 IF( analysis_interval <= 0 )CALL wrf_error_fatal('In grid FDDA, gfdda_interval_m must be > 0')
246 xtime_old = FLOOR(xtime/analysis_interval) * analysis_interval * 1.0
247 xtime_new = xtime_old + analysis_interval
248 IF( int4 == 1 ) THEN
249 coef = (xtime-xtime_old)/(xtime_new-xtime_old)
250 ELSE
251 coef = 1.0 ! Nudging toward a target value (*_ndg_new values)
252 ENDIF
254 IF ( wrf_dm_on_monitor()) THEN
256 CALL get_wrf_debug_level( dbg_level )
258 IF( xtime-xtime_old < 0.5*dt/60.0 ) THEN
260 IF( xtime < end_fdda_hour*60.0 ) THEN
261 WRITE(message,'(a,i1,a,f10.3,a)') &
262 'D0',id,' 3-D analysis nudging reads new data at time = ', xtime, ' min.'
263 CALL wrf_message( TRIM(message) )
264 WRITE(message,'(a,i1,a,2f8.2,a)') &
265 'D0',id,' 3-D analysis nudging bracketing times = ', xtime_old, xtime_new, ' min.'
266 CALL wrf_message( TRIM(message) )
267 ENDIF
269 actual_end_fdda_min = end_fdda_hour*60.0
270 IF( if_ramping == 1 .AND. dtramp_min > 0.0 ) &
271 actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
273 IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min ) THEN
274 ! Find the mid point of the tile and print out the sample values
275 i0 = (ite-its)/2+its
276 j0 = (jte-jts)/2+jts
278 IF( guv > 0.0 ) THEN
279 DO k = kts, kte
280 WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
281 ' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
282 ' u_ndg_old=', u_ndg_old(i0,k,j0), ' u_ndg_new=', u_ndg_new(i0,k,j0)
283 CALL wrf_message( TRIM(message) )
284 ENDDO
285 WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
286 ' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
287 ' mu_ndg_old=', mu_ndg_old(i0,j0), ' mu_ndg_new=', mu_ndg_new(i0,j0)
288 CALL wrf_message( TRIM(message) )
289 DO k = kts, kte
290 WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
291 ' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
292 ' v_ndg_old=', v_ndg_old(i0,k,j0), ' v_ndg_new=', v_ndg_new(i0,k,j0)
293 CALL wrf_message( TRIM(message) )
294 ENDDO
295 ENDIF
297 IF( gt > 0.0 ) THEN
298 DO k = kts, kte
299 WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
300 ' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
301 ' t_ndg_old=', t_ndg_old(i0,k,j0), ' t_ndg_new=', t_ndg_new(i0,k,j0)
302 CALL wrf_message( TRIM(message) )
303 ENDDO
304 ENDIF
306 IF( gq > 0.0 ) THEN
307 DO k = kts, kte
308 WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
309 ' D0',id,' sample 3-D analysis values at i,k,j=', i0, k, j0, &
310 ' q_ndg_old=', q_ndg_old(i0,k,j0), ' q_ndg_new=', q_ndg_new(i0,k,j0)
311 CALL wrf_message( TRIM(message) )
312 ENDDO
313 ENDIF
315 IF( int4 == 1 ) then
316 WRITE(message,'(a,i1,a)') ' D0',id, &
317 ' 3-D nudging towards the temporally interpolated analysis'
318 ELSE
319 WRITE(message,'(a,i1,a)') ' D0',id, &
320 ' 3-D nudging towards the target analysis'
321 ENDIF
323 ENDIF
324 ENDIF
325 ENDIF
327 jtsv=MAX0(jts,jds+1)
328 itsu=MAX0(its,ids+1)
330 jtf=MIN0(jte,jde-1)
331 ktf=MIN0(kte,kde-1)
332 itf=MIN0(ite,ide-1)
333 !
334 ! If the user-defined namelist switches (if_no_pbl_nudging_uv, if_no_pbl_nudging_t,
335 ! if_no_pbl_nudging_q switches) are set to 1, compute the weighting function, wpbl(:,k,:,:),
336 ! based on the PBL depth. wpbl = 1 above the PBL and wpbl = 0 in the PBL. If all
337 ! the switches are set to zero, wpbl = 1 (default value).
338 !
339 wpbl(:,:,:,:) = 1.0
341 IF( if_no_pbl_nudging_uv == 1 .OR. grid_sfdda >= 1 ) THEN
343 DO j=jts,jtf
344 DO i=itsu,itf
346 kpbl = 1
347 zpbl = 0.5 * ( pblh(i-1,j) + pblh(i,j) )
349 loop_ku: DO k=kts,ktf
350 ! zagl = 0.5 * ( z(i-1,k,j)-ht(i-1,j) + z(i,k,j)-ht(i,j) )
351 zagl_bot = 0.5 * ( z_at_w(i-1,k, j)-ht(i-1,j) + z_at_w(i,k, j)-ht(i,j) )
352 zagl_top = 0.5 * ( z_at_w(i-1,k+1,j)-ht(i-1,j) + z_at_w(i,k+1,j)-ht(i,j) )
353 IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
354 kpbl = k
355 EXIT loop_ku
356 ENDIF
357 ENDDO loop_ku
359 DO k=kts,ktf
360 IF( k <= kpbl ) wpbl(i, k, j, 1) = 0.0
361 IF( k == kpbl+1 ) wpbl(i, k, j, 1) = 0.1
362 IF( k > kpbl+1 ) wpbl(i, k, j, 1) = 1.0
363 ENDDO
365 ENDDO
366 ENDDO
368 DO i=its,itf
369 DO j=jtsv,jtf
371 kpbl = 1
372 zpbl = 0.5 * ( pblh(i,j-1) + pblh(i,j) )
374 loop_kv: DO k=kts,ktf
375 ! zagl = 0.5 * ( z(i,k,j-1)-ht(i,j-1) + z(i,k,j)-ht(i,j) )
376 zagl_bot = 0.5 * ( z_at_w(i,k, j-1)-ht(i,j-1) + z_at_w(i,k, j)-ht(i,j) )
377 zagl_top = 0.5 * ( z_at_w(i,k+1,j-1)-ht(i,j-1) + z_at_w(i,k+1,j)-ht(i,j) )
378 IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
379 kpbl = k
380 EXIT loop_kv
381 ENDIF
382 ENDDO loop_kv
384 DO k=kts,ktf
385 IF( k <= kpbl ) wpbl(i, k, j, 2) = 0.0
386 IF( k == kpbl+1 ) wpbl(i, k, j, 2) = 0.1
387 IF( k > kpbl+1 ) wpbl(i, k, j, 2) = 1.0
388 ENDDO
390 ENDDO
391 ENDDO
393 ENDIF
395 IF( if_no_pbl_nudging_t == 1 .OR. grid_sfdda >= 1 ) THEN
397 DO j=jts,jtf
398 DO i=its,itf
400 kpbl = 1
401 zpbl = pblh(i,j)
403 loop_kt: DO k=kts,ktf
404 ! zagl = z(i,k,j)-ht(i,j)
405 zagl_bot = z_at_w(i,k, j)-ht(i,j)
406 zagl_top = z_at_w(i,k+1,j)-ht(i,j)
407 IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
408 kpbl = k
409 EXIT loop_kt
410 ENDIF
411 ENDDO loop_kt
413 DO k=kts,ktf
414 IF( k <= kpbl ) wpbl(i, k, j, 3) = 0.0
415 IF( k == kpbl+1 ) wpbl(i, k, j, 3) = 0.1
416 IF( k > kpbl+1 ) wpbl(i, k, j, 3) = 1.0
417 ENDDO
419 ENDDO
420 ENDDO
422 ENDIF
424 IF( if_no_pbl_nudging_q == 1 .OR. grid_sfdda >= 1 ) THEN
426 DO j=jts,jtf
427 DO i=its,itf
429 kpbl = 1
430 zpbl = pblh(i,j)
432 loop_kq: DO k=kts,ktf
433 ! zagl = z(i,k,j)-ht(i,j)
434 zagl_bot = z_at_w(i,k, j)-ht(i,j)
435 zagl_top = z_at_w(i,k+1,j)-ht(i,j)
436 IF( zpbl >= zagl_bot .AND. zpbl < zagl_top ) THEN
437 kpbl = k
438 EXIT loop_kq
439 ENDIF
440 ENDDO loop_kq
442 DO k=kts,ktf
443 IF( k <= kpbl ) wpbl(i, k, j, 4) = 0.0
444 IF( k == kpbl+1 ) wpbl(i, k, j, 4) = 0.1
445 IF( k > kpbl+1 ) wpbl(i, k, j, 4) = 1.0
446 ENDDO
448 ENDDO
449 ENDDO
451 ENDIF
452 !
453 ! If the user-defined namelist switches (if_zfac_uv, if_zfac_t,
454 ! if_zfac_q swithes) are set to 1, compute the weighting function, wzfac(k,:),
455 ! based on the namelist specified k values (k_zfac_uv, k_zfac_t and k_zfac_q) below which analysis
456 ! nudging is turned off (wzfac = 1 above k_zfac_x and = 0 in below k_zfac_x). If all
457 ! the switche are set to zero, wzfac = 1 (default value).
458 !
459 wzfac(:,:) = 1.0
461 IF( if_zfac_uv == 1 ) THEN
463 DO j=jts,jtf
464 DO i=itsu,itf
465 DO k=kts,ktf
466 IF( k <= k_zfac_uv ) wzfac(k, 1:2) = 0.0
467 IF( k == k_zfac_uv+1 ) wzfac(k, 1:2) = 0.1
468 IF( k > k_zfac_uv+1 ) wzfac(k, 1:2) = 1.0
469 ENDDO
470 ENDDO
471 ENDDO
473 ENDIF
475 IF( if_zfac_t == 1 ) THEN
477 DO j=jts,jtf
478 DO i=itsu,itf
479 DO k=kts,ktf
480 IF( k <= k_zfac_t ) wzfac(k, 3) = 0.0
481 IF( k == k_zfac_t+1 ) wzfac(k, 3) = 0.1
482 IF( k > k_zfac_t+1 ) wzfac(k, 3) = 1.0
483 ENDDO
484 ENDDO
485 ENDDO
487 ENDIF
489 IF( if_zfac_q == 1 ) THEN
491 DO j=jts,jtf
492 DO i=itsu,itf
493 DO k=kts,ktf
494 IF( k <= k_zfac_q ) wzfac(k, 4) = 0.0
495 IF( k == k_zfac_q+1 ) wzfac(k, 4) = 0.1
496 IF( k > k_zfac_q+1 ) wzfac(k, 4) = 1.0
497 ENDDO
498 ENDDO
499 ENDDO
501 ENDIF
502 !
503 ! If if_ramping and dtramp_min are defined by user, compute a time weighting function, tfac,
504 ! for analysis nudging so that at the end of the nudging period (which has to be at a
505 ! analysis time) we ramp down the nudging coefficient, based on the use-defined sign of dtramp_min.
506 !
507 ! When dtramp_min is negative, ramping ends at end_fdda_hour and starts at
508 ! end_fdda_hour-ABS(dtramp_min).
509 !
510 ! When dtramp_min is positive, ramping starts at end_fdda_hour and ends at
511 ! end_fdda_hour+ABS(dtramp_min).
512 ! BPR BEGIN
513 ! In this case, during the rampdown we nudge towards the most recent past analysis and ignore
514 ! the future analysis (implemented by setting coef = 0.0)
515 ! THE FOLLOWING COMMENT IS NO LONGER APPLICABLE
516 ! In this case, the obs values are extrapolated using
517 ! the obs tendency saved from the previous FDDA window. More specifically for extrapolation,
518 ! coef (see codes below) is recalculated to reflect extrapolation during the ramping period.
519 ! THE PRECEDING COMMENT IS NOT LONGER APPLICABLE
520 ! BPR END
521 !
522 tfac = 1.0
523 !BPR BEGIN
524 tfac_sfc = 1.0
525 !BPR END
527 IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
529 IF( dtramp_min <= 0.0 ) THEN
530 actual_end_fdda_min = end_fdda_hour*60.0
531 ELSE
532 actual_end_fdda_min = end_fdda_hour*60.0 + dtramp_min
533 ENDIF
535 IF( xtime < actual_end_fdda_min-ABS(dtramp_min) )THEN
536 tfac = 1.0
537 ELSEIF( xtime >= actual_end_fdda_min-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min )THEN
538 tfac = ( actual_end_fdda_min - xtime ) / ABS(dtramp_min)
539 !BPR BEGIN
540 !The original method assumed that to be here in the code with dtramp_min>0
541 !meant that we were after the valid time of *_ndg_new and so used the
542 !values *_ndg_old and *_ndg_new to extrapolate a current analysis value.
543 !HOWEVER, in practice once we get to the valid time of *_ndg_new WRF will
544 !read in the next pair of *_ndg_old and *_ndg_new values, even if we are
545 !at the beginning of the rampdown period. Since the *_ndg_new values have
546 !a valid time after the beginning of the rampdown period we do not want
547 !to use them as we would be using analyses valid after the supposed end
548 !time of the analysis nudging.
549 !IF( dtramp_min > 0.0 ) coef = (xtime-xtime_old+analysis_interval)/(analysis_interval*1.0)
550 !Use only the *_ndg_old values
551 IF( dtramp_min > 0.0 ) coef = 0.0
552 !BPR END
553 ELSE
554 tfac = 0.0
555 ENDIF
557 !BPR BEGIN
558 !Now calculate the same quantities for surface analysis nudging
559 !Add actual_end_fdda_min_sfc, tfac_sfc
560 IF( dtramp_min <= 0.0 ) THEN
561 actual_end_fdda_min_sfc = end_fdda_hour_sfc*60.0
562 ELSE
563 actual_end_fdda_min_sfc = end_fdda_hour_sfc*60.0 + dtramp_min
564 ENDIF
566 IF( xtime < actual_end_fdda_min_sfc-ABS(dtramp_min) )THEN
567 tfac_sfc = 1.0
568 ELSEIF( xtime >= actual_end_fdda_min_sfc-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min_sfc )THEN
569 tfac_sfc = ( actual_end_fdda_min_sfc - xtime ) / ABS(dtramp_min)
570 ELSE
571 tfac_sfc = 0.0
572 ENDIF
573 !BPR END
575 ENDIF
577 !
578 ! Surface Analysis Nudging
579 !
580 IF( grid_sfdda >= 1 ) THEN
581 CALL SFDDAGD(itimestep,dx,dt,xtime, id, &
582 analysis_interval_sfc, end_fdda_hour_sfc, guv_sfc, gt_sfc, gq_sfc, &
583 rinblw, &
584 u3d,v3d,th3d,t3d, &
585 qv3d, &
586 p3d,pi3d, &
587 u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
588 rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
589 u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
590 rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
591 RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN,&
592 pblh, ht, regime, znt, z, z_at_w, &
593 ids,ide, jds,jde, kds,kde, &
594 ims,ime, jms,jme, kms,kme, &
595 its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
596 !BPR BEGIN
597 ! actual_end_fdda_min, tfac &
598 actual_end_fdda_min_sfc, tfac_sfc &
599 !BPR END
600 !
601 ! FASDAS
602 !
603 ,fasdas, SDA_HFX, SDA_QFX &
604 !
605 ! END FASDAS
606 !
607 )
608 ENDIF
609 !
610 ! Compute 3-D nudging tendencies for u, v, t and q
611 !
612 DO j=jts,jtf
613 DO k=kts,ktf
614 DO i=itsu,itf
615 val_analysis = u_ndg_old(i,k,j) *( 1.0 - coef ) + u_ndg_new(i,k,j) * coef
616 RUNDGDTEN(i,k,j) = RUNDGDTEN(i,k,j) + guv * wpbl(i,k,j,1) * wzfac(k,1) * tfac * &
617 ( val_analysis - u3d(i,k,j) )
618 ENDDO
619 ENDDO
620 ENDDO
622 DO j=jtsv,jtf
623 DO k=kts,ktf
624 DO i=its,itf
625 val_analysis = v_ndg_old(i,k,j) *( 1.0 - coef ) + v_ndg_new(i,k,j) * coef
626 RVNDGDTEN(i,k,j) = RVNDGDTEN(i,k,j) + guv * wpbl(i,k,j,2) * wzfac(k,2) * tfac * &
627 ( val_analysis - v3d(i,k,j) )
628 ENDDO
629 ENDDO
630 ENDDO
632 DO j=jts,jtf
633 DO k=kts,ktf
634 DO i=its,itf
637 val_analysis = t_ndg_old(i,k,j) *( 1.0 - coef ) + t_ndg_new(i,k,j) * coef
638 RTHNDGDTEN(i,k,j) = RTHNDGDTEN(i,k,j) + gt * wpbl(i,k,j,3) * wzfac(k,3) * tfac * &
639 ( val_analysis - th3d(i,k,j) + 300.0 )
640 !
641 !FASDAS
642 !
643 !TWG 2015 Pseudo Radiative Flux
645 rho_air(i,k,j) = p3d(i,k,j)/(287.0*th3d(i,k,j))
646 HFX_FDDA(i,k,j) = rho_air(i,k,j)*1004.0*RTHNDGDTEN(i,k,j)*dz8w(i,k,j)
648 !TWG 2015 END
649 !
650 !END FASDAS
651 !
652 val_analysis = q_ndg_old(i,k,j) *( 1.0 - coef ) + q_ndg_new(i,k,j) * coef
653 RQVNDGDTEN(i,k,j) = RQVNDGDTEN(i,k,j) + gq * wpbl(i,k,j,4) * wzfac(k,4) * tfac * &
654 ( val_analysis - qv3d(i,k,j) )
656 ENDDO
657 ENDDO
658 ENDDO
660 !BPR BEGIN
661 !Diagnostic print
662 IF ( wrf_dm_on_monitor()) THEN
663 IF( dbg_level .GE. 10 ) THEN
664 i0 = (ite-its)/2+its
665 j0 = (jte-jts)/2+jts
666 k0 = (kte-kts)/2+kts
667 WRITE(message,'(a,i1,a,f7.2,a,i4,a,i4,a,i4,a,f7.2,a,f4.2,a,f7.2,a,f4.2,a,f4.2 )') &
668 ' D0',id,' At xtime=',xtime,' FDDA sample pot. temp. analysis (i=',i0,', j=',j0,' k=', &
669 k0,') = (t_ndg_old=', t_ndg_old(i0,k0,j0),') * ',1.0-coef,' + (t_ndg_new=', &
670 t_ndg_new(i0,k0,j0),') * ',coef, ' where tfac=',tfac
671 CALL wrf_message( TRIM(message) )
672 WRITE(message,'(a,i1,a,f7.2,a,i4,a,f7.2,a,f7.2)') &
673 ' D0',id,' xtime_old=',xtime_old,' analysis_interval=',analysis_interval,' actual_end_fdda_min=', &
674 actual_end_fdda_min,' dtramp_min=',dtramp_min
675 CALL wrf_message( TRIM(message) )
676 END IF
677 END IF
678 !BPR END
680 END SUBROUTINE fddagd
682 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
683 SUBROUTINE sfddagd(itimestep,dx,dt,xtime, &
684 id, analysis_interval_sfc, end_fdda_hour_sfc, &
685 guv_sfc, gt_sfc, gq_sfc, rinblw, &
686 u3d,v3d,th3d,t3d, &
687 qv3d, &
688 p3d,pi3d, &
689 u10_ndg_old, v10_ndg_old, t2_ndg_old, th2_ndg_old, q2_ndg_old, &
690 rh_ndg_old, psl_ndg_old, ps_ndg_old, tob_ndg_old, odis_ndg_old, &
691 u10_ndg_new, v10_ndg_new, t2_ndg_new, th2_ndg_new, q2_ndg_new, &
692 rh_ndg_new, psl_ndg_new, ps_ndg_new, tob_ndg_new, odis_ndg_new, &
693 RUNDGDTEN,RVNDGDTEN,RTHNDGDTEN,RQVNDGDTEN,RMUNDGDTEN, &
694 pblh, ht, regime, znt, z, z_at_w, &
695 ids,ide, jds,jde, kds,kde, &
696 ims,ime, jms,jme, kms,kme, &
697 its,ite, jts,jte, kts,kte, wpbl, wzfac, if_ramping, dtramp_min, &
698 !BPR BEGIN
699 ! actual_end_fdda_min, tfac &
700 actual_end_fdda_min_sfc, tfac_sfc &
701 !BPR END
702 !
703 ! FASDAS
704 !
705 ,fasdas, SDA_HFX, SDA_QFX &
706 !
707 ! END FASDAS
708 !
709 )
710 !-------------------------------------------------------------------
711 USE module_model_constants
713 implicit none
715 !-------------------------------------------------------------------
716 !
717 ! This code was implemented by Aijun Deng (Penn State). The 3-D analysis nudging was comleted
718 ! and released in December 2006. The surface analysis nudging capability was added and
719 ! released in March 2009 with WRFV3.1.
720 !
721 !-- u3d 3d u-velocity staggered on u points
722 !-- v3d 3d v-velocity staggered on v points
723 !-- th3d 3d potential temperature (k)
724 !-- t3d temperature (k)
725 !-- qv3d 3d water vapor mixing ratio (kg/kg)
726 !-- p3d 3d pressure (pa)
727 !-- pi3d 3d exner function (dimensionless)
728 !-- rundgdten staggered u tendency due to
729 ! fdda grid nudging (m/s/s)
730 !-- rvndgdten staggered v tendency due to
731 ! fdda grid nudging (m/s/s)
732 !-- rthndgdten theta tendency due to
733 ! fdda grid nudging (K/s)
734 !-- rqvndgdten qv tendency due to
735 ! fdda grid nudging (kg/kg/s)
736 !-- rmundgdten mu tendency due to
737 ! fdda grid nudging (Pa/s)
738 !-- ids start index for i in domain
739 !-- ide end index for i in domain
740 !-- jds start index for j in domain
741 !-- jde end index for j in domain
742 !-- kds start index for k in domain
743 !-- kde end index for k in domain
744 !-- ims start index for i in memory
745 !-- ime end index for i in memory
746 !-- jms start index for j in memory
747 !-- jme end index for j in memory
748 !-- kms start index for k in memory
749 !-- kme end index for k in memory
750 !-- its start index for i in tile
751 !-- ite end index for i in tile
752 !-- jts start index for j in tile
753 !-- jte end index for j in tile
754 !-- kts start index for k in tile
755 !-- kte end index for k in tile
756 !-------------------------------------------------------------------
757 !
758 INTEGER, INTENT(IN) :: itimestep, analysis_interval_sfc, end_fdda_hour_sfc
760 INTEGER , INTENT(IN) :: id
761 REAL, INTENT(IN) :: dx,DT, xtime
763 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
764 ims,ime, jms,jme, kms,kme, &
765 its,ite, jts,jte, kts,kte
767 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
768 INTENT(IN) :: qv3d, &
769 p3d, &
770 pi3d, &
771 th3d, &
772 t3d, &
773 z, &
774 z_at_w
776 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
777 INTENT(INOUT) :: rundgdten, &
778 rvndgdten, &
779 rthndgdten, &
780 rqvndgdten
782 REAL, DIMENSION( ims:ime, jms:jme ), &
783 INTENT(INOUT) :: rmundgdten
785 REAL, DIMENSION( ims:ime, jms:jme ), &
786 INTENT(IN) :: u10_ndg_old, &
787 v10_ndg_old, &
788 t2_ndg_old, &
789 th2_ndg_old, &
790 q2_ndg_old, &
791 rh_ndg_old, &
792 psl_ndg_old, &
793 ps_ndg_old, &
794 u10_ndg_new, &
795 v10_ndg_new, &
796 t2_ndg_new, &
797 th2_ndg_new, &
798 q2_ndg_new, &
799 rh_ndg_new, &
800 psl_ndg_new, &
801 ps_ndg_new
803 REAL, DIMENSION( ims:ime, jms:jme ), &
804 INTENT(IN) :: tob_ndg_old, &
805 tob_ndg_new
807 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
808 INTENT(IN) :: u3d, &
809 v3d
811 REAL, DIMENSION( ims:ime, jms:jme ), &
812 INTENT(IN) :: odis_ndg_old, odis_ndg_new
814 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: pblh, &
815 ht, &
816 regime, &
817 znt
819 REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc, rinblw
821 INTEGER :: i, j, k, itsu, jtsv, itf, jtf, ktf, i0, j0
822 REAL :: xtime_old_sfc, xtime_new_sfc, coef, val_analysis, es
823 INTEGER :: kpbl, dbg_level
825 REAL :: zpbl, zagl, zagl_bot, zagl_top, tfac_sfc, actual_end_fdda_min_sfc
827 REAL, DIMENSION( its:ite, kts:kte, jts:jte, 4 ), &
828 INTENT(INout) :: wpbl ! 1: u, 2: v, 3: t, 4: q
829 REAL, DIMENSION( kts:kte, 4 ), &
830 INTENT(IN) :: wzfac ! 1: u, 2: v, 3: t, 4: q
831 REAL, DIMENSION( its:ite, jts:jte) :: wndcor_u, wndcor_v
832 REAL, DIMENSION( its-2:ite+2, jts-2:jte+2) :: blw_old, blw_new
833 REAL, DIMENSION( its:ite, kts:kte, jts:jte) :: qsat
835 REAL :: m, b=1.8, blw, rindx, x
837 REAL :: difz, wr14, wr1z, wr24, wr2z, wndfac, reg, znt0
838 INTEGER, INTENT(IN) :: if_ramping
839 REAL, INTENT(IN) :: dtramp_min
841 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
843 CHARACTER (LEN=256) :: message
844 INTEGER :: iwinds, idd, iqsat, int4
845 !
846 ! FASDAS
847 !
848 INTEGER, INTENT(IN) :: fasdas
849 REAL, DIMENSION( ims:ime, jms:jme ), &
850 INTENT( OUT) :: SDA_HFX, &
851 SDA_QFX
852 REAL :: stabFac, exf, DiffTN
853 REAL, DIMENSION( its:ite, kts:kte, jts:jte ) :: wpblfasdas
854 !
855 ! END FASDAS
856 !
857 iwinds = 1 ! 1: Scale the surface wind analysis to the lowest model level,
858 ! if the first model half-layer is greater than 10 meters
859 ! and we are in the free convection regime (REGIME=4.0). else: no
860 idd = 1 ! 1: Obs data density correction is applied. else: no
861 iqsat = 1 ! 1: Remove super saturation. eles: no
862 int4 = 1 ! 1: temporal ionterpolation. else: target nudging toward *_ndg_new values
864 IF( analysis_interval_sfc <= 0 )CALL wrf_error_fatal('In grid sfc FDDA, sgfdda_interval_m must be > 0')
865 xtime_old_sfc = FLOOR(xtime/analysis_interval_sfc) * analysis_interval_sfc * 1.0
866 xtime_new_sfc = xtime_old_sfc + analysis_interval_sfc
867 IF( int4 == 1 ) THEN
868 coef = (xtime-xtime_old_sfc)/(xtime_new_sfc-xtime_old_sfc) ! Temporal interpolation
869 ELSE
870 coef = 1.0 ! Nudging toward a target value (*_ndg_new values)
871 ENDIF
873 IF ( wrf_dm_on_monitor()) THEN
875 CALL get_wrf_debug_level( dbg_level )
877 IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
879 IF( xtime < end_fdda_hour_sfc*60.0 ) THEN
880 WRITE(message,'(a,i1,a,f10.3,a)') &
881 'D0',id,' surface analysis nudging reads new data at time = ', xtime, ' min.'
882 CALL wrf_message( TRIM(message) )
883 WRITE(message,'(a,i1,a,2f8.2,a)') &
884 'D0',id,' surface analysis nudging bracketing times = ', xtime_old_sfc, xtime_new_sfc, ' min.'
885 CALL wrf_message( TRIM(message) )
886 ENDIF
888 IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
889 ! Find the mid point of the tile and print out the sample values
890 i0 = (ite-its)/2+its
891 j0 = (jte-jts)/2+jts
893 IF( guv_sfc > 0.0 ) THEN
894 WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
895 ' D0',id,' sample surface analysis values at i,j=', i0, j0, &
896 ' u10_ndg_old=', u10_ndg_old(i0,j0), ' u10_ndg_new=', u10_ndg_new(i0,j0)
897 CALL wrf_message( TRIM(message) )
898 WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
899 ' D0',id,' sample surface analysis values at i,j=', i0, j0, &
900 ' v10_ndg_old=', v10_ndg_old(i0,j0), ' v10_ndg_new=', v10_ndg_new(i0,j0)
901 CALL wrf_message( TRIM(message) )
902 ENDIF
904 IF( gt_sfc > 0.0 ) THEN
905 WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
906 ' D0',id,' sample surface analysis values at i,j=', i0, j0, &
907 ' th2_ndg_old=', th2_ndg_old(i0,j0), ' th2_ndg_new=', th2_ndg_new(i0,j0)
908 CALL wrf_message( TRIM(message) )
909 ENDIF
911 IF( gq_sfc > 0.0 ) THEN
912 WRITE(message,'(a,i1,a,2i4,a,f10.4,a,f10.4)') &
913 ' D0',id,' sample surface analysis values at i,j=', i0, j0, &
914 ' q2_ndg_old=', q2_ndg_old(i0,j0), ' q2_ndg_new=', q2_ndg_new(i0,j0)
915 CALL wrf_message( TRIM(message) )
916 ENDIF
918 IF( iwinds == 1 ) &
919 WRITE(message,'(a,i1,a)') ' D0',id, &
920 ' surface wind analysis s scaled to the lowest model level, if dz1 > 10m and REGIME=4.'
922 IF( idd == 1 ) &
923 WRITE(message,'(a,i1,a)') ' D0',id, &
924 ' obs data density is used for additional weighting function'
926 IF( iqsat == 1 ) &
927 WRITE(message,'(a,i1,a)') ' D0',id, &
928 ' super saturation is not allowed for q analysis'
930 IF( int4 == 1 ) then
931 WRITE(message,'(a,i1,a)') ' D0',id, &
932 ' surface nudging towards the temporally interpolated analysis'
933 ELSE
934 WRITE(message,'(a,i1,a)') ' D0',id, &
935 ' surface nudging towards the target analysis'
936 ENDIF
938 ENDIF
939 ENDIF
940 ENDIF
943 jtsv=MAX0(jts,jds+1)
944 itsu=MAX0(its,ids+1)
946 jtf=MIN0(jte,jde-1)
947 ktf=MIN0(kte,kde-1)
948 itf=MIN0(ite,ide-1)
950 !
951 ! Compute the vertical weighting function to scale the surface wind analysis to
952 ! the lowest model level, if the first model half-layer is greater
953 ! than 10 meters and we are in the free convection regime (REGIME=4.0).
954 !
955 IF( iwinds == 1 ) THEN
956 wndcor_u(:,:) = 1.0
957 DO j=jts,jtf
958 DO i=itsu,itf
959 reg = 0.5 * ( regime(i-1, j) + regime(i, j) )
960 difz = 0.5 * ( z(i-1,1,j) - ht(i-1,j) &
961 + z(i, 1,j) - ht(i, j) )
962 IF( reg > 3.5 .AND. difz > 10.0 ) THEN
963 znt0 = 0.5 * ( znt(i-1, j) + znt(i, j) )
964 IF( znt0 <= 0.2) THEN
965 wndcor_u(i,j) = 1.0+0.320*znt0**0.2
966 ELSE
967 wndcor_u(i,j) = 1.169+0.315*znt0
968 ENDIF
970 wr14 = log(40.0/0.05)
971 wr1z = log(difz/0.05)
972 wr24 = log(40.0/1.0)
973 wr2z = log(difz/1.0)
974 wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
975 wndcor_u(i,j) = wndfac*wndcor_u(i,j)
976 ENDIF
977 ENDDO
978 ENDDO
980 IF ( wrf_dm_on_monitor()) THEN
981 IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
982 IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
983 i0 = (ite-its)/2+its
984 j0 = (jte-jts)/2+jts
985 WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
986 ' D0',id,' sample wndcor_u values at i,j=', i0, j0, &
987 ' wndcor_u=', wndcor_u(i0,j0)
988 CALL wrf_message( TRIM(message) )
989 ENDIF
990 ENDIF
991 ENDIF
992 ELSE
993 wndcor_u(:,:) = 1.0
994 ENDIF
996 IF( iwinds == 1 ) THEN
997 wndcor_v(:,:) = 1.0
998 DO j=jtsv,jtf
999 DO i=its,itf
1000 reg = 0.5 * ( regime(i, j-1) + regime(i, j) )
1001 difz = 0.5 * ( z(i,1,j-1) - ht(i,j-1) &
1002 + z(i,1,j ) - ht(i,j ) )
1003 IF( reg > 3.5 .AND. difz > 10.0 ) THEN
1004 znt0 = 0.5 * ( znt(i, j-1) + znt(i, j) )
1005 IF( znt0 <= 0.2) THEN
1006 wndcor_v(i,j) = 1.0+0.320*znt0**0.2
1007 ELSE
1008 wndcor_v(i,j) = 1.169+0.315*znt0
1009 ENDIF
1011 wr14 = log(40.0/0.05)
1012 wr1z = log(difz/0.05)
1013 wr24 = log(40.0/1.0)
1014 wr2z = log(difz/1.0)
1015 wndfac = 0.5*(WR1Z/WR14+WR2Z/WR24)
1016 wndcor_v(i,j) = wndfac*wndcor_v(i,j)
1017 ENDIF
1018 ENDDO
1019 ENDDO
1020 IF ( wrf_dm_on_monitor()) THEN
1021 IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
1022 IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
1023 i0 = (ite-its)/2+its
1024 j0 = (jte-jts)/2+jts
1025 WRITE(message,'(a,i1,a,2i4,a,f10.4)') &
1026 ' D0',id,' sample wndcor_v values at i,j=', i0, j0, &
1027 ' wndcor_v=', wndcor_v(i0,j0)
1028 CALL wrf_message( TRIM(message) )
1029 ENDIF
1030 ENDIF
1031 ENDIF
1032 ELSE
1033 wndcor_v(:,:) = 1.0
1034 ENDIF
1035 !
1036 ! Compute saturation mixing ratio so that nudging to a super-saturated state
1037 ! is not allowed.
1038 !
1039 IF( iqsat == 1 ) THEN
1040 DO j=jts,jtf
1041 DO k=kts,ktf
1042 DO i=its,itf
1043 es = SVP1*EXP(SVP2*(t3d(i,k,j)-SVPT0)/(t3d(i,k,j)-SVP3)) * 10.0 ! mb
1044 qsat(i,k,j) = EP_2*es/(p3d(i,k,j)/100.0-es)
1045 ENDDO
1046 ENDDO
1047 ENDDO
1049 IF ( wrf_dm_on_monitor()) THEN
1050 IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
1051 IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
1052 i0 = (ite-its)/2+its
1053 j0 = (jte-jts)/2+jts
1054 DO k = kts, kte
1055 WRITE(message,'(a,i1,a,3i4,a,f10.4,a,f10.4)') &
1056 ' D0',id,' sample moisture values (g/kg) at i,k,j=', i0, k, j0, &
1057 ' qv3d=', qv3d(i0,k,j0)*1000.0, ' qsat=', qsat(i0,k,j0)*1000.0
1058 CALL wrf_message( TRIM(message) )
1059 ENDDO
1060 ENDIF
1061 ENDIF
1062 ENDIF
1063 ENDIF
1064 !
1065 ! Obs data density weighting.
1066 !
1067 IF( idd == 1 ) THEN
1069 IF( rinblw < 0.001 ) THEN
1070 IF ( wrf_dm_on_monitor()) THEN
1071 WRITE(message,'(a)') 'Error in rinblw, please specify a reasonable value ***'
1072 CALL wrf_message( TRIM(message) )
1073 ENDIF
1074 CALL wrf_error_fatal('In grid FDDA')
1075 ENDIF
1077 rindx = rinblw*1000.0/dx
1078 m = -0.8*2.0/rindx
1080 DO j=MAX(jts-2,jds),MIN(jte+2,jde-1)
1081 DO i=MAX(its-2,ids),MIN(ite+2,ide-1)
1082 IF( odis_ndg_old(i,j) < 0.5*rinblw ) THEN
1083 blw_old(i,j) = 1.0
1084 ELSE
1085 x = min( odis_ndg_old(i,j)*1000./dx, rindx )
1086 blw_old(i,j) = m * x + b
1087 ENDIF
1089 IF( odis_ndg_new(i,j) < 0.5*rinblw ) THEN
1090 blw_new(i,j) = 1.0
1091 ELSE
1092 x = min( odis_ndg_new(i,j)*1000./dx, rindx )
1093 blw_new(i,j) = m * x + b
1094 ENDIF
1095 ENDDO
1096 ENDDO
1098 ! Smoother applies one point outside the tile, but one point in from boundaries
1099 CALL smther(blw_old, its-2,ite+2, jts-2,jte+2, 1, &
1100 MAX(its-1,ids+1), MIN(ite+1,ide-2), MAX(jts-1,jds+1), MIN(jte+1,jde-2))
1101 CALL smther(blw_new, its-2,ite+2, jts-2,jte+2, 1, &
1102 MAX(its-1,ids+1), MIN(ite+1,ide-2), MAX(jts-1,jds+1), MIN(jte+1,jde-2))
1104 WHERE ( blw_old > 1.0)
1105 blw_old = 1.0
1106 END WHERE
1107 WHERE ( blw_new > 1.0)
1108 blw_new = 1.0
1109 END WHERE
1110 WHERE ( blw_old < 0.0)
1111 blw_old = 0.0
1112 END WHERE
1113 WHERE ( blw_new < 0.0)
1114 blw_new = 0.0
1115 END WHERE
1117 IF ( wrf_dm_on_monitor()) THEN
1118 IF( xtime-xtime_old_sfc < 0.5*dt/60.0 ) THEN
1119 IF( dbg_level .GE. 10 .AND. xtime <= actual_end_fdda_min_sfc ) THEN
1120 i0 = (ite-its)/2+its
1121 j0 = (jte-jts)/2+jts
1122 WRITE(message,'(a,i1,a,2i4,4(a,f10.4))') &
1123 ' D0',id,' sample blw values at i,j=', i0, j0, &
1124 ' odis_ndg_old=', odis_ndg_old(i0,j0), ' km odis_ndg_new=', odis_ndg_new(i0,j0), &
1125 ' km blw_old=', blw_old(i0,j0), ' blw_new=', blw_new(i0,j0)
1126 CALL wrf_message( TRIM(message) )
1127 ENDIF
1128 ENDIF
1129 ENDIF
1131 ENDIF
1132 !
1133 ! tfac_sfc for surface analysis nudging
1134 !
1135 IF( xtime >= actual_end_fdda_min_sfc-ABS(dtramp_min) .AND. xtime <= actual_end_fdda_min_sfc &
1136 .AND. dtramp_min > 0.0 .AND. if_ramping == 1 ) &
1138 !BPR BEGIN
1139 !The original method assumed that to be here in the code with dtramp_min>0
1140 !meant that we were after the valid time of *_ndg_new and so used the
1141 !values *_ndg_old and *_ndg_new to extrapolate a current analysis value.
1142 !HOWEVER, in practice once we get to the valid time of *_ndg_new WRF will
1143 !read in the next pair of *_ndg_old and *_ndg_new values, even if we are
1144 !at the beginning of the rampdown period. Since the *_ndg_new values have
1145 !a valid time after the beginning of the rampdown period we do not want
1146 !to use them as we would be using analyses valid after the supposed end
1147 !time of the analysis nudging.
1148 !coef = (xtime-xtime_old_sfc+analysis_interval_sfc)/(analysis_interval_sfc*1.0)
1149 !Use only the *_ndg_old values
1150 coef = 0.0
1151 !BPR END
1153 ! print*, 'coef =', xtime_old_sfc, xtime, xtime_new_sfc, coef
1154 !
1155 ! Compute surface analysis nudging tendencies for u, v, t and q
1156 !
1157 ! FASDAS
1158 !
1159 IF( fasdas == 1 ) THEN
1161 !Edit TWG2015 Add new variable wpblfasdas to avoid altering global wpbl
1162 !field when using the FASDAS option
1164 wpblfasdas = 1.0
1165 DO j=jts,jtf
1166 DO k=kts,ktf
1167 DO i=its,itf
1168 if( k == 1 ) then
1169 wpblfasdas(i, k, j) = 0.0
1170 else
1171 wpblfasdas(i, k, j) = 1.0
1172 endif
1173 ENDDO
1174 ENDDO
1175 ENDDO
1177 ENDIF
1178 !
1179 ! END FASDAS
1180 !
1181 DO j=jts,jtf
1182 DO k=kts,ktf
1183 DO i=itsu,itf
1184 IF( idd == 1 ) THEN
1185 blw = 0.5* (blw_old(i-1,j)+blw_old(i,j)) * ( 1.0 - coef ) &
1186 + 0.5* (blw_new(i-1,j)+blw_new(i,j)) * coef
1187 ELSE
1188 blw = 1.0
1189 ENDIF
1190 val_analysis = u10_ndg_old(i,j) *( 1.0 - coef ) + u10_ndg_new(i,j) * coef
1191 val_analysis = val_analysis * wndcor_u(i,j)
1192 RUNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,1)) * wzfac(k,1) * tfac_sfc * blw * &
1193 ( val_analysis - u3d(i,1,j) )
1194 ENDDO
1195 ENDDO
1196 ENDDO
1198 DO j=jtsv,jtf
1199 DO k=kts,ktf
1200 DO i=its,itf
1202 IF( idd == 1 ) THEN
1203 blw = 0.5* (blw_old(i,j-1)+blw_old(i,j)) * ( 1.0 - coef ) &
1204 + 0.5* (blw_new(i,j-1)+blw_new(i,j)) * coef
1205 ELSE
1206 blw = 1.0
1207 ENDIF
1208 val_analysis = v10_ndg_old(i,j) *( 1.0 - coef ) + v10_ndg_new(i,j) * coef
1209 val_analysis = val_analysis * wndcor_v(i,j)
1210 RVNDGDTEN(i,k,j) = guv_sfc * (1.0-wpbl(i,k,j,2)) * wzfac(k,2) * tfac_sfc * blw * &
1211 ( val_analysis - v3d(i,1,j) )
1212 ENDDO
1213 ENDDO
1214 ENDDO
1216 DO j=jts,jtf
1217 DO k=kts,ktf
1218 DO i=its,itf
1219 IF( idd == 1 ) THEN
1220 blw = blw_old(i,j) * ( 1.0 - coef ) + blw_new(i,j) * coef
1221 ELSE
1222 blw = 1.0
1223 ENDIF
1225 !ckay2015 and TWG2015 add stability factor for stable regime, ensure
1226 !consistency between Direct and Indirect nudging, and Convert potential
1227 !temperature tendency to temperature tendency
1228 !
1229 ! FASDAS
1230 !
1231 IF( fasdas == 1 ) THEN
1232 if(regime(i,j).le.1.1) then
1233 stabFac = 1.0/3.0
1234 else
1235 stabFac = 1.0
1236 end if
1238 val_analysis = th2_ndg_old(i,j) *( 1.0 - coef ) + th2_ndg_new(i,j) * coef
1239 !BPR BEGIN
1240 !RTHNDGDTEN(i,k,j) = stabFac*gt_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,3) * tfac * blw * &
1241 RTHNDGDTEN(i,k,j) = stabFac*gt_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,3) * tfac_sfc * blw * &
1242 !BPR END
1243 ( val_analysis - th3d(i,1,j))
1245 DiffTN = val_analysis - th3d(i,1,j)
1246 if(k.eq.1) then
1247 exf = (1.0E05/p3d(i,k,j))**(287./1004.)
1248 SDA_HFX(i,j) = RTHNDGDTEN(i,k,j)/exf !TWG 2015
1249 else
1250 RTHNDGDTEN(i,k,j) = 0.0
1251 end if
1253 val_analysis = q2_ndg_old(i,j) *( 1.0 - coef ) + q2_ndg_new(i,j) * coef
1254 IF( iqsat == 1 .AND. val_analysis > qsat(i,k,j) ) val_analysis = qsat(i,k,j)
1255 !BPR BEGIN
1256 !RQVNDGDTEN(i,k,j) = stabFac*gq_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,4) * tfac * blw * &
1257 RQVNDGDTEN(i,k,j) = stabFac*gq_sfc * (1.0-wpblfasdas(i,k,j)) * wzfac(k,4) * tfac_sfc * blw * &
1258 !BPR END
1259 ( val_analysis - qv3d(i,k,j) )
1260 if(k.eq.1) then
1261 SDA_QFX(i,j) = RQVNDGDTEN(i,k,j)
1262 else
1263 RQVNDGDTEN(i,k,j) = 0.0
1264 end if
1266 ELSE
1268 val_analysis = th2_ndg_old(i,j) *( 1.0 - coef ) + th2_ndg_new(i,j) * coef
1269 !BPR BEGIN
1270 !RTHNDGDTEN(i,k,j) = gt_sfc * (1.0-wpbl(i,k,j,3)) * wzfac(k,3) * tfac * blw * &
1271 RTHNDGDTEN(i,k,j) = gt_sfc * (1.0-wpbl(i,k,j,3)) * wzfac(k,3) * tfac_sfc * blw * &
1272 !BPR END
1273 ( val_analysis - th3d(i,1,j))
1275 val_analysis = q2_ndg_old(i,j) *( 1.0 - coef ) + q2_ndg_new(i,j) * coef
1276 IF( iqsat == 1 .AND. val_analysis > qsat(i,k,j) ) val_analysis = qsat(i,k,j)
1277 !BPR BEGIN
1278 !RQVNDGDTEN(i,k,j) = gq_sfc * (1.0-wpbl(i,k,j,4)) * wzfac(k,4) * tfac * blw * &
1279 RQVNDGDTEN(i,k,j) = gq_sfc * (1.0-wpbl(i,k,j,4)) * wzfac(k,4) * tfac_sfc * blw * &
1280 !BPR END
1281 ( val_analysis - qv3d(i,k,j) )
1283 ENDIF
1284 !
1285 ! END FASDAS
1286 !
1287 ENDDO
1288 ENDDO
1289 ENDDO
1291 !BPR BEGIN
1292 !Diagnostic print
1293 IF ( wrf_dm_on_monitor()) THEN
1294 IF( dbg_level .GE. 10 ) THEN
1295 i0 = (ite-its)/2+its
1296 j0 = (jte-jts)/2+jts
1297 WRITE(message,'(a,i1,a,f7.2,a,i4,a,i4,a,f7.2,a,f4.2,a,f7.2,a,f4.2,a,f4.2 )') &
1298 ' D0',id,' At xtime=',xtime,' SFC FDDA sample TH2 analysis (i=',i0,', j=',j0,') = (th2_ndg_old=', &
1299 th2_ndg_old(i0,j0),') * ',1.0-coef,' + (th2_ndg_new=',th2_ndg_new(i0,j0),') * ',coef, &
1300 ' where tfac_sfc=',tfac_sfc
1301 CALL wrf_message( TRIM(message) )
1302 WRITE(message,'(a,i1,a,f7.2,a,i4,a,f7.2,a,f7.2)') &
1303 ' D0',id,' xtime_old_sfc=',xtime_old_sfc,' analysis_interval_sfc=',analysis_interval_sfc, &
1304 ' actual_end_fdda_min_sfc=', actual_end_fdda_min_sfc,' dtramp_min=',dtramp_min
1305 CALL wrf_message( TRIM(message) )
1306 END IF
1307 END IF
1308 !BPR END
1310 END SUBROUTINE sfddagd
1312 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1314 SUBROUTINE fddagdinit(id,rundgdten,rvndgdten,rthndgdten,rqvndgdten,rmundgdten,&
1315 SDA_HFX, SDA_QFX, QNORM, HFX_BOTH, QFX_BOTH, fasdas,& !fasdas
1316 HFX_FDDA, & !fasdas
1317 run_hours, &
1318 if_no_pbl_nudging_uv, if_no_pbl_nudging_t, if_no_pbl_nudging_q, &
1319 if_zfac_uv, k_zfac_uv, if_zfac_t, k_zfac_t, if_zfac_q, k_zfac_q, &
1320 guv, gt, gq, if_ramping, dtramp_min, end_fdda_hour, &
1321 grid_sfdda, guv_sfc, gt_sfc, gq_sfc, &
1322 restart, allowed_to_read, &
1323 ids, ide, jds, jde, kds, kde, &
1324 ims, ime, jms, jme, kms, kme, &
1325 its, ite, jts, jte, kts, kte )
1326 !-------------------------------------------------------------------
1327 IMPLICIT NONE
1328 !-------------------------------------------------------------------
1329 !
1330 INTEGER , INTENT(IN) :: id
1331 LOGICAL, INTENT(IN) :: restart, allowed_to_read
1332 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1333 ims, ime, jms, jme, kms, kme, &
1334 its, ite, jts, jte, kts, kte
1335 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
1336 rundgdten, &
1337 rvndgdten, &
1338 rthndgdten, &
1339 rqvndgdten
1340 !
1341 ! FASDAS
1342 !
1343 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: &
1344 SDA_HFX, &
1345 SDA_QFX, &
1346 QNORM, HFX_BOTH, QFX_BOTH
1347 INTEGER, INTENT(IN ) :: fasdas
1348 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(OUT) :: &
1349 HFX_FDDA
1350 !
1351 ! END FASDAS
1352 !
1353 INTEGER, INTENT(IN) :: run_hours
1354 INTEGER, INTENT(IN) :: if_no_pbl_nudging_uv, if_no_pbl_nudging_t, &
1355 if_no_pbl_nudging_q, end_fdda_hour
1356 INTEGER, INTENT(IN) :: if_zfac_uv, if_zfac_t, if_zfac_q
1357 INTEGER, INTENT(IN) :: k_zfac_uv, k_zfac_t, k_zfac_q
1358 INTEGER, INTENT(IN) :: if_ramping, grid_sfdda
1359 REAL, INTENT(IN) :: dtramp_min
1360 REAL, INTENT(IN) :: guv, gt, gq
1361 REAL, INTENT(IN) :: guv_sfc, gt_sfc, gq_sfc
1362 REAL :: actual_end_fdda_min
1364 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT) :: rmundgdten
1365 INTEGER :: i, j, k
1367 LOGICAL , EXTERNAL :: wrf_dm_on_monitor
1369 CHARACTER (LEN=256) :: message
1371 IF ( wrf_dm_on_monitor() ) THEN
1373 IF( guv > 0.0 ) THEN
1374 WRITE(message,'(a,i1,a,e12.4)') &
1375 'D0',id,' 3-D analysis nudging for wind is applied and Guv= ', guv
1376 CALL wrf_message(TRIM(message))
1377 ELSE IF( guv < 0.0 ) THEN
1378 CALL wrf_error_fatal('In grid FDDA, Guv must be positive.')
1379 ELSE
1380 WRITE(message,'(a,i1,a,e12.4)') &
1381 'D0',id,' 3-D analysis nudging for wind is not applied and Guv= ', guv
1382 CALL wrf_message(TRIM(message))
1383 ENDIF
1385 IF( gt > 0.0 ) THEN
1386 WRITE(message,'(a,i1,a,e12.4)') &
1387 'D0',id,' 3-D analysis nudging for temperature is applied and Gt= ', gt
1388 CALL wrf_message(TRIM(message))
1389 ELSE IF( gt < 0.0 ) THEN
1390 CALL wrf_error_fatal('In grid FDDA, Gt must be positive.')
1391 ELSE
1392 WRITE(message,'(a,i1,a,e12.4)') &
1393 'D0',id,' 3-D analysis nudging for temperature is not applied and Gt= ', gt
1394 CALL wrf_message(TRIM(message))
1395 ENDIF
1397 IF( gq > 0.0 ) THEN
1398 WRITE(message,'(a,i1,a,e12.4)') &
1399 'D0',id,' 3-D analysis nudging for water vapor mixing ratio is applied and Gq= ', gq
1400 CALL wrf_message(TRIM(message))
1401 ELSE IF( gq < 0.0 ) THEN
1402 CALL wrf_error_fatal('In grid FDDA, Gq must be positive.')
1403 ELSE
1404 WRITE(message,'(a,i1,a,e12.4)') &
1405 'D0',id,' 3-D analysis nudging for water vapor mixing ratio is not applied and Gq= ', gq
1406 CALL wrf_message(TRIM(message))
1407 ENDIF
1409 IF( guv > 0.0 .AND. if_no_pbl_nudging_uv == 1 ) THEN
1410 WRITE(message,'(a,i1,a)') &
1411 'D0',id,' 3-D analysis nudging for wind is turned off within the PBL.'
1412 CALL wrf_message(TRIM(message))
1413 ENDIF
1415 IF( gt > 0.0 .AND. if_no_pbl_nudging_t == 1 ) THEN
1416 WRITE(message,'(a,i1,a)') &
1417 'D0',id,' 3-D analysis nudging for temperature is turned off within the PBL.'
1418 CALL wrf_message(TRIM(message))
1419 ENDIF
1421 IF( gq > 0.0 .AND. if_no_pbl_nudging_q == 1 ) THEN
1422 WRITE(message,'(a,i1,a)') &
1423 'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off within the PBL.'
1424 CALL wrf_message(TRIM(message))
1425 ENDIF
1427 IF( guv > 0.0 .AND. if_zfac_uv == 1 ) THEN
1428 WRITE(message,'(a,i1,a,i3)') &
1429 'D0',id,' 3-D analysis nudging for wind is turned off below layer', k_zfac_uv
1430 CALL wrf_message(TRIM(message))
1431 ENDIF
1433 IF( gt > 0.0 .AND. if_zfac_t == 1 ) THEN
1434 WRITE(message,'(a,i1,a,i3)') &
1435 'D0',id,' 3-D analysis nudging for temperature is turned off below layer', k_zfac_t
1436 CALL wrf_message(TRIM(message))
1437 ENDIF
1439 IF( gq > 0.0 .AND. if_zfac_q == 1 ) THEN
1440 WRITE(message,'(a,i1,a,i3)') &
1441 'D0',id,' 3-D analysis nudging for water vapor mixing ratio is turned off below layer', &
1442 k_zfac_q
1443 CALL wrf_message(TRIM(message))
1444 ENDIF
1446 IF( grid_sfdda >=1 ) THEN
1447 IF( guv_sfc > 0.0 ) THEN
1448 WRITE(message,'(a,i1,a,e12.4)') &
1449 'D0',id,' surface analysis nudging for wind is applied and Guv_sfc= ', guv_sfc
1450 CALL wrf_message(TRIM(message))
1451 ELSE IF( guv_sfc < 0.0 ) THEN
1452 CALL wrf_error_fatal('In grid FDDA, Guv_sfc must be positive.')
1453 ELSE
1454 WRITE(message,'(a,i1,a,e12.4)') &
1455 'D0',id,' surface analysis nudging for wind is not applied and Guv_sfc= ', guv_sfc
1456 CALL wrf_message(TRIM(message))
1457 ENDIF
1459 IF( gt_sfc > 0.0 ) THEN
1460 WRITE(message,'(a,i1,a,e12.4)') &
1461 'D0',id,' surface analysis nudging for temperature is applied and Gt_sfc= ', gt_sfc
1462 CALL wrf_message(TRIM(message))
1463 ELSE IF( gt_sfc < 0.0 ) THEN
1464 CALL wrf_error_fatal('In grid FDDA, Gt_sfc must be positive.')
1465 ELSE
1466 WRITE(message,'(a,i1,a,e12.4)') &
1467 'D0',id,' surafce analysis nudging for temperature is not applied and Gt_sfc= ', gt_sfc
1468 CALL wrf_message(TRIM(message))
1469 ENDIF
1471 IF( gq_sfc > 0.0 ) THEN
1472 WRITE(message,'(a,i1,a,e12.4)') &
1473 'D0',id,' surface analysis nudging for water vapor mixing ratio is applied and Gq_sfc= ', gq_sfc
1474 CALL wrf_message(TRIM(message))
1475 ELSE IF( gq_sfc < 0.0 ) THEN
1476 CALL wrf_error_fatal('In grid FDDA, Gq_sfc must be positive.')
1477 ELSE
1478 WRITE(message,'(a,i1,a,e12.4)') &
1479 'D0',id,' surface analysis nudging for water vapor mixing ratio is not applied and Gq_sfc= ', gq_sfc
1480 CALL wrf_message(TRIM(message))
1481 ENDIF
1483 ENDIF
1485 IF( if_ramping == 1 .AND. ABS(dtramp_min) > 0.0 ) THEN
1486 IF( dtramp_min <= 0.0 ) THEN
1487 actual_end_fdda_min = end_fdda_hour*60.0
1488 ELSE
1489 actual_end_fdda_min = end_fdda_hour*60.0 + ABS(dtramp_min)
1490 ENDIF
1492 IF( actual_end_fdda_min <= run_hours*60. ) THEN
1493 WRITE(message,'(a,i1,a)') &
1494 'D0',id,' analysis nudging is ramped down near the end of the nudging period,'
1495 CALL wrf_message(TRIM(message))
1497 WRITE(message,'(a,f6.2,a,f6.2,a)') &
1498 ' starting at ', (actual_end_fdda_min - ABS(dtramp_min))/60.0, &
1499 'h, ending at ', actual_end_fdda_min/60.0,'h.'
1500 CALL wrf_message(TRIM(message))
1501 ENDIF
1502 ENDIF
1504 ENDIF
1506 IF(.not.restart) THEN
1507 DO j = jts,jte
1508 DO k = kts,kte
1509 DO i = its,ite
1510 rundgdten(i,k,j) = 0.
1511 rvndgdten(i,k,j) = 0.
1512 rthndgdten(i,k,j) = 0.
1513 rqvndgdten(i,k,j) = 0.
1514 ! TWG 2015 Psuedo radiative flux
1515 HFX_FDDA(i,k,j) = 0.
1516 ! TWG END
1517 if(k.eq.kts) rmundgdten(i,j) = 0.
1518 ENDDO
1519 ENDDO
1520 ENDDO
1521 !
1522 ! FASDAS
1523 !
1524 IF( fasdas == 1 ) THEN
1525 DO j = jts,jte
1526 DO i = its,ite
1527 SDA_HFX(I,J) = 0.0
1528 SDA_QFX(I,J) = 0.0
1529 QNORM(I,J) = 0.0
1530 HFX_BOTH(I,J) = 0.0
1531 QFX_BOTH(I,J) = 0.0
1532 ENDDO
1533 ENDDO
1534 ENDIF
1535 !
1536 ! END FASDAS
1537 !
1538 ENDIF
1540 END SUBROUTINE fddagdinit
1542 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1543 SUBROUTINE smther(slab, idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd)
1544 !
1545 ! PURPOSE: SPATIALLY SMOOTH DATA IN SLAB TO DAMPEN SHORT
1546 ! WAVELENGTH COMPONENTS.
1547 !
1548 ! Implemented based on the same smoothing subroutine used in MM5, with modifications to
1549 ! remove the extra loop that causes unneccessary desmoothing. Aijun Deng (Penn State),
1550 ! December 2008
1551 !
1552 IMPLICIT NONE
1554 INTEGER :: idimst, idimnd, jdimst, jdimnd, npass, ist, ind, jst, jnd
1555 INTEGER :: i, j, k, kk
1556 REAL :: avg
1557 REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB
1558 REAL, DIMENSION(idimst:idimnd, jdimst:jdimnd) :: SLAB_ORIG
1559 REAL, DIMENSION(2) :: XNU
1561 IF(NPASS.EQ.0)RETURN
1562 XNU(1)=0.50
1563 XNU(2)=-0.52
1564 DO K=1,NPASS
1565 KK = 2 - MOD(K,2)
1567 DO J=JDIMST,JDIMND
1568 DO I=IDIMST,IDIMND
1569 SLAB_ORIG(I,J) = SLAB(I,J)
1570 END DO
1571 END DO
1573 DO J=JST,JND
1574 DO I=IST,IND
1575 AVG = ( SLAB_ORIG(I+1,J ) + &
1576 SLAB_ORIG(I-1,J ) + &
1577 SLAB_ORIG(I ,J+1) + &
1578 SLAB_ORIG(I ,J-1) ) * 0.25
1579 SLAB(I,J)=SLAB_ORIG(I,J)+XNU(KK)*(AVG - SLAB_ORIG(I,J))
1580 ENDDO
1581 ENDDO
1583 ENDDO
1585 END SUBROUTINE smther
1587 !-------------------------------------------------------------------
1588 END MODULE module_fdda_psufddagd