| blob | 628fb93acc18f9a4a91e5a1c3612701ae7b34024 |
1 #if (NMM_CORE == 1)
2 MODULE module_diag_afwa
3 CONTAINS
4 SUBROUTINE diag_afwa_stub
5 END SUBROUTINE diag_afwa_stub
6 END MODULE module_diag_afwa
7 #else
9 !WRF:MEDIATION_LAYER:PHYSICS
11 MODULE module_diag_afwa
13 USE diag_functions
15 CONTAINS
17 SUBROUTINE afwa_diagnostics_driver ( grid , config_flags &
18 , moist &
19 , scalar &
20 , chem &
21 , th_phy , pi_phy , p_phy &
22 , u_phy , v_phy &
23 , dz8w , p8w , t8w , rho &
24 , ids, ide, jds, jde, kds, kde &
25 , ims, ime, jms, jme, kms, kme &
26 , ips, ipe, jps, jpe, kps, kpe &
27 , its, ite, jts, jte &
28 , k_start, k_end )
30 USE module_domain, ONLY : domain , domain_clock_get
31 USE module_configure, ONLY : grid_config_rec_type, model_config_rec
32 USE module_state_description
33 USE module_model_constants
34 USE module_utility
35 USE module_streams, ONLY: history_alarm, auxhist2_alarm
36 #ifdef DM_PARALLEL
37 USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval
38 #endif
40 IMPLICIT NONE
42 TYPE ( domain ), INTENT(INOUT) :: grid
43 TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags
45 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
46 ims, ime, jms, jme, kms, kme, &
47 ips, ipe, jps, jpe, kps, kpe
48 INTEGER :: k_start , k_end, its, ite, jts, jte
50 REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_moist), &
51 INTENT(IN ) :: moist
53 REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_scalar), &
54 INTENT(IN ) :: scalar
56 REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_chem), &
57 INTENT(IN ) :: chem
59 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
60 INTENT(IN ) :: th_phy &
61 , pi_phy &
62 , p_phy &
63 , u_phy &
64 , v_phy &
65 , dz8w &
66 , p8w &
67 , t8w &
68 , rho
70 ! Local
71 ! -----
72 CHARACTER*512 :: message
73 CHARACTER*256 :: timestr
74 INTEGER :: i,j,k,nz,ostat
75 INTEGER :: icing_opt
76 REAL :: bdump
77 INTEGER :: i_start, i_end, j_start, j_end
78 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qrain &
79 , qsnow &
80 , qgrpl &
81 , qvapr &
82 , qcloud &
83 , qice &
84 , ncloud &
85 , ngraup &
86 , rh &
87 , rh_cld &
88 , ptot &
89 , z_e &
90 , zagl
92 REAL, DIMENSION( ims:ime, jms:jme, 5 ) :: dustc
93 REAL, DIMENSION( ims:ime, jms:jme ) :: rh2m &
94 , rh20m &
95 , tv2m &
96 , tv20m &
97 , wind10m &
98 , wup_mask &
99 , wind125m &
100 , llws &
101 , pwater
103 LOGICAL :: do_buoy_calc
104 REAL :: zlfc_msl, dum1, dum2, dum3, wind_vel, wind_blend
105 REAL :: prate_mm_per_hr, factor
106 REAL :: u1km, v1km, ublend, vblend, u2000, v2000, us, vs
107 LOGICAL :: is_target_level
109 ! Timing
110 ! ------
111 TYPE(WRFU_Time) :: hist_time, aux2_time, CurrTime, StartTime
112 TYPE(WRFU_TimeInterval) :: dtint, histint, aux2int
113 LOGICAL :: is_after_history_dump, is_output_timestep, is_first_timestep
114 INTEGER , PARAMETER :: DEBUG_LEVEL = 1
116 ! Chirp the routine name for debugging purposes
117 ! ---------------------------------------------
118 write ( message, * ) 'inside afwa_diagnostics_driver'
119 CALL wrf_debug( 100 , message )
121 ! Get timing info
122 ! Want to know if when the last history output was
123 ! Check history and auxhist2 alarms to check last ring time and how often
124 ! they are set to ring
125 ! -----------------------------------------------------------------------
126 CALL WRFU_ALARMGET( grid%alarms( HISTORY_ALARM ), prevringtime=hist_time, &
127 ringinterval=histint)
128 CALL WRFU_ALARMGET( grid%alarms( AUXHIST2_ALARM ), prevringtime=aux2_time, &
129 ringinterval=aux2int)
131 ! Get domain clock
132 ! ----------------
133 CALL domain_clock_get ( grid, current_time=CurrTime, &
134 simulationStartTime=StartTime, &
135 current_timestr=timestr, time_step=dtint )
137 ! Set some booleans for use later
138 ! Following uses an overloaded .lt.
139 ! ---------------------------------
140 is_after_history_dump = ( Currtime .lt. hist_time + dtint )
142 ! Following uses an overloaded .ge.
143 ! ---------------------------------
144 is_output_timestep = (Currtime .ge. hist_time + histint - dtint .or. &
145 Currtime .ge. aux2_time + aux2int - dtint )
146 write ( message, * ) 'is output timestep? ', is_output_timestep
147 CALL wrf_debug( 100 , message )
149 ! Following uses an overloaded .eq.
150 ! ---------------------------------
151 is_first_timestep = ( Currtime .eq. StartTime + dtint )
153 ! Here is an optional check for bad data in case the model has gone
154 ! off the hinges unchecked until now. This happens under certain
155 ! configurations and can lead to very wrong data writes. This is just
156 ! a simple check for sane winds and potential temperatures.
157 ! --------------------------------------------------------------------
158 IF ( config_flags%afwa_bad_data_check .GT. 0 ) THEN
159 IF ( ( is_output_timestep ) .AND. ( .NOT. is_first_timestep ) ) THEN
160 DO i=its, MIN( ide-1, ite )
161 DO k=k_start, k_end
162 DO j=jts, MIN( jde-1, jte )
163 IF ( ( u_phy(i,k,j) .GT. 300. ) .OR. &
164 ( u_phy(i,k,j) .LT. -300. ) .OR. &
165 ( v_phy(i,k,j) .GT. 300. ) .OR. &
166 ( v_phy(i,k,j) .LT. -300. ) .OR. &
167 ( th_phy(i,k,j) .GT. 9999. ) .OR. &
168 ( th_phy(i,k,j) .LT. 99. ) ) THEN
169 write ( message, * ) "AFWA Diagnostics: ERROR - Model winds and/or " // &
170 "potential temperature appear to be bad. If you do not want this check, " // &
171 "set afwa_bad_data_check=0. i=",i,", j=",j,", k=",k,", u_phy=",u_phy(i,k,j), &
172 ", v_phy=", v_phy(i,k,j),", th_phy=",th_phy(i,k,j)
173 CALL wrf_error_fatal( message )
174 ENDIF
175 ENDDO
176 ENDDO
177 ENDDO
178 ENDIF
179 ENDIF
181 ! 3-D arrays for moisture variables
182 ! ---------------------------------
184 IF ( F_QV ) THEN
185 DO i=ims, ime
186 DO k=kms, kme
187 DO j=jms, jme
188 qvapr(i,k,j) = moist(i,k,j,P_QV)
189 ENDDO
190 ENDDO
191 ENDDO
192 END IF
193 IF ( F_QR ) THEN
194 DO i=ims, ime
195 DO k=kms, kme
196 DO j=jms, jme
197 qrain(i,k,j) = moist(i,k,j,P_QR)
198 ENDDO
199 ENDDO
200 ENDDO
201 END IF
202 IF ( F_QS ) THEN
203 DO i=ims, ime
204 DO k=kms, kme
205 DO j=jms, jme
206 qsnow(i,k,j) = moist(i,k,j,P_QS)
207 ENDDO
208 ENDDO
209 ENDDO
210 END IF
211 IF ( F_QG ) THEN
212 DO i=ims, ime
213 DO k=kms, kme
214 DO j=jms, jme
215 qgrpl(i,k,j) = moist(i,k,j,P_QG)
216 ENDDO
217 ENDDO
218 ENDDO
219 END IF
220 IF ( F_QC ) THEN
221 DO i=ims, ime
222 DO k=kms, kme
223 DO j=jms, jme
224 qcloud(i,k,j) = moist(i,k,j,P_QC)
225 ENDDO
226 ENDDO
227 ENDDO
228 END IF
229 IF ( F_QI ) THEN
230 DO i=ims, ime
231 DO k=kms, kme
232 DO j=jms, jme
233 qice(i,k,j) = moist(i,k,j,P_QI)
234 ENDDO
235 ENDDO
236 ENDDO
237 END IF
238 IF ( F_QNC ) THEN
239 DO i=ims, ime
240 DO k=kms, kme
241 DO j=jms, jme
242 ncloud(i,k,j) = scalar(i,k,j,P_QNC)
243 ENDDO
244 ENDDO
245 ENDDO
246 END IF
247 IF ( F_QNG ) THEN
248 DO i=ims, ime
249 DO k=kms, kme
250 DO j=jms, jme
251 ngraup(i,k,j) = scalar(i,k,j,P_QNG)
252 ENDDO
253 ENDDO
254 ENDDO
255 END IF
257 ! Total pressure
258 ! --------------
259 DO i=ims, ime
260 DO k=kms, kme
261 DO j=jms, jme
262 ptot(i,k,j)=grid%pb(i,k,j)+grid%p(i,k,j)
263 ENDDO
264 ENDDO
265 ENDDO
267 ! ZAGL (height above ground)
268 ! --------------------------
269 DO i=ims, ime
270 DO k=kms, kme
271 DO j=jms, jme
272 zagl(i,k,j)=grid%z(i,k,j)-grid%ht(i,j)
273 ENDDO
274 ENDDO
275 ENDDO
277 ! Calculate relative humidity
278 ! ---------------------------
279 IF ( F_QV ) THEN
280 DO i=ims,ime
281 DO k=kms,kme
282 DO j=jms,jme
283 rh(i,k,j)=calc_rh(ptot(i,k,j),grid%t_phy(i,k,j), qvapr(i,k,j))
284 ENDDO
285 ENDDO
286 ENDDO
287 ELSE
288 CALL wrf_debug ( DEBUG_LEVEL, 'Option rh requires: QV' )
289 END IF
290 IF ( F_QV .AND. F_QC .AND. F_QI ) THEN
291 DO i=ims,ime
292 DO k=kms,kme
293 DO j=jms,jme
294 rh_cld(i,k,j)=calc_rh(ptot(i,k,j),grid%t_phy(i,k,j), qvapr(i,k,j)+qcloud(i,k,j)+qice(i,k,j))
295 ENDDO
296 ENDDO
297 ENDDO
298 ELSE
299 CALL wrf_debug ( DEBUG_LEVEL, 'Option rh_cld requires: QV, QC, QI' )
300 END IF
302 ! Time-step precipitation (convective + nonconvective)
303 ! --------------------------------------------------------------
304 DO i=ims,ime
305 DO j=jms,jme
306 grid % afwa_precip(i,j) = grid%raincv(i,j) + grid%rainncv(i,j)
307 grid % afwa_totprecip(i,j) = grid%rainc(i,j) + grid%rainnc(i,j)
308 ENDDO
309 ENDDO
311 ! Calculate precipitable water
312 ! ----------------------------
313 IF ( F_QV .AND. F_QC ) THEN
314 nz=kme-kms+1
315 DO i=ims,ime
316 DO j=jms,jme
317 grid % afwa_pwat ( i, j ) = Pwat( nz, &
318 qvapr(i,kms:kme,j), &
319 qcloud(i,kms:kme,j), &
320 dz8w(i,kms:kme,j), &
321 rho(i,kms:kme,j) )
322 ENDDO
323 ENDDO
324 ELSE
325 CALL wrf_debug ( DEBUG_LEVEL, 'Option pwat requires: QV, QC' )
326 END IF
328 ! After each history dump, reset max/min value arrays
329 ! ----------------------------------------------------------------------
330 IF ( is_after_history_dump ) THEN
331 IF ( config_flags % afwa_severe_opt == 1 ) THEN
332 DO j = jms, jme
333 DO i = ims, ime
334 grid % wspd10max(i,j) = 0.
335 grid % afwa_llws(i,j) = 0.
336 ENDDO
337 ENDDO
338 ENDIF
339 ENDIF
341 ! Calculate the max 10 m wind speed between output times
342 ! ------------------------------------------------------
343 ! UPDATE 20150112 - GAC
344 ! Diagnose from model 10 m winds, and blend with 1 km AGL
345 ! winds when precipitation rate is > 50 mm/hr to account
346 ! for increased surface wind gust potential when precip
347 ! is heavy and when winds aloft are strong. Will use the
348 ! higher of the surface and the blended winds. Blending
349 ! is linear weighted between 50-150 mm/hr precip rates.
350 ! -------------------------------------------------------
351 IF ( config_flags % afwa_severe_opt == 1 ) THEN
352 DO j = jms, jme
353 DO i = ims, ime
354 wind_vel = uv_wind ( grid % u10(i,j) , grid % v10(i,j) )
355 prate_mm_per_hr = ( grid % afwa_precip(i,j) / grid % dt ) * 3600.
357 ! Is this an area of heavy precip? Calculate 1km winds to blend down
358 ! -------------------------------------------------------------------
359 IF ( prate_mm_per_hr .GT. 50. ) THEN
360 is_target_level=.false.
361 DO k=kms,kme
362 IF ( ( zagl(i,k,j) >= 1000. ) .and. &
363 ( .NOT. is_target_level ) .and. &
364 ( k .ne. kms ) ) THEN
365 is_target_level = .true.
366 u1km = u_phy(i,k-1,j) + (1000. - (zagl(i,k-1,j))) &
367 * ((u_phy(i,k,j) - u_phy(i,k-1,j))/(zagl(i,k,j)))
368 v1km = v_phy(i,k-1,j) + (1000. - (zagl(i,k-1,j))) &
369 * ((v_phy(i,k,j) - v_phy(i,k-1,j))/(zagl(i,k,j)))
370 EXIT ! We've found our level, break the loop
371 ENDIF
372 ENDDO
374 ! Compute blended wind
375 ! --------------------
376 factor = MAX ( ( ( 150. - prate_mm_per_hr ) / 100. ), 0. )
377 ublend = grid % u10(i,j) * factor + u1km * (1. - factor)
378 vblend = grid % v10(i,j) * factor + v1km * (1. - factor)
379 wind_blend = uv_wind ( ublend, vblend )
381 ! Set the surface wind to the blended wind if higher
382 ! --------------------------------------------------
383 IF ( wind_blend .GT. wind_vel ) THEN
384 wind_vel = wind_blend
385 ENDIF
386 ENDIF
388 IF ( wind_vel .GT. grid % wspd10max(i,j) ) THEN
389 grid % wspd10max(i,j) = wind_vel
390 ENDIF
391 ENDDO
392 ENDDO
393 ENDIF
395 ! Calculate 0-2000 foot (0 - 609.6 meter) shear.
396 ! ----------------------------------------------
397 IF ( config_flags % afwa_severe_opt == 1 ) THEN
398 DO j = jts, jte
399 DO i = its, ite
400 is_target_level=.false.
401 DO k=kms,kme
402 IF ( ( zagl(i,k,j) >= 609.6 ) .and. &
403 ( .NOT. is_target_level ) .and. &
404 ( k .ne. kms ) ) THEN
405 is_target_level = .true.
406 u2000 = u_phy(i,k-1,j) + (609.6 - (zagl(i,k-1,j))) &
407 * ((u_phy(i,k,j) - u_phy(i,k-1,j))/(zagl(i,k,j)))
408 v2000 = v_phy(i,k-1,j) + (609.6 - (zagl(i,k-1,j))) &
409 * ((v_phy(i,k,j) - v_phy(i,k-1,j))/(zagl(i,k,j)))
410 us = u2000 - grid % u10(i,j)
411 vs = v2000 - grid % v10(i,j)
412 llws(i,j) = uv_wind ( us , vs )
413 IF ( llws(i,j) .gt. grid % afwa_llws(i,j) ) THEN
414 grid % afwa_llws(i,j) = llws(i,j)
415 ENDIF
416 EXIT ! We've found our level, break the loop
417 ENDIF
418 ENDDO
419 ENDDO
420 ENDDO
421 ENDIF
423 #if ( WRF_CHEM == 1 )
424 ! Surface dust concentration array (ug m-3)
425 ! -----------------------------------------
426 DO i=ims, ime
427 DO j=jms, jme
428 dustc(i,j,1)=chem(i,k_start,j,p_dust_1)*rho(i,k_start,j)
429 dustc(i,j,2)=chem(i,k_start,j,p_dust_2)*rho(i,k_start,j)
430 dustc(i,j,3)=chem(i,k_start,j,p_dust_3)*rho(i,k_start,j)
431 dustc(i,j,4)=chem(i,k_start,j,p_dust_4)*rho(i,k_start,j)
432 dustc(i,j,5)=chem(i,k_start,j,p_dust_5)*rho(i,k_start,j)
433 ENDDO
434 ENDDO
435 #else
436 dustc(ims:ime,jms:jme,:)=0.
437 #endif
439 ! Calculate severe weather diagnostics. These variables should only be
440 ! output at highest frequency output. (e.g. auxhist2)
441 ! ---------------------------------------------------------------------
442 IF ( config_flags % afwa_severe_opt == 1 ) THEN
444 ! After each history dump, reset max/min value arrays
445 ! Note: This resets up_heli_max which is currently calculated within
446 ! rk_first_rk_step_part2.F, may want to move to this diagnostics package
447 ! later
448 ! ----------------------------------------------------------------------
449 IF ( is_after_history_dump ) THEN
450 DO j = jms, jme
451 DO i = ims, ime
452 grid%w_up_max(i,j) = 0.
453 grid%w_dn_max(i,j) = 0.
454 grid%tcoli_max(i,j) = 0.
455 grid%grpl_flx_max(i,j) = 0.
456 grid%up_heli_max(i,j) = 0.
457 grid%afwa_tornado(i,j) = 0.
458 grid%midrh_min_old(i,j) = grid%midrh_min(i,j) ! Save old midrh_min
459 grid%midrh_min(i,j) = 999.
460 grid%afwa_hail(i,j) = 0.
461 ENDDO
462 ENDDO
463 ENDIF ! is_after_history_dump
465 IF ( ( is_first_timestep ) .OR. ( is_output_timestep ) ) THEN
466 do_buoy_calc = .true.
467 ELSE
468 do_buoy_calc = .false.
469 ENDIF
471 IF ( F_QV .AND. F_QR .AND. F_QC .AND. F_QI .AND. F_QS .AND. F_QG .AND. F_QNG ) THEN
472 !-->RAS
473 ! We need to do some neighboring gridpoint comparisons in this next function;
474 ! set these values so we don't go off the edges of the domain. Updraft
475 ! duration on domain edges will always be 0.
476 ! ----------------------------------------------------------------------
477 i_start = its
478 i_end = ite
479 j_start = jts
480 j_end = jte
482 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
483 config_flags%nested) i_start = MAX( ids+1, its )
484 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
485 config_flags%nested) i_end = MIN( ide-1, ite )
486 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
487 config_flags%nested) j_start = MAX( jds+1, jts )
488 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
489 config_flags%nested) j_end = MIN( jde-1, jte )
490 IF ( config_flags%periodic_x ) i_start = its
491 IF ( config_flags%periodic_x ) i_end = ite
493 CALL severe_wx_diagnostics ( grid % wspd10max &
494 , grid % w_up_max &
495 , grid % w_dn_max &
496 , grid % up_heli_max &
497 , grid % tcoli_max &
498 , grid % midrh_min_old &
499 , grid % midrh_min &
500 , grid % afwa_hail &
501 , grid % afwa_cape &
502 , grid % afwa_cin &
503 , grid % afwa_zlfc &
504 , grid % afwa_plfc &
505 , grid % afwa_lidx &
506 , llws &
507 , grid % afwa_tornado &
508 , grid % grpl_flx_max &
509 , grid % u10 &
510 , grid % v10 &
511 , grid % w_2 &
512 , grid % uh &
513 , grid % t_phy &
514 , grid % t2 &
515 , grid % z &
516 , grid % ht &
517 , grid % tornado_mask &
518 , grid % tornado_dur &
519 , grid % dt &
520 , grid % afwa_pwat &
521 , u_phy &
522 , v_phy &
523 , ptot &
524 , qice &
525 , qsnow &
526 , qgrpl &
527 , ngraup &
528 , qvapr, qrain, qcloud &
529 , rho &
530 , dz8w &
531 , rh &
532 , do_buoy_calc &
533 , ims, ime, jms, jme, kms, kme &
534 , its, ite, jts, jte &
535 , k_start, k_end &
536 , j_start, j_end, i_start, i_end )
538 ELSE
539 CALL wrf_debug ( DEBUG_LEVEL, 'Option severe_wx_diagnostics requires: QV, QR, QC, QI, QS, QG, QNG' )
540 END IF
541 ENDIF ! afwa_severe_opt == 1
543 ! Calculate precipitation type diagnostics
544 ! ----------------------------------------
545 IF ( config_flags % afwa_ptype_opt == 1 ) THEN
547 ! First initialize precip buckets
548 ! -------------------------------
549 IF ( grid % itimestep .eq. 1) THEN
550 DO i=ims,ime
551 DO j=jms,jme
552 grid % afwa_rain(i,j)=0.
553 grid % afwa_snow(i,j)=0.
554 grid % afwa_ice(i,j)=0.
555 grid % afwa_fzra(i,j)=0.
556 grid % afwa_snowfall(i,j)=0.
557 ENDDO
558 ENDDO
559 ENDIF
561 ! Diagnose precipitation type
562 ! ---------------------------
563 CALL precip_type_diagnostics ( grid % t_phy &
564 , grid % t2 &
565 , rh &
566 , grid % z &
567 , dz8w &
568 , grid % ht &
569 , grid % afwa_precip &
570 , grid % swdown &
571 , grid % afwa_rain &
572 , grid % afwa_snow &
573 , grid % afwa_ice &
574 , grid % afwa_fzra &
575 , grid % afwa_snowfall &
576 , grid % afwa_ptype_ccn_tmp &
577 , grid % afwa_ptype_tot_melt &
578 , ids, ide, jds, jde, kds, kde &
579 , ims, ime, jms, jme, kms, kme &
580 , ips, ipe, jps, jpe, kps, kpe )
581 ENDIF ! afwa_ptype_opt == 1
583 ! --------------------------------------------------------------
584 ! The following packages are calculated only on output timesteps
585 ! --------------------------------------------------------------
586 IF ( is_output_timestep ) THEN
588 ! Calculate mean sea level pressure
589 ! ---------------------------------
590 DO j = jms, jme
591 DO i = ims, ime
592 grid % afwa_mslp ( i, j ) = MSLP ( grid % ht ( i, j ) &
593 , grid % psfc ( i, j ) &
594 , grid % z ( i, kms, j ) &
595 , qvapr ( i, kms, j ) &
596 , grid % t_phy ( i, kms, j ) )
597 ENDDO
598 ENDDO
600 ! Calculate 10 meter winds
601 ! ------------------------
602 DO i=ims,ime
603 DO j=jms,jme
604 wind10m(i,j)=uv_wind(grid%u10(i,j),grid%v10(i,j))
605 ENDDO
606 ENDDO
608 ! Calculate 2 meter relative humidity/Tv
609 ! --------------------------------------
610 DO i=ims,ime
611 DO j=jms,jme
612 rh2m(i,j)=calc_rh(grid%psfc(i,j), grid%t2(i,j), grid%q2(i,j))
613 tv2m(i,j)=grid%t2(i,j) * (1 + 0.61 * grid%q2(i,j))
614 ENDDO
615 ENDDO
617 ! Calculate buoyancy parameters.
618 ! ------------------------------
619 IF ( config_flags % afwa_buoy_opt == 1 ) THEN
620 nz = k_end - k_start + 1
622 ! Calculate buoyancy for surface parcel
623 ! -------------------------------------
624 DO j = jts, jte
625 DO i = its, ite
626 ostat = Buoyancy ( nz &
627 , grid%t_phy(i,kms:kme,j) &
628 , rh(i,kms:kme,j) &
629 , ptot(i,kms:kme,j) &
630 , grid % z(i,kms:kme,j) &
631 , 1 &
632 , grid % afwa_cape(i,j) &
633 , grid % afwa_cin(i,j) &
634 , zlfc_msl &
635 , grid % afwa_plfc(i,j) &
636 , grid % afwa_lidx(i,j) &
637 , 3 ) !Surface
639 ! Subtract terrain height to convert ZLFC from MSL to AGL
640 ! -------------------------------------------------------
641 IF ( zlfc_msl .ge. grid % ht ( i, j ) ) THEN
642 grid % afwa_zlfc ( i, j ) = zlfc_msl - grid % ht ( i, j )
643 ELSE
644 grid % afwa_zlfc( i, j ) = -1.
645 ENDIF
647 ! Add 273.15 to lifted index per some standard I don't understand
648 ! ---------------------------------------------------------------
649 IF ( grid % afwa_lidx ( i, j ) .ne. 999. ) THEN
650 grid % afwa_lidx ( i, j ) = grid % afwa_lidx ( i, j ) + 273.15
651 ENDIF
653 ! Calculate buoyancy again for most unstable layer
654 ! ------------------------------------------------
655 ostat = Buoyancy ( nz &
656 , grid%t_phy(i,kms:kme,j) &
657 , rh(i,kms:kme,j) &
658 , ptot(i,kms:kme,j) &
659 , grid % z(i,kms:kme,j) &
660 , 1 &
661 , grid % afwa_cape_mu(i,j) &
662 , grid % afwa_cin_mu(i,j) &
663 , dum1 &
664 , dum2 &
665 , dum3 &
666 , 1 ) !Most unstable
668 ENDDO
669 ENDDO
670 ENDIF ! afwa_buoy_opt == 1
672 IF ( config_flags % afwa_therm_opt == 1 ) THEN
673 write ( message, * ) 'Calculating thermal indices'
674 CALL wrf_debug( 100 , message )
675 CALL thermal_diagnostics ( grid % t2 &
676 , grid % psfc &
677 , rh2m &
678 , wind10m &
679 , grid % afwa_heatidx &
680 , grid % afwa_wchill &
681 , grid % afwa_fits &
682 , ids, ide, jds, jde, kds, kde &
683 , ims, ime, jms, jme, kms, kme &
684 , ips, ipe, jps, jpe, kps, kpe )
685 ENDIF ! afwa_therm_opt == 1
687 IF ( config_flags % afwa_turb_opt == 1 ) THEN
688 write ( message, * ) 'Calculating turbulence indices'
690 !~ For now, hard code turbulence layer bottom and top AGL heights
691 ! --------------------------------------------------------------
692 grid % afwa_tlyrbot = (/ 1500., 3000., 4600., 6100., 7600., 9100., &
693 10700. /)
694 grid % afwa_tlyrtop = (/ 3000., 4600., 6100., 7600., 9100., 10700., &
695 12700. /)
696 call turbulence_diagnostics ( u_phy &
697 , v_phy &
698 , grid % t_phy &
699 , ptot &
700 , zagl &
701 , grid % defor11 &
702 , grid % defor12 &
703 , grid % defor22 &
704 , grid % afwa_turb &
705 , grid % afwa_llturb &
706 , grid % afwa_llturblgt &
707 , grid % afwa_llturbmdt &
708 , grid % afwa_llturbsvr &
709 !, config_flags % num_turb_layers &
710 , 7 &
711 , grid % afwa_tlyrbot &
712 , grid % afwa_tlyrtop &
713 , ids, ide, jds, jde, kds, kde &
714 , ims, ime, jms, jme, kms, kme &
715 , ips, ipe, jps, jpe, kps, kpe )
716 ENDIF ! afwa_turb_opt == 1
718 ! Calculate equivalent radar reflectivity factor (z_e) using
719 ! old RIP code (2004) if running radar or VIL packages.
720 ! ----------------------------------------------------------
721 IF ( config_flags % afwa_radar_opt == 1 .or. &
722 config_flags % afwa_vil_opt == 1 ) THEN
723 write ( message, * ) 'Calculating Radar'
724 CALL wrf_debug( 100 , message )
725 CALL wrf_dbzcalc ( rho &
726 , grid%t_phy &
727 , qrain &
728 , qsnow &
729 , qgrpl &
730 , z_e &
731 , ids, ide, jds, jde, kds, kde &
732 , ims, ime, jms, jme, kms, kme &
733 , ips, ipe, jps, jpe, kps, kpe )
734 ENDIF ! afwa_radar_opt == 1 .or. afwa_vil_opt == 1
736 ! Calculate derived radar variables
737 ! ---------------------------------
738 ! UPDATE: removed refd_max calculation, which was not working correctly.
739 ! To correctly calculate refd_max, this routine could be called every
740 ! time step, but instead, we are only going to calculate reflectivity
741 ! on output time steps and avoid cost of calculating refd_max. GAC2014
742 ! ----------------------------------------------------------------------
743 IF ( config_flags % afwa_radar_opt == 1 ) THEN
744 write ( message, * ) 'Calculating derived radar variables'
745 CALL wrf_debug( 100 , message )
746 CALL radar_diagnostics ( grid % refd &
747 , grid % refd_com &
748 ! , grid % refd_max &
749 , grid % echotop &
750 , grid % z &
751 , z_e &
752 , ids, ide, jds, jde, kds, kde &
753 , ims, ime, jms, jme, kms, kme &
754 , ips, ipe, jps, jpe, kps, kpe )
755 ENDIF ! afwa_radar_opt == 1
757 ! Calculate VIL and reflectivity every history output timestep
758 ! ------------------------------------------------------------
759 IF ( config_flags % afwa_vil_opt == 1 ) THEN
760 write ( message, * ) 'Calculating VIL'
761 CALL wrf_debug( 100 , message )
762 CALL vert_int_liquid_diagnostics ( grid % vil &
763 , grid % radarvil &
764 , grid % t_phy &
765 , qrain &
766 , qsnow &
767 , qgrpl &
768 , z_e &
769 , dz8w &
770 , rho &
771 , ids, ide, jds, jde, kds, kde &
772 , ims, ime, jms, jme, kms, kme &
773 , ips, ipe, jps, jpe, kps, kpe )
774 ENDIF ! afwa_vil_opt ==1
776 IF ( F_QR .AND. F_QC .AND. F_QNC ) THEN
777 ! Calculate icing and freezing level
778 ! ----------------------------------
779 IF ( config_flags % afwa_icing_opt == 1 ) THEN
781 ! Determine icing option from microphysics scheme
782 ! -----------------------------------------------
784 IF ( config_flags % mp_physics == GSFCGCESCHEME ) THEN
785 icing_opt=1
786 ELSEIF ( config_flags % mp_physics == ETAMPNEW ) THEN
787 icing_opt=2
788 ELSEIF ( config_flags % mp_physics == THOMPSON ) THEN
789 icing_opt=3
790 ELSEIF ( config_flags % mp_physics == WSM5SCHEME .OR. &
791 config_flags % mp_physics == WSM6SCHEME ) THEN
792 icing_opt=4
793 ELSEIF ( config_flags % mp_physics == MORR_TWO_MOMENT .OR. &
794 config_flags % mp_physics == MORR_TM_AERO ) THEN !TWG add 2017
796 ! Is this run with prognostic cloud droplets or no?
797 ! -------------------------------------------------
798 IF (config_flags % progn > 0) THEN
799 icing_opt=6
800 ELSE
801 icing_opt=5
802 ENDIF
803 ELSEIF ( config_flags % mp_physics == WDM6SCHEME ) THEN
804 icing_opt=7
805 ELSE
806 icing_opt=0 ! Not supported
807 ENDIF
809 IF ( icing_opt .NE. 0 ) THEN
810 write ( message, * ) 'Calculating Icing with icing opt ',icing_opt
811 CALL wrf_debug( 100 , message )
812 CALL icing_diagnostics ( icing_opt &
813 , grid % fzlev &
814 , grid % icing_lg &
815 , grid % icing_sm &
816 , grid % qicing_lg_max &
817 , grid % qicing_sm_max &
818 , grid % qicing_lg &
819 , grid % qicing_sm &
820 , grid % icingtop &
821 , grid % icingbot &
822 , grid % t_phy &
823 , grid % z &
824 , dz8w &
825 , rho &
826 , qrain &
827 , qcloud &
828 , ncloud &
829 , ids, ide, jds, jde, kds, kde &
830 , ims, ime, jms, jme, kms, kme &
831 , ips, ipe, jps, jpe, kps, kpe )
832 ELSE
833 CALL wrf_debug ( DEBUG_LEVEL, 'Icing diagnostics not processed due to unknown MP scheme' )
834 END IF
835 ENDIF ! afwa_icing_opt
836 ELSE
837 CALL wrf_debug ( DEBUG_LEVEL, 'Option icing_diagnostics requires: QC, QR, QNC' )
838 END IF
840 ! Calculate visiblility diagnostics
841 ! ---------------------------------
842 IF ( config_flags % afwa_vis_opt == 1 ) THEN
844 ! Interpolate 20m temperature and RH
845 ! ----------------------------------
846 DO i=ims,ime
847 DO j=jms,jme
848 tv20m(i,j) = -999.
849 rh20m(i,j) = -999.
850 DO k = kps, MIN(kpe,kde-1)
851 IF (tv20m (i,j) .eq. -999. .AND. zagl (i,k,j) .ge. 20.) THEN
853 ! If the lowest model level > 20 m AGL, interpolate using 2-m temperature/RH
854 ! --------------------------------------------------------------------------
855 IF (k .eq. kps) THEN
856 tv20m(i,j) = tv2m(i,j) + &
857 (20. - 2.) / &
858 (zagl(i,k,j) - 2.) * &
859 (grid%t_phy(i,k,j) * (1 + 0.61 * qvapr(i,k,j)) - tv2m(i,j))
860 rh20m(i,j) = rh2m(i,j) + &
861 (20. - 2.) / &
862 (zagl(i,k,j) - 2.) * &
863 (rh(i,k,j) - rh2m(i,j))
864 ELSE
865 tv20m(i,j) = grid%t_phy(i,k-1,j) * (1 + 0.61 * qvapr(i,k-1,j)) + &
866 ((20. - zagl(i,k-1,j)) / &
867 (zagl(i,k,j) - zagl(i,k-1,j))) * &
868 (grid%t_phy(i,k,j) * (1 + 0.61 * qvapr(i,k,j)) - &
869 grid%t_phy(i,k-1,j) * (1 + 0.61 * qvapr(i,k-1,j)))
870 rh20m(i,j) = rh (i,k-1,j) + &
871 ((20. - zagl (i,k-1,j)) / &
872 (zagl (i,k,j) - zagl (i,k-1,j))) * &
873 (rh (i,k,j) - rh (i,k-1,j))
874 ENDIF
875 ENDIF
876 ENDDO
877 ENDDO
878 ENDDO
880 ! Calculate 125 meter winds
881 ! -------------------------
882 DO i=ims,ime
883 DO j=jms,jme
884 wind125m(i,j) = -999.
885 DO k = kps, MIN(kpe,kde-1)
886 IF (wind125m (i,j) .eq. -999. .AND. zagl (i,k,j) .ge. 125.) THEN
888 ! If the lowest model level > 125 m AGL, use level 1 wind
889 ! -------------------------------------------------------
890 IF (k .eq. kps) THEN
891 wind125m(i,j) = uv_wind(u_phy(i,k,j),v_phy(i,k,j))
892 ELSE
893 wind125m(i,j) = uv_wind(u_phy(i,k-1,j),v_phy(i,k-1,j)) + &
894 ((125. - zagl(i,k-1,j)) / &
895 (zagl(i,k,j) - zagl(i,k-1,j))) * &
896 (uv_wind(u_phy(i,k,j),v_phy(i,k,j)) - &
897 uv_wind(u_phy(i,k-1,j),v_phy(i,k-1,j)))
898 ENDIF
899 ENDIF
900 ENDDO
901 ENDDO
902 ENDDO
904 IF ( F_QC .AND. F_QR .AND. F_QI .AND. F_QS .AND. F_QG ) THEN
905 write ( message, * ) 'Calculating visibility'
906 CALL wrf_debug( 100 , message )
907 CALL vis_diagnostics ( qcloud(ims:ime,k_start,jms:jme) &
908 , qrain(ims:ime,k_start,jms:jme) &
909 , qice(ims:ime,k_start,jms:jme) &
910 , qsnow(ims:ime,k_start,jms:jme) &
911 , qgrpl(ims:ime,k_start,jms:jme) &
912 , rho(ims:ime,k_start,jms:jme) &
913 , wind10m &
914 , wind125m &
915 , grid % afwa_pwat &
916 , grid % q2 &
917 , rh2m &
918 , rh20m &
919 , tv2m &
920 , tv20m &
921 , dustc &
922 , grid % afwa_vis &
923 , grid % afwa_vis_dust &
924 , grid % afwa_vis_alpha &
925 , ids, ide, jds, jde, kds, kde &
926 , ims, ime, jms, jme, kms, kme &
927 , ips, ipe, jps, jpe, kps, kpe )
928 ELSE
929 CALL wrf_debug ( DEBUG_LEVEL, 'Option vis_diagnostics requires: QC, QR, QI, QS, QG' )
930 END IF
931 ENDIF
933 IF ( F_QC .AND. F_QI .AND. F_QS ) THEN
934 ! Calculate cloud diagnostics
935 ! ---------------------------
936 IF ( config_flags % afwa_cloud_opt == 1 ) THEN
937 write ( message, * ) 'Calculating cloud'
938 CALL wrf_debug( 100 , message )
939 CALL cloud_diagnostics (qcloud &
940 , qice &
941 , qsnow &
942 , rh_cld &
943 , dz8w &
944 , rho &
945 , grid % z &
946 , grid % ht &
947 , grid % afwa_cloud &
948 , grid % afwa_cloud_ceil &
949 , ids, ide, jds, jde, kds, kde &
950 , ims, ime, jms, jme, kms, kme &
951 , ips, ipe, jps, jpe, kps, kpe )
952 ENDIF
953 ELSE
954 CALL wrf_debug ( DEBUG_LEVEL, 'Option cloud_diagnostics requires: QC, QI, QS' )
955 END IF
957 ENDIF ! is_output_timestep
959 END SUBROUTINE afwa_diagnostics_driver
963 SUBROUTINE severe_wx_diagnostics ( wspd10max &
964 , w_up_max &
965 , w_dn_max &
966 , up_heli_max &
967 , tcoli_max &
968 , midrh_min_old &
969 , midrh_min &
970 , afwa_hail &
971 , cape &
972 , cin &
973 , zlfc &
974 , plfc &
975 , lidx &
976 , llws &
977 , afwa_tornado &
978 , grpl_flx_max &
979 , u10 &
980 , v10 &
981 , w_2 &
982 , uh &
983 , t_phy &
984 , t2 &
985 , z &
986 , ht &
987 , tornado_mask &
988 , tornado_dur &
989 , dt &
990 , pwat &
991 , u_phy &
992 , v_phy &
993 , p &
994 , qi &
995 , qs &
996 , qg &
997 , ng &
998 , qv, qr, qc &
999 , rho &
1000 , dz8w &
1001 , rh &
1002 , do_buoy_calc &
1003 , ims, ime, jms, jme, kms, kme &
1004 , its, ite, jts, jte &
1005 , k_start, k_end &
1006 , j_start, j_end, i_start, i_end )
1009 INTEGER, INTENT(IN) :: its, ite, jts, jte, k_start, k_end &
1010 , ims, ime, jms, jme, kms, kme &
1011 , j_start, j_end, i_start, i_end
1014 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1015 INTENT(IN ) :: p &
1016 , w_2 &
1017 , t_phy &
1018 , u_phy &
1019 , v_phy &
1020 , qi &
1021 , qs &
1022 , qg &
1023 , ng &
1024 , qv, qr, qc &
1025 , rho &
1026 , z &
1027 , dz8w &
1028 , rh
1031 REAL, DIMENSION( ims:ime, jms:jme ), &
1032 INTENT(IN ) :: u10 &
1033 , v10 &
1034 , wspd10max &
1035 , uh &
1036 , t2 &
1037 , ht &
1038 , midrh_min_old &
1039 , up_heli_max &
1040 , llws &
1041 , pwat
1044 REAL, DIMENSION( ims:ime, jms:jme ), &
1045 INTENT(INOUT) :: w_up_max &
1046 , w_dn_max &
1047 , tcoli_max &
1048 , midrh_min &
1049 , afwa_hail &
1050 , afwa_tornado &
1051 , grpl_flx_max &
1052 , tornado_mask &
1053 , tornado_dur
1056 REAL, DIMENSION( ims:ime, jms:jme ), &
1057 INTENT(INOUT) :: cape &
1058 , cin &
1059 , zlfc &
1060 , plfc &
1061 , lidx
1063 REAL, INTENT(IN) :: dt
1064 LOGICAL, INTENT(IN) :: do_buoy_calc
1066 ! Local
1067 ! -----
1068 INTEGER :: i,j,k
1069 INTEGER :: kts,kte
1070 REAL :: zagl, zlfc_msl, melt_term, midrh_term, hail, midrh
1071 REAL :: dum1, dum2, dum3
1072 REAL :: tornado, lfc_term, shr_term, midrh2_term, uh_term
1073 REAL :: wind_vel, p_tot, tcoli, grpl_flx, w_n15, qg_n15
1074 INTEGER :: nz, ostat
1075 REAL, DIMENSION( ims:ime, jms:jme ) :: w_up &
1076 , w_dn &
1077 , tornado_mask_prev &
1078 , tornado_dur_prev
1079 REAL :: time_factor, time_factor_prev
1080 LOGICAL :: is_target_level
1082 ! Calculate midlevel relative humidity minimum
1083 ! --------------------------------------------
1084 DO i=ims,ime
1085 DO j=jms,jme
1086 is_target_level=.false.
1087 DO k=kms,kme
1088 zagl = z(i,k,j) - ht(i,j)
1089 IF ( ( zagl >= 3500. ) .and. &
1090 ( .NOT. is_target_level ) .and. &
1091 ( k .ne. kms ) ) THEN
1092 is_target_level = .true.
1093 midrh = rh(i,k-1,j) + (3500. - (z(i,k-1,j) - ht(i,j))) &
1094 * ((rh(i,k,j) - rh(i,k-1,j))/(z(i,k,j) - z(i,k-1,j)))
1095 IF ( midrh .lt. midrh_min(i,j) ) THEN
1096 midrh_min(i,j) = midrh
1097 ENDIF
1098 ENDIF
1099 ENDDO
1100 ENDDO
1101 ENDDO
1103 ! Calculate the max 10 m wind speed between output times
1104 ! ------------------------------------------------------
1105 !DO j = jts, jte
1106 ! DO i = its, ite
1107 ! wind_vel = uv_wind ( u10(i,j) , v10(i,j) )
1108 ! IF ( wind_vel .GT. wspd10max(i,j) ) THEN
1109 ! wspd10max(i,j) = wind_vel
1110 ! ENDIF
1111 ! ENDDO
1112 !ENDDO
1114 ! Vertical velocity quantities between output times
1115 ! -------------------------------------------------
1116 w_up=0.
1117 w_dn=0.
1118 DO j = jts, jte
1119 DO k = k_start, k_end
1120 DO i = its, ite
1121 p_tot = p(i,k,j) / 100.
1123 ! Check vertical velocity field below 400 mb
1124 IF ( p_tot .GT. 400. .AND. w_2(i,k,j) .GT. w_up(i,j) ) THEN
1125 w_up(i,j) = w_2(i,k,j)
1126 IF ( w_up(i,j) .GT. w_up_max(i,j) ) THEN
1127 w_up_max(i,j) = w_up(i,j)
1128 ENDIF
1129 ENDIF
1130 IF ( p_tot .GT. 400. .AND. w_2(i,k,j) .LT. w_dn(i,j) ) THEN
1131 w_dn(i,j) = w_2(i,k,j)
1132 IF ( w_dn(i,j) .LT. w_dn_max(i,j) ) THEN
1133 w_dn_max(i,j) = w_dn(i,j)
1134 ENDIF
1135 ENDIF
1136 ENDDO
1137 ENDDO
1138 ENDDO
1140 ! Hail diameter in millimeters (Weibull distribution)
1141 ! ---------------------------------------------------
1142 DO j = jts, jte
1143 DO i = its, ite
1144 melt_term=max(t2(i,j)-288.15,0.)
1145 midrh_term=max(2*(min(midrh_min(i,j),midrh_min_old(i,j))-70.),0.)
1146 ! Change exponent to 1.1 to reduce probabilities for large sizes
1147 !hail=max((w_up(i,j)/1.4)**1.25-melt_term-midrh_term,0.)
1148 hail=max((w_up(i,j)/1.4)**1.1-melt_term-midrh_term,0.)
1149 hail=hail*((uh(i,j)/100)+0.25)
1150 IF ( hail .gt. afwa_hail(i,j) ) THEN
1151 afwa_hail(i,j)=hail
1152 ENDIF
1153 ENDDO
1154 ENDDO
1156 ! Lightning (total column-integrated cloud ice)
1157 ! Note this formula is basically stolen from the VIL calculation.
1158 ! ---------------------------------------------------------------
1159 DO j = jts, jte
1160 DO i = its, ite
1161 tcoli=0.
1162 DO k = k_start, k_end
1163 tcoli = tcoli + &
1164 (qi (i,k,j) + &
1165 qs (i,k,j) + &
1166 qg (i,k,j)) &
1167 *dz8w (i,k,j) * rho(i,k,j)
1168 ENDDO
1169 IF ( tcoli .GT. tcoli_max(i,j) ) THEN
1170 tcoli_max(i,j) = tcoli
1171 ENDIF
1172 ENDDO
1173 ENDDO
1175 ! Lighting (pseudo graupel flux calculation)
1176 ! Model graupel mixing ration (g/kg) times w (m/s) at the -15C level
1177 ! Values should range from around 25 in marginal lightning situations,
1178 ! to over 400 in situations with very frequent lightning.
1179 ! --------------------------------------------------------------------
1180 DO j = jts, jte
1181 DO i = its, ite
1182 grpl_flx=0
1183 w_n15=-999.
1184 DO k = k_start, k_end
1186 ! Interpolate w and qg to -15C level and calculate graupel flux
1187 ! as simply graupel x vertical velocity at -15C
1188 ! -------------------------------------------------------------
1189 IF ( t_phy (i,k,j) .LE. 258.15 .AND. w_n15 .EQ. -999. .AND. &
1190 k .GT. k_start .AND. qg (i,k,j) .GT. 1.E-20 ) THEN
1191 w_n15 = w_2 (i,k,j)
1192 qg_n15 = 1000. * qg (i,k,j) ! g/kg
1193 grpl_flx = qg_n15 * w_n15
1194 ENDIF
1195 ENDDO
1196 IF ( grpl_flx .GT. grpl_flx_max(i,j) ) THEN
1197 grpl_flx_max(i,j) = grpl_flx
1198 ENDIF
1199 ENDDO
1200 ENDDO
1202 ! Update buoyancy parameters.
1203 ! ---------------------------
1204 IF ( do_buoy_calc ) THEN
1205 nz = k_end - k_start + 1
1206 DO j = jts, jte
1207 DO i = its, ite
1208 ostat = Buoyancy ( nz &
1209 , t_phy(i,kms:kme ,j) &
1210 , rh(i,kms:kme ,j) &
1211 , p(i,kms:kme ,j) &
1212 , z(i,kms:kme ,j) &
1213 , 1 &
1214 , cape(i,j) &
1215 , cin(i,j) &
1216 , zlfc_msl &
1217 , plfc(i,j) &
1218 , lidx(i,j) &
1219 , 3 ) !Surface
1221 ! Add 273.15 to lifted index per some standard I don't understand
1222 ! ---------------------------------------------------------------
1223 IF ( lidx ( i, j ) .ne. 999. ) lidx ( i, j ) = lidx ( i, j ) + 273.15
1225 ! Subtract terrain height to convert ZLFC from MSL to AGL
1226 ! -------------------------------------------------------
1227 IF ( zlfc_msl .ge. 0. ) THEN
1228 zlfc ( i, j ) = zlfc_msl - ht ( i, j )
1229 ELSE
1230 zlfc( i, j ) = -1.
1231 ENDIF
1233 ENDDO
1234 ENDDO
1235 ENDIF
1237 ! Maximum tornado wind speed in ms-1.
1238 ! First, save off old tornado mask and duration arrays
1239 ! ----------------------------------------------------
1240 tornado_dur_prev(:,:)=tornado_dur(:,:)
1241 tornado_mask_prev(:,:)=tornado_mask(:,:)
1243 ! Initialize all tornado variables
1244 ! --------------------------------
1245 tornado_mask(:,:)=0.
1246 tornado_dur(:,:)=0.
1248 DO j = j_start, j_end
1249 DO i = i_start, i_end
1250 tornado = 0.
1252 ! Current tornado algorithm
1253 ! -------------------------
1254 IF ( zlfc(i,j) .ge. 0. ) THEN
1255 uh_term = min(max((uh(i,j) - 25.) / 50., 0.), 1.)
1256 shr_term = min(max((llws(i,j) - 2.) / 10., 0.), 1.)
1257 lfc_term = min(max((3000. - zlfc(i,j)) / 1500., 0.), 1.)
1258 midrh2_term = min(max((90. - &
1259 min(midrh_min(i,j),midrh_min_old(i,j))) / 30., 0.), 1.)
1260 tornado = 50. * uh_term * shr_term * lfc_term * midrh2_term
1261 ENDIF
1263 ! Does tornado algorithm indicate all possible ingredients
1264 ! for a minimum tornado, if so turn on mask
1265 ! --------------------------------------------------------
1266 IF (tornado .GT. 0.) THEN
1267 tornado_mask(i,j) = 1.
1268 ENDIF
1270 ! Update the duration of this tornado if there was previously
1271 ! a tornado mask at this or an adjacent point
1272 ! -----------------------------------------------------------
1273 IF ( ( tornado_mask(i,j) .GT. 0.5 ) .OR. &
1274 ( MAXVAL(tornado_mask_prev(i-1:i+1,j-1:j+1)) .GT. 0.5 ) ) THEN
1275 tornado_dur(i,j) = MAXVAL(tornado_dur_prev(i-1:i+1,j-1:j+1)) + dt
1276 ENDIF
1278 ! Ramp up value of tornado beta value linearly in time until &
1279 ! it has been sustained 5 minutes
1280 ! ------------------------------------------------------------
1281 time_factor = MIN(tornado_dur(i,j)/900.,1.)
1282 tornado = tornado * time_factor
1284 ! OPTIONAL
1285 ! Adjust previous max tornado wind upward to account for longer
1286 ! duration - if after 5 minutes, no adjustment made as
1287 ! time_factor/time_factor_prev=1
1288 ! -------------------------------------------------------------
1289 !time_factor_prev = MIN((tornado_dur(i,j) - dt)/900.,1.)
1290 !IF ( ( time_factor_prev .GT. 0. ) .AND. &
1291 ! ( time_factor_prev .LT. 1. ) ) THEN
1292 ! afwa_tornado(i,j) = afwa_tornado(i,j) * time_factor/time_factor_prev
1293 !ENDIF
1295 ! Now that we are comparing apples to apples, see if current tornado
1296 ! wind is higher than previous highest value for this point
1297 ! ------------------------------------------------------------------
1298 IF ( tornado .GT. afwa_tornado(i,j) ) THEN
1299 afwa_tornado(i,j) = tornado
1300 ENDIF
1301 ENDDO
1302 ENDDO
1304 END SUBROUTINE severe_wx_diagnostics
1308 SUBROUTINE vert_int_liquid_diagnostics ( vil &
1309 , radarvil &
1310 , t_phy &
1311 , qrain &
1312 , qsnow &
1313 , qgrpl &
1314 , z_e &
1315 , dz8w &
1316 , rho &
1317 , ids, ide, jds, jde, kds, kde &
1318 , ims, ime, jms, jme, kms, kme &
1319 , ips, ipe, jps, jpe, kps, kpe )
1321 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1322 ims, ime, jms, jme, kms, kme, &
1323 ips, ipe, jps, jpe, kps, kpe
1325 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1326 INTENT(IN) :: rho &
1327 , qrain &
1328 , qsnow &
1329 , qgrpl &
1330 , t_phy &
1331 , z_e &
1332 , dz8w
1334 REAL, DIMENSION( ims:ime, jms:jme ), &
1335 INTENT(INOUT) :: vil &
1336 , radarvil
1338 ! Local
1339 ! -----
1340 INTEGER :: i,j,k,ktime
1342 ! Calculate vertically integrated liquid water (though its mostly not
1343 ! "liquid" now is it?) [kg/m^2]
1344 ! -------------------------------------------------------------------
1345 DO i = ips, MIN(ipe,ide-1)
1346 DO j = jps, MIN(jpe,jde-1)
1347 vil (i,j) = 0.0
1348 DO k = kps, MIN(kpe,kde-1)
1349 vil (i,j) = vil (i,j) + &
1350 (qrain (i,k,j) + &
1351 qsnow (i,k,j) + &
1352 qgrpl (i,k,j)) &
1353 *dz8w (i,k,j) * rho(i,k,j)
1354 ENDDO
1355 ENDDO
1356 ENDDO
1358 ! Diagnose "radar-derived VIL" from equivalent radar reflectivity
1359 ! radarVIL = (integral of LW*dz )/1000.0 (in kg/m^2)
1360 ! LW = 0.00344 * z_e** (4/7) in g/m^3
1361 ! ---------------------------------------------------------------
1362 DO i = ips, MIN(ipe,ide-1)
1363 DO j = jps, MIN(jpe,jde-1)
1364 radarvil (i,j) = 0.0
1365 DO k = kps, MIN(kpe,kde-1)
1366 radarvil (i,j) = radarvil (i,j) + &
1367 0.00344 * z_e(i,k,j)**0.57143 &
1368 *dz8w (i,k,j)/1000.0
1369 END DO
1370 END DO
1371 END DO
1373 END SUBROUTINE vert_int_liquid_diagnostics
1377 SUBROUTINE icing_diagnostics ( icing_opt &
1378 , fzlev &
1379 , icing_lg &
1380 , icing_sm &
1381 , qicing_lg_max &
1382 , qicing_sm_max &
1383 , qicing_lg &
1384 , qicing_sm &
1385 , icingtop &
1386 , icingbot &
1387 , t_phy &
1388 , z &
1389 , dz8w &
1390 , rho &
1391 , qrain &
1392 , qcloud &
1393 , ncloud &
1394 , ids, ide, jds, jde, kds, kde &
1395 , ims, ime, jms, jme, kms, kme &
1396 , ips, ipe, jps, jpe, kps, kpe )
1398 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1399 ims, ime, jms, jme, kms, kme, &
1400 ips, ipe, jps, jpe, kps, kpe
1402 INTEGER, INTENT(IN) :: icing_opt
1404 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1405 INTENT(IN) :: z &
1406 , qrain &
1407 , qcloud &
1408 , ncloud &
1409 , rho &
1410 , dz8w &
1411 , t_phy
1413 REAL, DIMENSION( ims:ime, jms:jme ), &
1414 INTENT( OUT) :: fzlev &
1415 , icing_lg &
1416 , icing_sm &
1417 , qicing_lg_max &
1418 , qicing_sm_max &
1419 , icingtop &
1420 , icingbot
1422 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1423 INTENT( OUT) :: qicing_lg &
1424 , qicing_sm
1427 ! Local
1428 ! -----
1429 INTEGER :: i,j,k,ktime,ktop,kbot
1430 REAL :: qcfrac_lg, qcfrac_sm, qc, qr, small, all
1432 ! Initializations
1433 ! ---------------
1434 fzlev (ips:ipe,jps:jpe) = -999. ! Arbitrary unset/initial value
1435 icingtop (ips:ipe,jps:jpe) = -999. ! Arbitrary unset/initial value
1436 icingbot (ips:ipe,jps:jpe) = -999. ! Arbitrary unset/initial value
1437 icing_lg (ips:ipe,jps:jpe) = 0.
1438 icing_sm (ips:ipe,jps:jpe) = 0.
1439 qicing_lg_max (ips:ipe,jps:jpe) = 0.
1440 qicing_sm_max (ips:ipe,jps:jpe) = 0.
1441 qicing_sm(ips:ipe,kps:kpe,jps:jpe)=0.
1442 qicing_lg(ips:ipe,kps:kpe,jps:jpe)=0.
1444 ! Loop through i and j
1445 ! --------------------
1446 DO i = ips, MIN(ipe,ide-1)
1447 DO j = jps, MIN(jpe,jde-1)
1449 ! Go up the column and look for sub freezing temperatures
1450 ! -------------------------------------------------------
1451 ktop=-1
1452 kbot=-1
1453 DO k = kps, MIN(kpe,kde-1)
1454 IF (t_phy(i,k,j) .lt. 273.15) THEN
1456 ! Any cloud water we find will be supercooled.
1457 ! Based on microphysics scheme, determine the fraction of
1458 ! large (>50 um) supercooled cloud water drops.
1459 ! Source: Becky Selin, 16WS
1460 ! -------------------------------------------------------
1461 qc = qcloud (i,k,j)
1462 qr = qrain (i,k,j)
1463 nc = ncloud(i,k,j)
1464 den = rho(i,k,j)
1465 qcfrac_lg = 0.
1466 qcfrac_sm = 0.
1468 ! Eta (Ferrier)
1469 ! -------------
1470 IF (icing_opt .eq. 2) THEN
1471 IF (qc .lt. 2.5E-4) THEN
1472 qcfrac_lg = 395000. * qc**2. + 102.9 * qc
1473 ELSEIF (qc .lt. 1.4E-3) THEN
1474 qcfrac_lg = 276.1 * qc - 0.01861
1475 ELSE
1476 qcfrac_lg = 0.3 * log(641.789 * qc) + 0.4
1477 ENDIF
1479 ! Thompson
1480 ! --------
1481 ! RAS Per James McCormick's stats, more large supercooled
1482 ! drops are needed from the Thompson members. Changing
1483 ! calculation to be like WSM5/6 members.
1484 !ELSEIF (icing_opt .eq. 3) THEN
1485 ! IF (qc .lt. 1.0E-3) THEN
1486 ! qcfrac_lg = 2205.0 * qc**2. + 3.232 * qc
1487 ! ELSEIF (qc .lt. 3.0E-3) THEN
1488 ! qcfrac_lg = 24.1 * qc - 0.01866
1489 ! ELSE
1490 ! qcfrac_lg = 0.127063 * log(550.0 * qc) - 0.01
1491 ! ENDIF
1493 ! Thompson or WSM5/6
1494 ! ------------------
1495 !ELSEIF (icing_opt .eq. 4) THEN
1496 ELSEIF ((icing_opt .eq. 3) .OR. (icing_opt .eq. 4)) THEN
1497 IF (qc .lt. 5.E-4) THEN
1498 qcfrac_lg = 50420.0 * qc**2. + 29.39 * qc
1499 ELSEIF (qc .lt. 1.4E-3) THEN
1500 qcfrac_lg = 97.65 * qc - 0.02152
1501 ELSE
1502 qcfrac_lg = 0.2 * log(646.908 * qc) + 0.135
1503 ENDIF
1505 ! Morrison 2-moment, constant CCN
1506 ! -------------------------------
1507 ELSEIF (icing_opt .eq. 5) THEN
1508 IF (qc .lt. 1.4E-3) THEN
1509 qcfrac_lg = 28000. * qc**2. + 0.1 * qc
1510 ELSEIF (qc .lt. 2.6E-3) THEN
1511 qcfrac_lg = 112.351 * qc - 0.102272
1512 ELSE
1513 qcfrac_lg = 0.3 * log(654.92 * qc) * 0.301607
1514 ENDIF
1516 ! WDM6 or Morrison 2-moment w/ prognostic CCN
1517 ! -------------------------------------------
1518 ELSEIF ((icing_opt .eq. 6) .OR. (icing_opt .eq. 7)) THEN
1519 IF ((qc .gt. 1.0E-12) .and. (nc .gt. 1.0E-12)) THEN
1520 small = -nc * exp(-nc*3141.59265*(5.E-5)**3./(6000.*den*qc))+nc
1521 all = -nc * exp(-nc*3141.59265*(2.)**3./(6000.*den*qc))+nc
1522 qcfrac_lg = 1. - (small / all)
1523 ELSE
1524 qcfac_lg = 0.
1525 ENDIF
1526 ENDIF
1527 qcfrac_lg = max(qcfrac_lg, 0.)
1529 ! Small (<50 um) supercooled cloud water drop fraction (1 - large).
1530 ! -----------------------------------------------------------------
1531 IF (icing_opt .ne. 0 ) THEN
1532 qcfrac_sm = 1 - qcfrac_lg
1533 ENDIF
1535 ! Supercooled drop mixing ratio
1536 ! -----------------------------
1537 qicing_lg (i,k,j) = max(qr + qcfrac_lg * qc, 0.)
1538 qicing_sm (i,k,j) = max(qcfrac_sm * qc, 0.)
1540 ! Column integrated icing
1541 ! -----------------------
1542 icing_lg (i,j) = icing_lg (i,j) + qicing_lg (i,k,j) &
1543 * dz8w (i,k,j) * rho(i,k,j)
1544 icing_sm (i,j) = icing_sm (i,j) + qicing_sm (i,k,j) &
1545 * dz8w (i,k,j) * rho(i,k,j)
1547 ! Column maximum supercooled drop mixing ratio
1548 ! --------------------------------------------
1549 IF ( qicing_lg(i,k,j) .gt. qicing_lg_max(i,j) ) THEN
1550 qicing_lg_max (i,j) = qicing_lg(i,k,j)
1551 ENDIF
1552 IF ( qicing_sm(i,k,j) .gt. qicing_sm_max(i,j) ) THEN
1553 qicing_sm_max (i,j) = qicing_sm(i,k,j)
1554 ENDIF
1556 ! Freezing level calculation
1557 ! --------------------------
1558 IF (fzlev (i,j) .eq. -999.) THEN ! At freezing level
1559 IF (k .ne. kps) THEN ! If not at surface, interpolate.
1560 fzlev (i,j) = z (i,k-1,j) + &
1561 ((273.15 - t_phy (i,k-1,j)) &
1562 /(t_phy (i,k,j) - t_phy (i,k-1,j))) &
1563 *(z (i,k,j) - z (i,k-1,j))
1564 ELSE ! If at surface, use first level.
1565 fzlev(i,j) = z (i,k,j)
1566 ENDIF
1567 ENDIF
1569 ! Icing layer top and bottom indices (where icing > some arbitrary
1570 ! small value). Set bottom index of icing layer to current k index
1571 ! if not yet set. Set top index of icing layer to current k index.
1572 ! ----------------------------------------------------------------
1573 IF ((qicing_lg (i,k,j) + qicing_sm (i,k,j)) .ge. 1.E-5) THEN
1574 IF (kbot .eq. -1) kbot = k
1575 ktop=k
1576 ENDIF
1577 ENDIF
1578 END DO
1580 ! Interpolate bottom of icing layer from kbot (bottom index of icing
1581 ! layer). Icing bottom should not go below freezing level.
1582 ! ------------------------------------------------------------------
1583 IF (kbot .ne. -1) THEN
1584 IF (kbot .ne. kps) THEN ! If not at surface, interpolate
1585 icingbot (i,j) = z (i,kbot-1,j) + ((1.E-5 - &
1586 (qicing_lg (i,kbot-1,j) + qicing_sm (i,kbot-1,j))) &
1587 / ((qicing_lg (i,kbot,j) + qicing_sm (i,kbot,j)) &
1588 - (qicing_lg (i,kbot-1,j) + qicing_sm (i,kbot-1,j)))) &
1589 * (z (i,kbot,j) - z (i,kbot-1,j))
1590 icingbot (i,j) = MAX(icingbot (i,j), fzlev (i,j))
1591 ELSE ! If at surface use first level.
1592 icingbot (i,j) = z(i,kbot,j)
1593 ENDIF
1594 ENDIF
1596 ! Interpolate top of icing layer from ktop (top index of icing layer).
1597 ! Icing top should not go below icing bottom (obviously).
1598 ! --------------------------------------------------------------------
1599 IF (ktop .ne. -1 .and. ktop .ne. kpe) THEN ! If not undefined or model top
1600 icingtop (i,j) = z (i,ktop,j) + ((1.E-5 - &
1601 (qicing_lg (i,ktop,j) + qicing_sm (i,ktop,j))) &
1602 / ((qicing_lg (i,ktop+1,j) + qicing_sm (i,ktop+1,j)) &
1603 - (qicing_lg (i,ktop,j) + qicing_sm (i,ktop,j)))) &
1604 * (z (i,ktop+1,j) - z (i,ktop,j))
1605 icingtop (i,j) = MAX(icingtop (i,j), icingbot (i,j))
1606 ENDIF
1607 END DO
1608 END DO
1610 END SUBROUTINE icing_diagnostics
1614 SUBROUTINE radar_diagnostics ( refd &
1615 , refd_com &
1616 ! , refd_max &
1617 , echotop &
1618 , z &
1619 , z_e &
1620 , ids, ide, jds, jde, kds, kde &
1621 , ims, ime, jms, jme, kms, kme &
1622 , ips, ipe, jps, jpe, kps, kpe )
1624 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1625 ims, ime, jms, jme, kms, kme, &
1626 ips, ipe, jps, jpe, kps, kpe
1628 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1629 INTENT(IN) :: z &
1630 , z_e
1632 REAL, DIMENSION( ims:ime, jms:jme ), &
1633 INTENT(INOUT) :: refd &
1634 , refd_com &
1635 ! , refd_max &
1636 , echotop
1638 ! Local
1639 ! -----
1640 INTEGER :: i,j,k,ktime
1642 DO j = jps, MIN(jpe,jde-1)
1643 DO i = ips, MIN(ipe,ide-1)
1644 ktop = -1 ! Undefined
1645 echotop (i,j) = 0.
1646 refd_com (i,j) = 0.
1647 refd (i,j) = 0.
1648 DO k = kps, MIN(kpe,kde-1)
1649 IF (z_e(i,k,j) .gt. 1.e-20) THEN
1651 ! Reflectivity (first level)
1652 ! --------------------------
1653 IF (k == kps) refd(i,j) = MAX(10.0 * log10(z_e(i,k,j)),0.)
1655 ! ! Max reflectivity over the output interval
1656 ! ! -----------------------------------------
1657 ! IF (refd(i,j) .gt. refd_max(i,j)) refd_max(i,j) = refd(i,j)
1659 ! Composite reflectivity calc (max reflectivity in the column)
1660 ! ------------------------------------------------------------
1661 IF (10.0 * log10(z_e(i,k,j)) .gt. refd_com(i,j)) THEN
1662 refd_com(i,j) = 10.0 * log10(z_e(i,k,j))
1663 ENDIF
1664 ENDIF
1666 ! Echo top - the highest level w/ dBZ > 18 (z_e > 63.0957)
1667 ! --------------------------------------------------------
1668 IF ( z_e(i,k,j) .gt. 63.0957) THEN
1669 ktop = k
1670 ENDIF
1671 END DO
1672 IF ( ktop .ne. -1 ) THEN ! Interpolate to echo top height (GAC)
1673 echotop (i,j) = z (i,ktop,j) + &
1674 ((63.0957 - z_e (i,ktop,j)) &
1675 /(z_e (i,ktop+1,j) - z_e (i,ktop,j))) &
1676 *(z (i,ktop+1,j) - z (i,ktop,j))
1677 ENDIF
1678 END DO
1679 END DO
1681 END SUBROUTINE radar_diagnostics
1685 SUBROUTINE wrf_dbzcalc( rho &
1686 , t_phy &
1687 , qr &
1688 , qs &
1689 , qg &
1690 , z_e &
1691 , ids, ide, jds, jde, kds, kde &
1692 , ims, ime, jms, jme, kms, kme &
1693 , ips, ipe, jps, jpe, kps, kpe )
1695 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1696 ims, ime, jms, jme, kms, kme, &
1697 ips, ipe, jps, jpe, kps, kpe
1699 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1700 INTENT(IN ) :: rho &
1701 , t_phy &
1702 , qr &
1703 , qs &
1704 , qg
1706 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1707 INTENT( OUT) :: z_e
1709 REAL :: factor_r, factor_s, factor_g, factorb_s, factorb_g, ronv, sonv, gonv
1710 REAL :: temp_c, rhoair, qgr, qra, qsn
1711 INTEGER :: i, j, k
1713 INTEGER, PARAMETER :: iBrightBand = 1
1714 REAL, PARAMETER :: T_0 = 273.15
1715 REAL, PARAMETER :: PI = 3.1415926536
1716 REAL, PARAMETER :: rgas=287.04, gamma_seven = 720.0, alpha2 = 0.224
1718 ! Densities of rain, snow, graupel, and cloud ice.
1719 ! ------------------------------------------------
1720 REAL, PARAMETER :: rho_w = 1000.0, rho_r = 1000.0, rho_s = 100.0
1721 REAL, PARAMETER :: rho_g = 400.0, rho_i = 890.0
1722 REAL, PARAMETER :: ron=8.e6, son=2.e7, gon=5.e7, r1=1.e-15
1723 REAL, PARAMETER :: ron_min = 8.e6, ron2=1.e10
1724 REAL, PARAMETER :: ron_qr0 = 0.0001, ron_delqr0 = 0.25*ron_qr0
1725 REAL, PARAMETER :: ron_const1r = (ron2-ron_min)*0.5
1726 REAL, PARAMETER :: ron_const2r = (ron2+ron_min)*0.5
1728 ! Constant intercepts
1729 ! -------------------
1730 ronv = 8.e6 ! m^-4
1731 sonv = 2.e7 ! m^-4
1732 gonv = 4.e6 ! m^-4
1734 factor_r = gamma_seven * 1.e18 * (1./(pi*rho_r))**1.75
1735 factor_s = gamma_seven * 1.e18 * (1./(pi*rho_s))**1.75 &
1736 * (rho_s/rho_w)**2 * alpha2
1737 factor_g = gamma_seven * 1.e18 * (1./(pi*rho_g))**1.75 &
1738 * (rho_g/rho_w)**2 * alpha2
1740 ! For each grid point
1741 ! -------------------
1742 DO j = jps, jpe
1743 DO k = kps, kpe
1744 DO i = ips, ipe
1746 factorb_s = factor_s
1747 factorb_g = factor_g
1749 ! In this case snow or graupel particle scatters like liquid
1750 ! water because it is assumed to have a liquid skin
1751 ! ----------------------------------------------------------
1752 IF( iBrightBand == 1 ) THEN
1753 IF (t_phy(i,k,j) > T_0) THEN
1754 factorb_s = factor_s /alpha2
1755 factorb_g = factor_g /alpha2
1756 ENDIF
1757 ENDIF
1759 ! Calculate variable intercept parameters
1760 ! ---------------------------------------
1761 temp_c = amin1(-0.001, t_phy(i,k,j)- T_0)
1762 sonv = amin1(2.0e8, 2.0e6*exp(-0.12*temp_c))
1763 gonv = gon
1764 qgr = QG(i,k,j)
1765 qra = QR(i,k,j)
1766 qsn = QS(i,k,j)
1767 IF (qgr.gt.r1) THEN
1768 gonv = 2.38*(pi*rho_g/(rho(i,k,j)*qgr))**0.92
1769 gonv = max(1.e4, min(gonv,gon))
1770 ENDIF
1771 ronv = ron2
1772 IF (qra.gt. r1) THEN
1773 ronv = ron_const1r*tanh((ron_qr0-qra)/ron_delqr0) + ron_const2r
1774 ENDIF
1776 IF (qra < 0.0 ) qra = 0.0
1777 IF (qsn < 0.0 ) qsn = 0.0
1778 IF (qgr < 0.0 ) qgr = 0.0
1779 z_e(i,k,j) = factor_r * (rho(i,k,j) * qra)**1.75 / ronv**.75 + &
1780 factorb_s * (rho(i,k,j) * qsn)**1.75 / sonv**.75 + &
1781 factorb_g * (rho(i,k,j) * qgr)**1.75 / gonv**.75
1783 IF ( z_e(i,k,j) < 0.0 ) z_e(i,k,j) = 0.0
1785 END DO
1786 END DO
1787 END DO
1789 END SUBROUTINE wrf_dbzcalc
1793 SUBROUTINE precip_type_diagnostics ( t_phy &
1794 , t2 &
1795 , rh &
1796 , z &
1797 , dz8w &
1798 , ht &
1799 , precip &
1800 , swdown &
1801 , rain &
1802 , snow &
1803 , ice &
1804 , frz_rain &
1805 , snowfall &
1806 , ccn_tmp &
1807 , total_melt &
1808 , ids, ide, jds, jde, kds, kde &
1809 , ims, ime, jms, jme, kms, kme &
1810 , ips, ipe, jps, jpe, kps, kpe )
1812 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
1813 ims, ime, jms, jme, kms, kme, &
1814 ips, ipe, jps, jpe, kps, kpe
1816 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
1817 INTENT(IN ) :: t_phy &
1818 , rh &
1819 , z &
1820 , dz8w
1821 REAL, DIMENSION( ims:ime, jms:jme ), &
1822 INTENT(IN ) :: t2 &
1823 , ht &
1824 , precip &
1825 , swdown
1826 REAL, DIMENSION( ims:ime, jms:jme ), &
1827 INTENT(INOUT) :: snowfall &
1828 , rain &
1829 , frz_rain &
1830 , snow &
1831 , ice
1832 REAL, INTENT(IN) :: ccn_tmp
1833 REAL, INTENT(IN) :: total_melt
1835 ! Local
1836 ! -----
1837 REAL, DIMENSION( ims:ime, jms:jme ) :: &
1838 melt &
1839 , mod_2m_tmp &
1840 , cloud_top_tmp &
1841 , maxtmp
1843 INTEGER, DIMENSION( ims:ime, jms:jme ) :: &
1844 cloud_top_k_index &
1845 , precip_type
1847 LOGICAL, DIMENSION (ims:ime, jms:jme ) :: &
1848 saturation
1850 REAL, PARAMETER :: snow_ratio=5.0
1852 ! Loop through all points
1853 ! Search vertically twice--first to find the cloud top temperature and the
1854 ! maximum temperature. Second, determine if any melting or re-freezing will
1855 ! occur to make ice pellets or freezing rain
1856 ! -------------------------------------------------------------------------
1857 DO i=ips,ipe
1858 DO j=jps,jpe
1860 saturation(i,j)=.false.
1861 melt(i,j)=0.0
1862 precip_type(i,j)=0
1864 ! Modify surface temperature for solar insolation (W/m2)
1865 ! Set max temperature in the atmopshere
1866 ! ------------------------------------------------------
1867 mod_2m_tmp(i,j)=t2(i,j)+(swdown(i,j)/100.0)
1868 maxtmp(i,j)=mod_2m_tmp(i,j)
1870 ! Only look at points that have precip and are not warm at the surface
1871 ! --------------------------------------------------------------------
1872 IF (precip(i,j) .gt. 0.0) THEN
1873 !IF (mod_2m_tmp(i,j) .gt. 277.15) THEN
1874 IF (mod_2m_tmp(i,j) .gt. 275.15) THEN
1875 precip_type(i,j)=1 ! Rain
1876 ELSE
1878 ! Check sounding from top for saturation (RH-water gt 80%)--this is
1879 ! the cloud top. Erase saturation if RH lt 70% (spurious moist layer
1880 ! aloft)
1881 ! ------------------------------------------------------------------
1882 cloud_top_k_index(i,j)=kpe
1883 DO k=kpe,kps,-1
1884 IF ((z(i,k,j)-ht(i,j)) .gt. 0.0) THEN
1885 IF (t_phy(i,k,j) .gt. maxtmp(i,j)) THEN
1886 maxtmp(i,j)=t_phy(i,k,j)
1887 ENDIF
1888 IF (rh(i,k,j) .gt. 80 .and. saturation(i,j) .eqv. .false.) THEN
1889 cloud_top_tmp(i,j)=t_phy(i,k,j)
1890 cloud_top_k_index(i,j)=k
1891 saturation(i,j)=.true.
1892 precip_type(i,j)=2 ! Snow
1893 ENDIF
1894 IF (rh(i,k,j) .le. 70 .and. saturation(i,j) .eqv. .true.) THEN
1895 saturation(i,j)=.false.
1896 ENDIF
1897 ENDIF
1898 ENDDO
1900 ! Perform simple check to assign types with no melting layer
1901 ! shenanigans going on
1902 ! ----------------------------------------------------------------
1903 IF (cloud_top_tmp(i,j) .le. ccn_tmp .and. &
1904 maxtmp(i,j) .le. 273.15) THEN
1905 precip_type(i,j)=2 ! Snow
1906 ENDIF
1908 ! ELSE, have to go through the profile again to see if snow melts,
1909 ! and if anything re-freezes
1910 ! ----------------------------------------------------------------
1911 DO k=cloud_top_k_index(i,j),kps,-1
1912 IF ((z(i,k,j)-ht(i,j)) .gt. 0.0) THEN
1914 ! Condition 0--assign falling rain when we get to the
1915 ! supercooled temperature if too warm
1916 ! ---------------------------------------------------
1917 IF (cloud_top_tmp(i,j) .eq. t_phy(i,k,j) .and. &
1918 cloud_top_tmp(i,j) .gt. ccn_tmp) THEN
1919 precip_type(i,j)=1 ! Rain
1920 ENDIF
1922 ! Condition 1--falling frozen precip that will start to melt
1923 ! Add up melting energy over warm layers--if enough, turn to
1924 ! liquid
1925 ! ----------------------------------------------------------
1926 IF ((precip_type(i,j) .eq. 2 .or. precip_type(i,j) .eq. 3) .and. &
1927 t_phy(i,k,j) .gt. 273.15) THEN
1928 melt(i,j)=melt(i,j)+9.8*(((t_phy(i,k,j)-273.15)/273.15)* &
1929 (dz8w(i,k,j)))
1930 IF (melt(i,j) .gt. total_melt) THEN
1931 precip_type(i,j)=1 ! Rain
1932 melt(i,j)=0.0 ! Reset melting energy in case it re-freezes
1933 ENDIF
1934 ENDIF
1936 ! Condition 2--falling partially melted precip encounters
1937 ! sub-freezing air. Snow will be converted to ice pellets if
1938 ! at least 1/4 of it melted. Instantaneous freeze-up, simplistic
1939 ! --------------------------------------------------------------
1940 IF (t_phy(i,k,j) .le. 273.15 .and. &
1941 melt(i,j) .gt. total_melt/4.0 .and. &
1942 (precip_type(i,j) .eq. 2 .or. precip_type(i,j) .eq. 3)) THEN
1943 precip_type(i,j)=3 ! Ice
1944 melt(i,j)=0.0
1945 ENDIF
1947 ! Condition 3--falling liquid that will re-freeze--must reach
1948 ! nucleation temperature
1949 ! -----------------------------------------------------------
1950 IF (precip_type(i,j) .eq. 1) THEN
1951 IF (t_phy(i,k,j) .le. ccn_tmp) THEN
1952 precip_type(i,j)=3 ! Ice
1953 ENDIF
1954 ENDIF
1955 ENDIF ! End if (z-ht)>0
1956 ENDDO ! End do k=kpe,kps,-1
1957 ENDIF ! End if mod_2m_tmp>273.15
1959 ! Accumulate precip according to precip_type
1960 ! ------------------------------------------
1961 IF (precip_type(i,j) .eq. 3) THEN
1962 ice(i,j)=ice(i,j)+precip(i,j)
1963 ENDIF
1964 IF (precip_type(i,j) .eq. 2) THEN
1965 snow(i,j)=snow(i,j)+precip(i,j)
1966 snowfall(i,j)=snowfall(i,j)+snow_ratio*precip(i,j) &
1967 *(5.-mod_2m_tmp(i,j)+273.15)**0.4
1968 ENDIF
1969 IF (precip_type(i,j) .eq. 1) THEN
1970 IF (mod_2m_tmp(i,j) .gt. 273.15) THEN
1971 rain(i,j)=rain(i,j)+precip(i,j)
1972 ELSE
1973 frz_rain(i,j)=frz_rain(i,j)+precip(i,j)
1974 ENDIF
1975 ENDIF
1977 ENDIF ! End if precip>0
1979 ENDDO ! End do j=jps,jpe
1980 ENDDO ! End do i=ips,ipe
1982 END SUBROUTINE precip_type_diagnostics
1986 SUBROUTINE vis_diagnostics ( qcloud &
1987 , qrain &
1988 , qice &
1989 , qsnow &
1990 , qgrpl &
1991 , rho &
1992 , wind10m &
1993 , wind125m &
1994 , pwater &
1995 , q2m &
1996 , rh2m &
1997 , rh20m &
1998 , tv2m &
1999 , tv20m &
2000 , dustc &
2001 , vis &
2002 , vis_dust &
2003 , vis_alpha &
2004 , ids, ide, jds, jde, kds, kde &
2005 , ims, ime, jms, jme, kms, kme &
2006 , ips, ipe, jps, jpe, kps, kpe )
2008 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2009 ims, ime, jms, jme, kms, kme, &
2010 ips, ipe, jps, jpe, kps, kpe
2012 INTEGER, PARAMETER :: ndust=5
2014 REAL, DIMENSION( ims:ime, jms:jme ), &
2015 INTENT(IN ) :: qcloud &
2016 , qrain &
2017 , qice &
2018 , qsnow &
2019 , qgrpl &
2020 , rho &
2021 , wind10m &
2022 , wind125m &
2023 , pwater &
2024 , rh2m &
2025 , q2m &
2026 , rh20m &
2027 , tv2m &
2028 , tv20m
2029 REAL, DIMENSION( ims:ime, jms:jme, ndust ), &
2030 INTENT(IN ) :: dustc
2031 REAL, DIMENSION( ims:ime, jms:jme ), &
2032 INTENT( OUT) :: vis &
2033 , vis_dust &
2034 , vis_alpha
2036 ! Local
2037 ! -----
2038 INTEGER :: i,j,k,d
2039 REAL, PARAMETER :: visfactor=3.912
2040 REAL, DIMENSION (ndust) :: dustfact
2041 REAL :: bc, br, bi, bs, dust_extcoeff, hydro_extcoeff, extcoeff, vis_haze
2042 REAL :: tvd, rh, prob_ext_coeff_gt_p29, haze_ext_coeff
2043 REAL :: vis_hydlith, alpha_haze
2045 ! Dust factor based on 5 bin AFWA dust scheme. This is a simplification
2046 ! of the scheme in WRFPOST. More weight is applied to smaller particles.
2047 ! -----------------------------------------------------------------------
2048 dustfact=(/1.470E-6,7.877E-7,4.623E-7,2.429E-7,1.387E-7/)
2050 DO i=ims,ime
2051 DO j=jms,jme
2053 ! Hydrometeor extinction coefficient
2054 ! For now, lump graupel in with rain
2055 ! -------------------------------------------------------------------
2056 ! Update: GAC 20131213 Our test results with surface cloud and ice
2057 ! are very unfavorable. Model doesn't do a great job handling clouds
2058 ! at the model surface. Therefore, we will not trust surface
2059 ! cloud water/ice. (Commented out below).
2060 ! -------------------------------------------------------------------
2061 !br=2.240*qrain(i,j)**0.75
2062 !bs=10.36*qsnow(i,j)**0.78
2063 !bc=144.7*qcloud(i,j)**0.88
2064 !bi=327.8*qice(i,j)
2065 !hydro_extcoeff=bc+br+bi+bs
2066 !br=2.240*(qrain(i,j)+qgrpl(i,j))**0.75
2067 !bs=10.36*(qsnow(i,j)*rho(i,j))**0.78
2068 ! Update: moisture variables should be in mass concentration (g m^-3)
2069 br=1.1*(1000.*rho(i,j)*(qrain(i,j)+qgrpl(i,j)))**0.75
2070 bs=10.36*(1000.*rho(i,j)*qsnow(i,j))**0.78
2071 hydro_extcoeff=(br+bs)/1000. ! m^-1
2073 ! Dust extinction coefficient
2074 ! ---------------------------
2075 dust_extcoeff=0.
2076 DO d=1,ndust
2077 dust_extcoeff=dust_extcoeff+dustfact(d)*dustc(i,j,d)
2078 ENDDO
2080 ! UPDATE: GAC 20131213 Old algorithm commented out below
2081 !! Visibility due to haze obscuration
2082 !! -------------------------------------------------------
2083 !vis_haze=1500.*(105.-rh2m(i,j)+wind10m(i,j))
2084 !
2085 !! Calculate total visibility
2086 !! Take minimum visibility from hydro/lithometeors and haze
2087 !! Define maximum visibility as 20 km (UPDATE: 999.999 km)
2088 !! --------------------------------------------------------
2089 !extcoeff=hydro_extcoeff+dust_extcoeff
2090 !IF (extcoeff .gt. 0.) THEN
2091 ! vis(i,j)=MIN(visfactor/extcoeff,vis_haze)
2092 !ELSE
2093 ! vis(i,j)=999999.
2094 !ENDIF
2096 ! Update: GAC 20131213 New haze/fog visibility algorithm
2097 ! Start with relative humidity predictor. Increase
2098 ! predicted visibility as mixing ratio decreases (as
2099 ! there is less water to condense).
2100 ! -------------------------------------------------------
2101 vis_haze=999999.
2102 IF (q2m(i,j) .gt. 0.) THEN
2103 !vis_haze=1500.*(105.-rh2m(i,j))*(15./(1000.*q2m(i,j)))
2104 vis_haze=1500.*(105.-rh2m(i,j))*(5./min(1000.*q2m(i,j),5.))
2105 ENDIF
2107 ! Calculate a Weibull function "alpha" term. This can be
2108 ! used later with visibility (which acts as the "beta" term
2109 ! in the Weibull function) to create a probability distribution
2110 ! for visibility. Alpha can be thought of as the "level of
2111 ! certainty" that beta (model visibility) is correct. Fog is
2112 ! notoriously difficult to model. In the below algorithm,
2113 ! the alpha value (certainty) decreases as PWAT, mixing ratio,
2114 ! or winds decrease (possibly foggy conditions), but increases
2115 ! if RH decreases (more certainly not foggy). If PWAT is lower
2116 ! then there is a higher chance of radiation fog because there
2117 ! is less insulating cloud above.
2118 ! -------------------------------------------------------------
2119 alpha_haze=3.6
2120 IF (q2m(i,j) .gt. 0.) THEN
2121 alpha_haze=0.1 + pwater(i,j)/25. + wind125m(i,j)/3. + &
2122 (100.-rh2m(i,j))/10. + 1./(1000.*q2m(i,j))
2123 alpha_haze=min(alpha_haze,3.6)
2124 ENDIF
2126 ! Calculate visibility from hydro/lithometeors
2127 ! Maximum visibility -> 999999 meters
2128 ! --------------------------------------------
2129 extcoeff=hydro_extcoeff+dust_extcoeff
2130 IF (extcoeff .gt. 0.) THEN
2131 vis_hydlith=min(visfactor/extcoeff, 999999.)
2132 ELSE
2133 vis_hydlith=999999.
2134 ENDIF
2136 ! Calculate total visibility
2137 ! Take minimum visibility from hydro/lithometeors and haze
2138 ! Set alpha to be alpha_haze if haze dominates, or 3.6
2139 ! (a Guassian distribution) when hydro/lithometeors dominate
2140 ! ----------------------------------------------------------
2141 IF (vis_hydlith < vis_haze) THEN
2142 vis(i,j)=vis_hydlith
2143 vis_alpha(i,j)=3.6
2144 ELSE
2145 vis(i,j)=vis_haze
2146 vis_alpha(i,j)=alpha_haze
2147 ENDIF
2149 ! Calculate dust visibility
2150 ! Again, define maximum visibility as 20 km
2151 ! -----------------------------------------
2152 IF (dust_extcoeff .gt. 0.) THEN
2153 vis_dust(i,j)=MIN(visfactor/dust_extcoeff,999999.)
2154 ELSE
2155 vis_dust(i,j)=999999.
2156 ENDIF
2157 ENDDO
2158 ENDDO
2160 END SUBROUTINE vis_diagnostics
2164 SUBROUTINE cloud_diagnostics (qcloud &
2165 , qice &
2166 , qsnow &
2167 , rh &
2168 , dz8w &
2169 , rho &
2170 , z &
2171 , ht &
2172 , cloud &
2173 , cloud_ceil &
2174 , ids, ide, jds, jde, kds, kde &
2175 , ims, ime, jms, jme, kms, kme &
2176 , ips, ipe, jps, jpe, kps, kpe )
2178 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2179 ims, ime, jms, jme, kms, kme, &
2180 ips, ipe, jps, jpe, kps, kpe
2182 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
2183 INTENT(IN ) :: qcloud &
2184 , qice &
2185 , qsnow &
2186 , rh &
2187 , dz8w &
2188 , rho &
2189 , z
2191 REAL, DIMENSION( ims:ime, jms:jme ), &
2192 INTENT(IN ) :: ht
2194 REAL, DIMENSION( ims:ime, jms:jme ), &
2195 INTENT( OUT) :: cloud &
2196 , cloud_ceil
2198 ! Local
2199 ! -----
2200 INTEGER :: i, j, k
2201 REAL :: tot_cld_cond, maxrh, cld_frm_cnd, cld_frm_rh, z_maxrh
2202 REAL :: snow_extcoeff, vis_snow, cloud_lo, zagl_up, zagl_lo
2203 REAL, PARAMETER :: min_ceil = 125. ! Minimum ceiling of 125 m
2205 ! Calculate cloud cover based on total cloud condensate, or if none
2206 ! present, from maximum relative humidity in the column.
2207 ! -----------------------------------------------------------------
2208 DO i=ims,ime
2209 DO j=jms,jme
2211 ! Initialize some key variables
2212 ! -----------------------------
2213 tot_cld_cond = 0.
2214 cloud(i,j) = 0.
2215 maxrh = -9999.
2216 cloud_ceil(i,j) = -9999.
2217 cloud_lo = 0.
2219 ! Now go up the column to find our cloud base
2220 ! -------------------------------------------
2221 DO k=kms,kme
2223 ! Total cloud condensate
2224 ! Don't trust modeled cloud below 125 m. We
2225 ! let the alpha curve determine probabilities
2226 ! below 125 m...
2227 ! ---------------------------------------------
2228 IF ( z(i,k,j) - ht (i,j) .gt. min_ceil ) THEN
2230 ! Maximum column relative humidity
2231 ! --------------------------------
2232 IF (rh (i,k,j) .gt. maxrh) THEN
2233 maxrh = rh (i,k,j)
2234 z_maxrh = z(i,k,j)
2235 ENDIF
2237 ! ! Cloud cover parameterization. Take maximum value
2238 ! ! from condensate and relative humidity terms.
2239 ! ! ------------------------------------------------
2240 ! tot_cld_cond = tot_cld_cond + (qcloud (i,k,j) + qice (i,k,j) &
2241 ! + qsnow (i,k,j)) * dz8w (i,k,j) * rho(i,k,j)
2242 ! cld_frm_cnd = 50. * tot_cld_cond
2243 ! cld_frm_rh = MAX(((maxrh - 70.) / 30.),0.)
2244 ! cloud (i,j) = MAX(cld_frm_cnd,cld_frm_rh)
2246 ! Calculate cloud cover beta value by summing
2247 ! relative humidity above 70% as we go up the
2248 ! column. Assume a higher probability of a
2249 ! cloud if we have an accumulation of high
2250 ! relative humidity over a typical cloud
2251 ! depth of 500m. (Note dz8w is distance
2252 ! between full eta levels. Note also that rh
2253 ! is derived from the sum of qvapor, qcloud,
2254 ! and qcloud, with a maximum of 100%).
2255 ! -------------------------------------------
2256 cld_frm_rh = MAX(((rh (i,k,j) - 90.) / 10.),0.)
2257 cloud (i,j) = cloud (i,j) + ( cld_frm_rh * dz8w (i,k,j) ) / 250.
2259 ! Calculate cloud ceiling, the level at which
2260 ! parameterization of cloud cover exceeds 80%
2261 ! Once we exceed the 80% threshold, we will
2262 ! interpolate downward to find the ceiling.
2263 ! If this is the lowest level, then we will
2264 ! simply set ceiling to that level. After
2265 ! we interpolate, if ceiling is below the
2266 ! minimum we trust, we set it to the minimum.
2267 ! -------------------------------------------
2268 IF ( cloud_ceil (i,j) .eq. -9999. .and. cloud (i,j) .gt. 0.8 ) THEN
2269 zagl_up = z (i,k,j) - ht (i,j)
2270 IF ( k .EQ. kps ) THEN
2271 cloud_ceil (i,j) = zagl_up
2272 ELSE
2273 zagl_lo = z (i,k-1,j) - ht (i,j)
2274 cloud_ceil (i,j) = zagl_lo + &
2275 ((0.8 - cloud_lo) / &
2276 (cloud (i,j) - cloud_lo)) * &
2277 (zagl_up - zagl_lo)
2278 cloud_ceil (i,j) = MAX(cloud_ceil (i,j),ceil_min)
2279 ENDIF
2280 ENDIF
2282 ! Save cloud amount here to use for interpolation
2283 ! -----------------------------------------------
2284 cloud_lo=cloud(i,j)
2285 ENDIF
2286 ENDDO
2288 ! If we did not encounter any definitive cloud earlier, we set
2289 ! cloud ceiling to the level of maximum relative humidity. (If
2290 ! there is any cloud, this is our best guess as to where it
2291 ! will reside in the vertical). Height is AGL.
2292 ! ------------------------------------------------------------
2293 IF (cloud_ceil (i,j) .eq. -9999.) THEN
2294 cloud_ceil (i,j) = z_maxrh - ht (i,j)
2295 ENDIF
2297 ! Compare cloud ceiling to surface visibility reduction due to snow
2298 ! Note difference from horizontal visibility algorithm for snow
2299 ! -----------------------------------------------------------------
2300 IF (qsnow (i,1,j) .GT. 0. .AND. rho (i,1,j) .GT. 0.) THEN
2301 snow_extcoeff = 25. * (1000. * rho(i,1,j) * qsnow (i,1,j))**0.78
2302 snow_extcoeff = snow_extcoeff / 1000.
2303 vis_snow = 3.912 / snow_extcoeff
2304 IF (vis_snow .LT. cloud_ceil (i,j)) cloud_ceil (i,j) = vis_snow
2305 ENDIF
2307 ENDDO
2308 ENDDO
2310 END SUBROUTINE cloud_diagnostics
2314 SUBROUTINE thermal_diagnostics ( t2 &
2315 , psfc &
2316 , rh2m &
2317 , wind10m &
2318 , heatidx &
2319 , wchill &
2320 , fits &
2321 , ids, ide, jds, jde, kds, kde &
2322 , ims, ime, jms, jme, kms, kme &
2323 , ips, ipe, jps, jpe, kps, kpe )
2325 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2326 ims, ime, jms, jme, kms, kme, &
2327 ips, ipe, jps, jpe, kps, kpe
2329 REAL, DIMENSION( ims:ime, jms:jme ), &
2330 INTENT(IN ) :: t2 &
2331 , psfc &
2332 , rh2m &
2333 , wind10m
2335 REAL, DIMENSION( ims:ime, jms:jme ), &
2336 INTENT( OUT) :: heatidx &
2337 , wchill &
2338 , fits
2340 ! Local
2341 ! -----
2342 INTEGER :: i, j
2344 DO i=ims,ime
2345 DO j=jms,jme
2347 ! Heat Index
2348 ! ----------
2349 heatidx ( i, j ) = calc_hi ( t2 ( i, j ) &
2350 , rh2m ( i, j ) )
2352 ! Wind Chill
2353 ! ----------
2354 wchill ( i, j ) = calc_wc ( t2 ( i, j ) &
2355 , wind10m ( i, j ) )
2357 ! Fighter Index of Thermal Stress
2358 ! -------------------------------
2359 fits ( i, j ) = calc_fits ( psfc ( i, j ) &
2360 , t2 ( i, j ) &
2361 , rh2m ( i, j ) )
2363 ENDDO
2364 ENDDO
2366 END SUBROUTINE thermal_diagnostics
2370 SUBROUTINE turbulence_diagnostics ( u_phy &
2371 , v_phy &
2372 , t_phy &
2373 , p &
2374 , zagl &
2375 , defor11 &
2376 , defor12 &
2377 , defor22 &
2378 , turb &
2379 , llturb &
2380 , llturblgt &
2381 , llturbmdt &
2382 , llturbsvr &
2383 , nlyrs &
2384 , lyrbot &
2385 , lyrtop &
2386 , ids, ide, jds, jde, kds, kde &
2387 , ims, ime, jms, jme, kms, kme &
2388 , ips, ipe, jps, jpe, kps, kpe )
2390 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
2391 ims, ime, jms, jme, kms, kme, &
2392 ips, ipe, jps, jpe, kps, kpe
2394 INTEGER, INTENT(IN) :: nlyrs
2396 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
2397 INTENT(IN ) :: u_phy &
2398 , v_phy &
2399 , t_phy &
2400 , p &
2401 , zagl &
2402 , defor11 &
2403 , defor12 &
2404 , defor22
2406 REAL, DIMENSION( nlyrs ), &
2407 INTENT(IN ) :: lyrtop &
2408 , lyrbot
2410 REAL, DIMENSION( ims:ime, nlyrs, jms:jme ), &
2411 INTENT( OUT) :: turb
2413 REAL, DIMENSION( ims:ime, jms:jme ), &
2414 INTENT( OUT) :: llturb &
2415 , llturblgt &
2416 , llturbmdt &
2417 , llturbsvr
2419 ! Local
2420 ! -----
2421 INTEGER :: i, j, k, n, bot, top, nlayer
2423 REAL :: ugrdtop &
2424 , ugrdbot &
2425 , vgrdtop &
2426 , vgrdbot &
2427 , defor11top &
2428 , defor11bot &
2429 , defor12top &
2430 , defor12bot &
2431 , defor22top &
2432 , defor22bot
2434 REAL, DIMENSION( kms:kme ) :: this_zagl &
2435 , this_tK &
2436 , this_p &
2437 , this_u &
2438 , this_v
2440 REAL :: wind, therm, mtn_wave, tpd_wave
2442 !~ Initialize variables.
2443 ! ----------------------
2444 turb = REAL ( 0 )
2445 llturb = REAL ( 0 )
2446 llturblgt = REAL ( 0 )
2447 llturbsvr = REAL ( 0 )
2449 !~ Loop through the grid.
2450 ! ----------------------
2451 DO i=ims,ime
2452 DO j=jms,jme
2454 !~ Loop through the turbulence layers
2455 ! ----------------------------------
2456 DO n = 1, nlyrs
2458 !~ Interpolate relevent variables to turbulence layers
2459 ! ---------------------------------------------------
2460 ugrdtop = lin_interp ( zagl ( i, kms:kme-1, j ) &
2461 , u_phy ( i, kms:kme-1, j ) &
2462 , lyrtop ( n ) )
2463 ugrdbot = lin_interp ( zagl ( i, kms:kme-1, j ) &
2464 , u_phy ( i, kms:kme-1, j ) &
2465 , lyrbot ( n ) )
2466 vgrdtop = lin_interp ( zagl ( i, kms:kme-1, j ) &
2467 , v_phy ( i, kms:kme-1, j ) &
2468 , lyrtop ( n ) )
2469 vgrdbot = lin_interp ( zagl ( i, kms:kme-1, j ) &
2470 , v_phy ( i, kms:kme-1, j ) &
2471 , lyrbot ( n ) )
2472 defor11top = lin_interp ( zagl ( i, kms:kme-1, j ) &
2473 , defor11 ( i, kms:kme-1, j ) &
2474 , lyrtop ( n ) )
2475 defor11bot = lin_interp ( zagl ( i, kms:kme-1, j ) &
2476 , defor11 ( i, kms:kme-1, j ) &
2477 , lyrbot ( n ) )
2478 defor12top = lin_interp ( zagl ( i, kms:kme-1, j ) &
2479 , defor12 ( i, kms:kme-1, j ) &
2480 , lyrtop ( n ) )
2481 defor12bot = lin_interp ( zagl ( i, kms:kme-1, j ) &
2482 , defor12 ( i, kms:kme-1, j ) &
2483 , lyrbot ( n ) )
2484 defor22top = lin_interp ( zagl ( i, kms:kme-1, j ) &
2485 , defor22 ( i, kms:kme-1, j ) &
2486 , lyrtop ( n ) )
2487 defor22bot = lin_interp ( zagl ( i, kms:kme-1, j ) &
2488 , defor22 ( i, kms:kme-1, j ) &
2489 , lyrbot ( n ) )
2491 !~ Compute Knapp-Ellrod clear air turbulence
2492 ! -----------------------------------------
2493 turb ( i, n, j ) = CATTurbulence ( ugrdbot &
2494 , ugrdtop &
2495 , vgrdbot &
2496 , vgrdtop &
2497 , defor11bot &
2498 , defor11top &
2499 , defor12bot &
2500 , defor12top &
2501 , defor22bot &
2502 , defor22top &
2503 , lyrbot (n) &
2504 , lyrtop (n) )
2506 ENDDO
2508 !~ Get top and bottom index of 0-1500 m AGL layer
2509 ! ----------------------------------------------
2510 bot = kms
2511 top = kms
2512 DO k=kms+1,kme
2513 IF ( zagl ( i, k, j ) .gt. 1500. ) THEN
2514 top = k
2515 EXIT
2516 ENDIF
2517 ENDDO
2518 nlayer = top - bot + 1 !~ Number of layers from bottom to top
2520 !~ Copy current column at this i,j point into working arrays
2521 ! ---------------------------------------------------------
2522 this_zagl = zagl ( i, kms:kme, j )
2523 this_tK = t_phy ( i, kms:kme, j )
2524 this_p = p ( i, kms:kme, j )
2525 this_u = u_phy ( i, kms:kme, j )
2526 this_v = v_phy ( i, kms:kme, j )
2528 !~ Interpolate req'd vars to the 1500 m level
2529 !~ Overwrite the "top" index with these values
2530 ! -------------------------------------------
2531 this_zagl ( top ) = 1500.
2532 this_tK ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
2533 , t_phy ( i, kms:kme-1, j ) &
2534 , this_zagl (top) )
2535 this_p ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
2536 , p ( i, kms:kme-1, j ) &
2537 , this_zagl (top) )
2538 this_u ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
2539 , u_phy ( i, kms:kme-1, j ) &
2540 , this_zagl (top) )
2541 this_v ( top ) = lin_interp ( zagl ( i, kms:kme-1, j ) &
2542 , v_phy ( i, kms:kme-1, j ) &
2543 , this_zagl (top) )
2545 !~ Compute the low level turbulence index (from 0 - 1500 m AGL)
2546 !~ There are four components to this index: a wind speed term,
2547 !~ thermodynamic term, mountain wave term, and trapped wave term.
2548 !~ These terms will utilize the working arrays we have just
2549 !~ defined, using only the portion of the array valid in the
2550 !~ 0-1500 m layer, e.g. this_u (bot:top).
2551 !~ The algorithm itself was developed by:
2552 !~
2553 !~ Mr. James McCormick
2554 !~ Aviation Hazards Team
2555 !~ Air Force Weather Agency 16WS/WXN
2556 !~ DSN: 271-1689 Comm: (402) 294-1689
2557 !~ James.McCormick.ctr@offutt.af.mil
2558 ! -------------------------------------------------------------- !
2559 !~ Step 1: Compute the wind speed term ~!
2560 ! -------------------------------------------------------------- !
2561 wind = LLT_WindSpeed ( nlayer, this_u (bot:top) &
2562 , this_v (bot:top) )
2564 ! -------------------------------------------------------------- !
2565 !~ Step 2: Compute the thermodynamic term ~!
2566 ! -------------------------------------------------------------- !
2567 therm = LLT_Thermodynamic ( nlayer, this_tK(bot:top) &
2568 , this_zagl(bot:top) )
2570 ! -------------------------------------------------------------- !
2571 !~ Step 3: Compute the mountain wave term ~!
2572 ! -------------------------------------------------------------- !
2573 mtn_wave = LLT_MountainWave ( nlayer, terrain_dx, terrain_dy &
2574 , this_u(bot:top), this_v(bot:top) &
2575 , this_tK(bot:top), this_zagl(bot:top) )
2577 ! -------------------------------------------------------------- !
2578 !~ Step 4: Compute the trapped wave term ~!
2579 ! -------------------------------------------------------------- !
2580 tpd_wave = LLT_TrappedWave ( nlayer, this_u(bot:top) &
2581 , this_v(bot:top), this_p(bot:top) )
2583 ! -------------------------------------------------------------- !
2584 !~ Step 5: Combine the above and arrive at the turbulence index. ~!
2585 ! -------------------------------------------------------------- !
2586 llturb ( i,j ) = 1.-((1.-wind)*(1.-therm)*(1.-mtn_wave)*(1.-tpd_wave))
2588 !~ Compute probabilities of light, moderate, and severe LLT
2589 ! --------------------------------------------------------
2590 llturblgt ( i,j ) = (((((((llturb (i,j) * REAL (100))-REAL (30)) &
2591 *2.5)*.01)**2)*0.75)*REAL(100))
2592 IF ( llturb (i,j) < 0.3 ) llturblgt ( i,j ) = REAL ( 0 )
2593 IF ( llturblgt (i,j) > REAL (90) ) llturblgt ( i,j ) = REAL ( 90 )
2595 llturbmdt ( i,j ) = (((((((llturb (i,j) * REAL (100))-REAL (35)) &
2596 *2.22222)*.01)*0.75)**2)*88.88888)
2597 IF ( llturb (i,j) < 0.35 ) llturbmdt ( i,j ) = REAL ( 0 )
2598 IF ( llturbmdt (i,j) > REAL (70) ) llturbmdt ( i,j ) = REAL ( 70 )
2600 llturbsvr ( i,j ) = (((((((llturb (i,j) * REAL (100))-REAL (40)) &
2601 *REAL(2))*.01)*0.5)**2)*REAL(100))
2602 IF ( llturb (i,j) < 0.40 ) llturbsvr ( i,j ) = REAL ( 0 )
2603 IF ( llturbsvr (i,j) > REAL (35) ) llturbsvr ( i,j ) = REAL ( 35 )
2605 ENDDO
2606 ENDDO
2608 END SUBROUTINE turbulence_diagnostics
2610 END MODULE module_diag_afwa
2612 #endif