| blob | 098e9d9655e29051b7be9d451268ef62e03488df |
1 #if (NMM_CORE == 1)
2 MODULE module_diag_hailcast
3 CONTAINS
4 SUBROUTINE diag_hailcast_stub
5 END SUBROUTINE diag_hailcast_stub
6 END MODULE module_diag_hailcast
7 #else
9 MODULE module_diag_hailcast
11 CONTAINS
13 SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags &
14 , moist &
15 , rho &
16 , ids, ide, jds, jde, kds, kde &
17 , ims, ime, jms, jme, kms, kme &
18 , ips, ipe, jps, jpe, kps, kpe &
19 , its, ite, jts, jte &
20 , k_start, k_end &
21 , dt, itimestep &
22 , haildt,curr_secs &
23 , haildtacttime )
25 USE module_domain, ONLY : domain , domain_clock_get
26 USE module_configure, ONLY : grid_config_rec_type, model_config_rec
27 USE module_state_description
28 USE module_model_constants
29 USE module_utility
30 USE module_streams, ONLY: history_alarm, auxhist2_alarm
31 #ifdef DM_PARALLEL
32 USE module_dm, ONLY: wrf_dm_sum_real, wrf_dm_maxval
33 #endif
35 IMPLICIT NONE
37 TYPE ( domain ), INTENT(INOUT) :: grid
38 TYPE ( grid_config_rec_type ), INTENT(IN) :: config_flags
40 INTEGER, INTENT(IN) :: ids, ide, jds, jde, kds, kde, &
41 ims, ime, jms, jme, kms, kme, &
42 ips, ipe, jps, jpe, kps, kpe
43 INTEGER :: k_start , k_end, its, ite, jts, jte
45 REAL, DIMENSION( ims:ime, kms:kme, jms:jme , num_moist), &
46 INTENT(IN ) :: moist
48 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
49 INTENT(IN ) :: rho
51 REAL, INTENT(IN ),OPTIONAL :: haildt
52 REAL, INTENT(IN ),OPTIONAL :: curr_secs
53 REAL, DIMENSION( ims:ime, jms:jme ), &
54 INTENT(INOUT),OPTIONAL :: haildtacttime
55 INTEGER, INTENT(IN ) :: itimestep
56 REAL, INTENT(IN ) :: dt
59 ! Local
60 ! -----
61 CHARACTER*512 :: message
62 CHARACTER*256 :: timestr
63 INTEGER :: i,j,k,nz
64 INTEGER :: i_start, i_end, j_start, j_end
65 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) :: qr &
66 , qs &
67 , qg &
68 , qv &
69 , qc &
70 , qi &
71 , ptot
72 REAL, DIMENSION( ims:ime, jms:jme ) :: wup_mask_prev &
73 , wdur_prev
74 REAL :: dhail1,dhail2,dhail3,dhail4,dhail5
76 ! Timing
77 ! ------
78 TYPE(WRFU_Time) :: hist_time, aux2_time, CurrTime, StartTime
79 TYPE(WRFU_TimeInterval) :: dtint, histint, aux2int
80 LOGICAL :: is_after_history_dump, is_output_timestep, is_first_timestep
81 LOGICAL :: doing_adapt_dt, run_param, decided
82 INTEGER :: stephail, itimestep_basezero
85 ! Chirp the routine name for debugging purposes
86 ! ---------------------------------------------
87 write ( message, * ) 'inside hailcast_diagnostics_driver'
88 CALL wrf_debug( 100 , message )
90 ! Get timing info
91 ! Want to know if when the last history output was
92 ! Check history and auxhist2 alarms to check last ring time and how often
93 ! they are set to ring
94 ! -----------------------------------------------------------------------
95 CALL WRFU_ALARMGET( grid%alarms( HISTORY_ALARM ), prevringtime=hist_time, &
96 ringinterval=histint)
97 CALL WRFU_ALARMGET( grid%alarms( AUXHIST2_ALARM ), prevringtime=aux2_time, &
98 ringinterval=aux2int)
100 ! Get domain clock
101 ! ----------------
102 CALL domain_clock_get ( grid, current_time=CurrTime, &
103 simulationStartTime=StartTime, &
104 current_timestr=timestr, time_step=dtint )
106 ! Set some booleans for use later
107 ! Following uses an overloaded .lt.
108 ! ---------------------------------
109 is_after_history_dump = ( Currtime .lt. hist_time + dtint )
111 ! Following uses an overloaded .ge.
112 ! ---------------------------------
113 is_output_timestep = (Currtime .ge. hist_time + histint - dtint .or. &
114 Currtime .ge. aux2_time + aux2int - dtint )
115 write ( message, * ) 'is output timestep? ', is_output_timestep
116 CALL wrf_debug( 100 , message )
118 ! Following uses an overloaded .eq.
119 ! ---------------------------------
120 is_first_timestep = ( Currtime .eq. StartTime + dtint )
123 ! After each history dump, reset max/min value arrays
124 ! ----------------------------------------------------------------------
125 IF ( is_after_history_dump ) THEN
126 DO j = jts, jte
127 DO i = its, ite
128 grid%hailcast_dhail1(i,j) = 0.
129 grid%hailcast_dhail2(i,j) = 0.
130 grid%hailcast_dhail3(i,j) = 0.
131 grid%hailcast_dhail4(i,j) = 0.
132 grid%hailcast_dhail5(i,j) = 0.
133 ENDDO
134 ENDDO
135 ENDIF ! is_after_history_dump
138 ! We need to do some neighboring gridpoint comparisons for the updraft
139 ! duration calculations; set i,j start and end values so we don't go off
140 ! the edges of the domain. Updraft duration on domain edges will always be 0.
141 ! ----------------------------------------------------------------------
142 i_start = its
143 i_end = ite
144 j_start = jts
145 j_end = jte
147 IF ( config_flags%open_xs .OR. config_flags%specified .OR. &
148 config_flags%nested) i_start = MAX( ids+1, its )
149 IF ( config_flags%open_xe .OR. config_flags%specified .OR. &
150 config_flags%nested) i_end = MIN( ide-2, ite ) !-1
151 IF ( config_flags%open_ys .OR. config_flags%specified .OR. &
152 config_flags%nested) j_start = MAX( jds+1, jts )
153 IF ( config_flags%open_ye .OR. config_flags%specified .OR. &
154 config_flags%nested) j_end = MIN( jde-2, jte ) !-1
155 IF ( config_flags%periodic_x ) i_start = its
156 IF ( config_flags%periodic_x ) i_end = ite
159 ! Make a copy of the updraft duration, mask variables
160 ! ---------------------------------------------------
161 !wdur_prev(:,:) = grid%hailcast_wdur(:,:)
162 !wup_mask_prev(:,:) = grid%hailcast_wup_mask(:,:)
163 DO j = MAX( jts-1 , jds ), MIN( jte+1 , jde-1 )
164 DO i = MAX( its-1 , ids ), MIN( ite+1 , ide-1)
165 wdur_prev(i,j) = grid%hailcast_wdur(i,j)
166 wup_mask_prev(i,j) = grid%hailcast_wup_mask(i,j)
167 ENDDO
168 ENDDO
171 ! Determine updraft mask (where updraft greater than some threshold)
172 ! ---------------------------------------------------
173 DO j = jts, jte
174 DO i = its, ite
175 grid%hailcast_wup_mask(i,j) = 0
176 grid%hailcast_wdur(i,j) = 0
178 DO k = k_start, k_end
179 IF ( grid%w_2(i,k,j) .ge. 10. ) THEN
180 grid%hailcast_wup_mask(i,j) = 1
181 ENDIF
182 ENDDO
183 ENDDO
184 ENDDO
186 ! Determine updraft duration; make sure not to call point outside the domain
187 ! ---------------------------------------------------
188 DO j = j_start, j_end
189 DO i = i_start, i_end
191 ! Determine updraft duration using updraft masks
192 ! ---------------------------------------------------
193 IF ( (grid%hailcast_wup_mask(i,j).eq.1) .OR. &
194 (MAXVAL(wup_mask_prev(i-1:i+1,j-1:j+1)).eq.1) ) THEN
195 grid%hailcast_wdur(i,j) = &
196 MAXVAL(wdur_prev(i-1:i+1,j-1:j+1)) + grid%dt
197 ENDIF
198 ENDDO
199 ENDDO
202 ! Only run the scheme every "haildt" seconds as set in the namelist.
203 ! Code borrowed from module_pbl_driver, although here haildt is
204 ! in seconds, not minutes (bldt is in minutes).
205 ! ---------------------------------------------------
207 ! Are we using adaptive timestepping?
208 doing_adapt_dt = .FALSE.
209 IF ( (config_flags%use_adaptive_time_step) .and. &
210 ( (.not. grid%nested) .or. &
211 ( (grid%nested) .and. (abs(grid%dtbc) < 0.0001) ) ) )THEN
212 doing_adapt_dt = .TRUE.
214 !170519 - bug fix RAS
215 !Change haildtacttime to an array so can be set individually
216 ! by gridpoint to account for quilted processes
217 DO j = j_start, j_end
218 DO i = i_start, i_end
219 IF ( haildtacttime(i,j) .eq. 0. ) THEN
220 haildtacttime(i,j) = CURR_SECS + haildt
221 END IF
222 ENDDO
223 ENDDO
224 END IF
226 ! Calculate STEPHAIL
227 stephail = NINT(haildt / dt)
228 stephail = MAX(stephail,1)
230 ! itimestep starts at 1 - we want it to start at 0 so correctly
231 ! divisibile by stephail.
232 itimestep_basezero = itimestep -1
235 ! Do we run through this scheme or not?
236 ! Test 1: If this is the initial model time, then yes.
237 ! ITIMESTEP=1
238 ! Test 2: If the user asked for hailcast to be run every time step, then yes.
239 ! HAILDT=0 or STEPHAIL=1
240 ! Test 3: If not adaptive dt, and this is on the requested hail frequency,
241 ! then yes.
242 ! MOD(ITIMESTEP,STEPHAIL)=0
243 ! Test 4: If using adaptive dt and the current time is past the last
244 ! requested activate hailcast time, then yes.
245 ! CURR_SECS >= HAILDTACTTIME
247 ! If we do run through the scheme, we set the flag run_param to TRUE and we set
248 ! the decided flag to TRUE. The decided flag says that one of these tests was
249 ! able to say "yes", run the scheme. We only proceed to other tests if the
250 ! previous tests all have left decided as FALSE. If we set run_param to TRUE
251 ! and this is adaptive time stepping, we set the time to the next hailcast run.
253 run_param = .FALSE.
254 decided = .FALSE.
255 IF ( ( .NOT. decided ) .AND. &
256 ( itimestep_basezero .EQ. 0 ) ) THEN
257 run_param = .TRUE.
258 decided = .TRUE.
259 END IF
261 IF ( PRESENT(haildt) )THEN
262 IF ( ( .NOT. decided ) .AND. &
263 ( ( haildt .EQ. 0. ) .OR. ( stephail .EQ. 1 ) ) ) THEN
264 run_param = .TRUE.
265 decided = .TRUE.
266 END IF
267 ELSE
268 IF ( ( .NOT. decided ) .AND. &
269 ( stephail .EQ. 1 ) ) THEN
270 run_param = .TRUE.
271 decided = .TRUE.
272 END IF
273 END IF
275 IF ( ( .NOT. decided ) .AND. &
276 ( .NOT. doing_adapt_dt ) .AND. &
277 ( MOD(itimestep_basezero,stephail) .EQ. 0 ) ) THEN
278 run_param = .TRUE.
279 decided = .TRUE.
280 END IF
282 !170519 - bug fix RAS
283 !Changed haildtacttime to an array so can be set individually
284 ! by gridpoint to account for quilted processes
285 IF ( ( .NOT. decided ) .AND. &
286 ( doing_adapt_dt ) .AND. &
287 ( curr_secs .GE. haildtacttime(i_start, j_start) ) ) THEN
288 run_param = .TRUE.
289 decided = .TRUE.
290 DO j = j_start, j_end
291 DO i = i_start, i_end
292 haildtacttime(i,j) = curr_secs + haildt
293 ENDDO
294 ENDDO
295 END IF
297 write ( message, * ) 'timestep to run HAILCAST? ', run_param
298 CALL wrf_debug( 100 , message )
300 IF (run_param) THEN
302 ! 3-D arrays for moisture variables
303 ! ---------------------------------
304 DO i=its, ite !ims, ime
305 DO k=kms, kme
306 DO j=jts, jte !jms, jme
307 qv(i,k,j) = moist(i,k,j,P_QV)
308 qr(i,k,j) = moist(i,k,j,P_QR)
309 qs(i,k,j) = moist(i,k,j,P_QS)
310 qg(i,k,j) = moist(i,k,j,P_QG)
311 qc(i,k,j) = moist(i,k,j,P_QC)
312 qi(i,k,j) = moist(i,k,j,P_QI)
313 ENDDO
314 ENDDO
315 ENDDO
317 ! Total pressure
318 ! --------------
319 DO i=its, ite! ims, ime
320 DO k=kms, kme
321 DO j=jts, jte !jms, jme
322 ptot(i,k,j)=grid%pb(i,k,j)+grid%p(i,k,j)
323 ENDDO
324 ENDDO
325 ENDDO
327 ! Hail diameter in millimeters (HAILCAST)
328 ! ---------------------------------------------------
329 nz = k_end - k_start
330 DO j = jts, jte
331 DO i = its, ite
333 ! Only call hailstone driver if updraft has been
334 ! around longer than 15 min
335 ! ----------------------------------------------
336 IF (grid%hailcast_wdur(i,j) .gt. 900) THEN
337 CALL hailstone_driver ( grid%t_phy(i,kms:kme,j), &
338 grid%z(i,kms:kme,j), &
339 grid%ht(i, j), &
340 ptot(i,kms:kme,j), &
341 rho(i,kms:kme,j), &
342 qv(i,kms:kme,j), &
343 qi(i,kms:kme,j), &
344 qc(i,kms:kme,j), &
345 qr(i,kms:kme,j), &
346 qs(i,kms:kme,j), &
347 qg(i,kms:kme,j), &
348 grid%w_2(i,kms:kme,j), &
349 grid%hailcast_wdur(i,j), &
350 nz, &
351 dhail1, dhail2, &
352 dhail3, dhail4, &
353 dhail5 )
354 IF (dhail1 .gt. grid%hailcast_dhail1(i,j)) THEN
355 grid%hailcast_dhail1(i,j) = dhail1
356 ENDIF
357 IF (dhail2 .gt. grid%hailcast_dhail2(i,j)) THEN
358 grid%hailcast_dhail2(i,j) = dhail2
359 ENDIF
360 IF (dhail3 .gt. grid%hailcast_dhail3(i,j)) THEN
361 grid%hailcast_dhail3(i,j) = dhail3
362 ENDIF
363 IF (dhail4 .gt. grid%hailcast_dhail4(i,j)) THEN
364 grid%hailcast_dhail4(i,j) = dhail4
365 ENDIF
366 IF (dhail5 .gt. grid%hailcast_dhail5(i,j)) THEN
367 grid%hailcast_dhail5(i,j) = dhail5
368 ENDIF
369 ENDIF
370 ENDDO
371 ENDDO
373 ! Calculate the mean and standard deviation of the hail diameter
374 ! distribution over different embryo sizes
375 ! ----------------------------------------
376 DO j = jts, jte !jms, jme
377 DO i = its, ite !ims, ime
378 !mean
379 grid%hailcast_diam_mean(i,j)=(grid%hailcast_dhail1(i,j)+&
380 grid%hailcast_dhail2(i,j) +grid%hailcast_dhail3(i,j)+&
381 grid%hailcast_dhail4(i,j) +grid%hailcast_dhail5(i,j))/5.
382 !max
383 grid%hailcast_diam_max(i,j)=MAX(grid%hailcast_dhail1(i,j),&
384 grid%hailcast_dhail2(i,j), grid%hailcast_dhail3(i,j),&
385 grid%hailcast_dhail4(i,j), grid%hailcast_dhail5(i,j))
386 !sample standard deviation
387 grid%hailcast_diam_std(i,j) = SQRT( ( &
388 (grid%hailcast_dhail1(i,j)-grid%hailcast_diam_mean(i,j))**2.+&
389 (grid%hailcast_dhail2(i,j)-grid%hailcast_diam_mean(i,j))**2.+&
390 (grid%hailcast_dhail3(i,j)-grid%hailcast_diam_mean(i,j))**2.+&
391 (grid%hailcast_dhail4(i,j)-grid%hailcast_diam_mean(i,j))**2.+&
392 (grid%hailcast_dhail5(i,j)-grid%hailcast_diam_mean(i,j))**2.)&
393 / 4.0)
394 ENDDO
395 ENDDO
397 END IF
399 END SUBROUTINE hailcast_diagnostic_driver
401 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
402 !!!!
403 !!!! Hailstone driver, adapted from hailstone subroutine in HAILCAST
404 !!!! Inputs:
405 !!!! 1-d (nz)
406 !!!! TCA temperature (K)
407 !!!! h1d height above sea level (m)
408 !!!! PA total pressure (Pa)
409 !!!! rho1d density (kg/m3)
410 !!!! RA vapor mixing ratio (kg/kg)
411 !!!! qi1d cloud ice mixing ratio (kg/kg)
412 !!!! qc1d cloud water mixing ratio (kg/kg)
413 !!!! qr1d rain water mixing ratio (kg/kg)
414 !!!! qg1d graupel mixing ratio (kg/kg)
415 !!!! qs1d snow mixing ratio (kg/kg)
416 !!!! VUU updraft speed at each level (m/s)
417 !!!! Float
418 !!!! ht terrain height (m)
419 !!!! wdur duration of any updraft > 10 m/s within 1 surrounding
420 !!!! gridpoint
421 !!!! Integer
422 !!!! nz number of vertical levels
423 !!!!
424 !!!! Output:
425 !!!! dhail hail diameter in mm
426 !!!! 1st-5th rank-ordered hail diameters returned
427 !!!!
428 !!!! 13 Aug 2013 .................................Becky Adams-Selin AER
429 !!!! adapted from hailstone subroutine in SPC's HAILCAST
430 !!!! 18 Mar 2014 .................................Becky Adams-Selin AER
431 !!!! added variable rime layer density, per Ziegler et al. (1983)
432 !!!! 4 Jun 2014 ..................................Becky Adams-Selin AER
433 !!!! removed initial embryo size dependency on microphysics scheme
434 !!!! 5 Jun 2014 ..................................Becky Adams-Selin AER
435 !!!! used smaller initial embryo sizes
436 !!!! 25 Jun 2015..................................Becky Adams-Selin AER
437 !!!! Significant revamping. Fixed units bug in HEATBUD that caused
438 !!!! hailstone temperature instabilities. Similar issue fixed in BREAKUP
439 !!!! subroutine. Removed graupel from ice content. Changed initial
440 !!!! embryo size and location to better match literature. Added
441 !!!! enhanced melting when hailstone collides with liquid water
442 !!!! in regions above freezing. Final diameter returned is ice diameter
443 !!!! only. Added hailstone temperature averaging over previous timesteps
444 !!!! to decrease initial temperature instability at small embyro diameters.
445 !!!! 3 Sep 2015...................................Becky Adams-Selin AER
446 !!!! Insert embryos at -13C; interpret pressure and other variables to
447 !!!! that exact temperature level.
448 !!!! 16 Nov 2015...................................Becky Adams-Selin AER
449 !!!! Hailstone travels horizontally through updraft instead of being
450 !!!! locked in the center.
451 !!!! 9 Jul 2016...................................Becky Adams-Selin AER
452 !!!! Uses an adiabatic liquid cloud water profile generated from
453 !!!! the saturation vapor pressure using the model temperature
454 !!!! profile.
455 !!!! 04 Feb 2017...................................Becky Adams-Selin AER
456 !!!! Added adaptive time-stepping option.
457 !!!! ***Don't set any higher than 60 seconds***
458 !!!! 06 May 2017...................................Becky Adams-Selin AER
459 !!!! Bug fixes for some memory issues in the adiabatic liquid
460 !!!! water profile code.
461 !!!! 23 Apr 2019...................................Becky Adams-Selin AER
462 !!!! Added check to make sure embryo spends at least some time in
463 !!!! cloud before returning a positive hail diameter.
464 !!!! 16 Dec 2022...................................Becky Adams-Selin AER
465 !!!! Consolodate a number of bug fixes, including fixing behavior
466 !!!! of haildt when tiling is used, and non-SI units errors in
467 !!!! HEATBUD subroutine.
468 !!!!
469 !!!! See Adams-Selin and Ziegler 2016, MWR for further documentation.
470 !!!!
471 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
472 SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,&
473 RA, qi1d,qc1d,qr1d,qs1d,qg1d, &
474 VUU, wdur, &
475 nz,dhail1,dhail2,dhail3,dhail4, &
476 dhail5 )
477 IMPLICIT NONE
478 INTEGER, INTENT(IN ) :: nz
480 REAL, DIMENSION( nz ), &
481 INTENT(IN ) :: TCA & ! temperature (K)
482 , rho1d &
483 , h1d &
484 , PA & ! pressure (Pa)
485 , RA & ! vapor mixing ratio (kg/kg)
486 , VUU & ! updraft speed (m/s)
487 , qi1d,qc1d,qr1d &
488 , qs1d,qg1d
490 REAL, INTENT(IN ) :: ht &
491 , wdur
493 !Output: 1st-5th rank-ordered hail diameters returned
494 REAL, INTENT(INOUT) :: dhail1 & ! hail diameter (mm);
495 , dhail2 &
496 , dhail3 &
497 , dhail4 &
498 , dhail5
499 !Local variables
500 REAL ZBAS, TBAS, WBASP ! height, temp, pressure of cloud base
501 REAL RBAS ! mix ratio of cloud base
502 REAL cwitot ! total cloud water, ice mix ratio
503 INTEGER KBAS ! k of cloud base
504 REAL tk_embryo ! temperature at which initial embryo is inserted
505 REAL ZFZL, TFZL, WFZLP ! height, temp, pressure of embryo start point
506 REAL RFZL ! mix ratio of embryo start point
507 REAL VUFZL, DENSAFZL ! updraft speed, density of embryo start point
508 INTEGER KFZL ! k of embryo start point
509 INTEGER nofroze ! keeps track if hailstone has ever been frozen
510 INTEGER CLOUDON ! use to zero out cloud water, ice once past updraft duration
511 REAL RTIME ! real updraft duration (sec)
512 REAL TAU, TAU_1, TAU_2 ! upper time limit of simulation (sec)
513 REAL delTAU ! difference between TAU_2 and TAU_1 (sec)
514 REAL g ! gravity (m/s)
515 REAL r_d ! constant
516 !hailstone parameters
517 REAL*8 DD, D, D_ICE ! hail diameter (m)
518 REAL VT ! terminal velocity (m/s)
519 REAL V ! actual stone velocity (m/s)
520 REAL TS ! hailstone temperature (K)
521 !HAILSTONE temperature differencing
522 REAL TSm1, TSm2 ! hailstone temperature at previous 3 timesteps
523 REAL FW ! fraction of stone that is liquid
524 REAL WATER ! mass of stone that is liquid
525 REAL CRIT ! critical water mass allowed on stone surface
526 REAL DENSE ! hailstone density (kg/m3)
527 INTEGER ITYPE ! wet (2) or dry (1) growth regime
528 !1-d column arrays of updraft parameters
529 REAL, DIMENSION( nz ) :: &
530 RIA, & ! frozen content mix ratio (kg/kg)
531 RWA, & ! liquid content mix ratio (kg/kg)
532 VUU_pert ! perturbed updraft profile (m/s)
533 !1-d column arrays of updraft characteristics for adiabatic water profile
534 REAL, DIMENSION( nz ) :: &
535 RWA_adiabat, & ! adiabatic liquid content mixing ratio (kg/kg)
536 RWA_new, &
537 ESVA, & ! saturation vapor pressure
538 RSA ! saturation mixing ratio
539 !in-cloud updraft parameters at location of hailstone
540 REAL P ! in-cloud pressure (Pa)
541 REAL RS ! in-cloud saturation mixing ratio
542 REAL RI, RW ! ice, liquid water mix. ratio (kg/kg)
543 REAL XI, XW ! ice, liquid water content (kg/m3 air)
544 REAL PC ! in-cloud fraction of frozen water
545 REAL TC ! in-cloud temperature (K)
546 REAL VU ! in-cloud updraft speed (m/s)
547 REAL VUMAX ! in-cloud updraft speed read from WRF (max allowed)
548 REAL VUCORE ! perturbed in-cloud updraft speed
549 REAL DENSA ! in-cloud updraft density (kg/m3)
550 REAL Z ! height of hailstone (m)
551 REAL DELRW ! diff in sat vap. dens. between hail and air (kg/m3)
552 REAL, DIMENSION(5) :: dhails !hail diameters from 5 different embryo size
553 !mean sub-cloud layer variables
554 REAL TLAYER,RLAYER,PLAYER ! mean sub-cloud temp, mix ratio, pres
555 REAL TSUM,RSUM,PSUM ! sub-cloud layer T, R, P sums
556 REAL LDEPTH ! layer depth
557 !internal function variables
558 REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE
559 REAL dum, icefactor, adiabatic_factor
561 REAL sec, secdel ! time step, increment in seconds
562 INTEGER i, j, k, IFOUT, ind(1)
563 CHARACTER*256 :: message
566 !secdel = 0.05
567 secdel = 5.0
568 g=9.81
569 r_d = 287.
571 !!!!!!!!!!!!!!!! 1. UPDRAFT PROPERTIES !!!!!!!!!!!!!!!!!!!!!!!
572 !!! DEFINE UPDRAFT PROPERTIES !!!
573 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
574 ! Upper limit of simulation in seconds
575 TAU = 7200. !simulation ends
577 !Set initial updraft strength - you could reduce to simulate the embryo
578 ! hovering around the edges of the updraft, as in Heymsfield and Musil (1982)
579 DO i=1,nz
580 VUU_pert(i) = VUU(i) * 1.
581 ENDDO
584 ! Initialize diameters to 0.
585 DO i=1,5
586 dhails(i) = 0.
587 ENDDO
589 ! Cap updraft lifetime at 2000 sec.
590 IF (wdur .GT. 2000) THEN
591 RTIME = 2000.
592 ELSE
593 RTIME = wdur
594 ENDIF
597 !!!!!!!!!!!!!!!! 2. INITIAL EMBRYO !!!!!!!!!!!!!!!!!!!!!!!!!!!
598 !!! FIND PROPERTIES OF INITIAL EMBRYO LOCATION !!!
599 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
600 !Find the cloud base for end-of-algorithm purposes.
601 KBAS=nz
602 !KFZL=nz
603 DO k=1,nz
604 cwitot = qi1d(k) + qc1d(k)
605 !No longer include graupel in in-cloud ice amounts
606 !RIA(k) = qi1d(k) + qs1d(k) + qg1d(k)
607 RIA(k) = qi1d(k) + qs1d(k)
608 !RWA(k) = qc1d(k) + qr1d(k)
609 RWA(k) = qc1d(k)
610 !IF ((RIA(k) .ge. 0.0001) .and. (TCA(k).lt.260.155) .and. &
611 ! (k .lt. KFZL)) THEN
612 ! KFZL = k
613 !ENDIF
614 IF ((cwitot .ge. 1.E-12) .and. (k .lt. KBAS)) THEN
615 KBAS = k
616 ENDIF
617 ENDDO
618 !QC - our embryo can't start below the cloud base.
619 !IF (KFZL .lt. KBAS) THEN
620 ! KFZL = KBAS
621 !ENDIF
623 !Pull heights, etc. of these levels out of 1-d arrays.
624 ZBAS = h1d(KBAS)
625 TBAS = TCA(KBAS)
626 WBASP = PA(KBAS)
627 RBAS = RA(KBAS)
630 !At coarser resolutions WRF underestimates liquid water aloft.
631 !A fairer estimate is of the adiabatic liquid water profile, or
632 !the difference between the saturation mixing ratio at each level
633 !and the cloud base water vapor mixing ratio
634 DO k=1,nz
635 RWA_new(k) = RWA(k)
636 IF (k.LT.KBAS) THEN
637 RWA_adiabat(k) = 0.
638 CYCLE
639 ENDIF
640 !saturation vapor pressure
641 !ESVA(k) = 6.112*exp((2.453E6/461)*(1./273. - 1./TCA(k))) !hPa
642 ESVA(k) = 611.2*exp(17.67*(TCA(k)-273.155)/(TCA(k)-29.655)) !Pa
643 !saturation vapor mixing ratio
644 RSA(k) = 0.62197 * ESVA(k) / (PA(k) - ESVA(k))
645 !Up above -31, start converting to ice, entirely by -38
646 !(Rosenfeld and Woodley 2000)
647 IF (TCA(k).gt.242.155) THEN
648 icefactor = 1.
649 ELSE IF ((TCA(k).LE.242.155).AND.(TCA(k).GT.235.155)) THEN
650 icefactor = (1-(242.155-TCA(k))/5.)
651 ELSE
652 icefactor = 0.
653 ENDIF
654 !Don't want any negative liquid water values
655 IF (RBAS.GT.RSA(k)) THEN
656 RWA_adiabat(k) = (RBAS - RSA(k))*icefactor
657 ELSE
658 RWA_adiabat(k) = RWA(k)
659 ENDIF
660 !Remove cloud liquid water outputted at previous lower levels
661 ! -- bug fix 170506 -- !
662 IF (k.eq.KBAS) THEN
663 RWA_new(k) = RWA_adiabat(k)
664 ELSE IF ((k.ge.KBAS+1).AND.(RWA_adiabat(k).ge.1.E-12)) THEN
665 RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1)
666 IF (RWA_new(k).LT.0) RWA_new(k) = 0.
667 ENDIF
668 ENDDO
669 !remove the height factor from RWA_new
670 DO k=KBAS+1,nz !bug fix, ensure index start 1 higher
671 RWA_new(k) = RWA_new(k) / (h1d(k)-h1d(k-1))
672 ENDDO
676 !!!!!!!!!!!!!!!! 3. INITIAL EMBRYO SIZE !!!!!!!!!!!!!!!!!!!!!
677 !!! SET CONSTANT RANGE OF INITIAL EMBRYO SIZES !!!
678 !!! START LOOP !!!
679 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
680 ! See Adams-Selin and Ziegler 2016 MWR for explanation of why
681 ! these sizes were picked.
682 !Run each initial embryo size perturbation
683 DO i=1,5
684 SELECT CASE (i)
685 CASE (1)
686 !Initial hail embryo diameter in m, at cloud base
687 DD = 5.E-3
688 tk_embryo = 265.155 !-8C
689 CASE (2)
690 DD = 7.5E-3
691 tk_embryo = 265.155 !-8C
692 CASE (3)
693 DD = 5.E-3
694 tk_embryo = 260.155 !-13C
695 CASE (4)
696 DD = 7.5E-3
697 tk_embryo = 260.155 !-13C
698 CASE (5)
699 tk_embryo = 260.155 !-13C
700 DD = 1.E-2
701 END SELECT
704 RTIME = 2000.
705 IF (wdur .LT. RTIME) RTIME = wdur
707 TFZL = tk_embryo
708 CALL INTERPP(PA, WFZLP, TCA, tk_embryo, IFOUT, nz)
709 CALL INTERP(h1d, ZFZL, WFZLP, IFOUT, PA, nz)
710 CALL INTERP(RA, RFZL, WFZLP, IFOUT, PA, nz)
711 CALL INTERP(VUU_pert, VUFZL, WFZLP, IFOUT, PA, nz)
712 CALL INTERP(rho1d, DENSAFZL, WFZLP, IFOUT, PA, nz)
714 !Check if height of cloud base is above embryo insertion point
715 !RAS-20190423-!
716 IF (ZFZL < ZBAS-1000) GOTO 400
719 !Begin hail simulation time (seconds)
720 sec = 0.
723 !!!!!!!!!!!!!!!!!! 4. INITIAL PARAMETERS !!!!!!!!!!!!!!!!!
724 !!! PUT INITIAL PARAMETERS IN VARIABLES !!!
725 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
726 !Set initial values for parameters at freezing level
727 P = WFZLP
728 RS = RFZL
729 TC = TFZL
730 VU = VUFZL
731 Z = ZFZL - ht
732 LDEPTH = Z
733 DENSA = DENSAFZL
735 !Set initial hailstone parameters
736 nofroze=1 !Set test for embryo: 0 for never been frozen; 1 frozen
737 TS = TC
738 TSm1 = TS
739 TSm2 = TS
740 D = DD !hailstone diameter in m
741 FW = 0.0
742 DENSE = 500. !kg/m3
743 ITYPE=1. !Assume starts in dry growth.
744 CLOUDON=1 !we'll eventually turn cloud "off" once updraft past time limit
746 !Start time loop.
747 DO WHILE (sec .lt. TAU)
748 sec = sec + secdel
750 !!!!!!!!!!!!!!!!!! 5. CALCULATE PARAMETERS !!!!!!!!!!!!!!!!!
751 !!! CALCULATE UPDRAFT PROPERTIES !!!
752 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
753 !Intepolate vertical velocity to our new pressure level
754 CALL INTERP(VUU_pert,VUMAX,P,IFOUT,PA,nz)
755 !print *, 'INTERP VU: ', VU, P
757 !Outside pressure levels? If so, exit loop
758 IF (IFOUT.EQ.1) GOTO 100
760 !Sine wave multiplier on updraft strength
761 IF (SEC .GT. 0.0 .AND. SEC .LT. RTIME) THEN
762 VUCORE = VUMAX * (0.5*SIN((3.14159*SEC)/(RTIME))+0.5)*1.2
763 VU = VUCORE
764 ELSEIF (SEC .GE. RTIME) THEN
765 VU = 0.0
766 CLOUDON = 0
767 ENDIF
769 !Calculate terminal velocity of the hailstone
770 ! (use previous values)
771 CALL TERMINL(DENSA,DENSE,D,VT,TC)
773 !Actual velocity of hailstone (upwards positive)
774 V = VU - VT
776 !Use hydrostatic eq'n to calc height of next level
777 P = P - DENSA*g*V*secdel
778 Z = Z + V*secdel
780 !Interpolate cloud temp, qvapor at new p-level
781 CALL INTERP(TCA,TC,P,IFOUT,PA,nz)
782 CALL INTERP(RA,RS,P,IFOUT,PA,nz)
784 !New density of in-cloud air
785 DENSA=P/(r_d*(1.+0.609*RS/(1.+RS))*TC)
787 !Interpolate liquid, frozen water mix ratio at new level
788 CALL INTERP(RIA,RI,P,IFOUT,PA,nz)
789 CALL INTERP(RWA_new,RW,P,IFOUT,PA,nz)
790 XI = RI * DENSA * CLOUDON
791 XW = RW * DENSA * CLOUDON
792 IF( (XW+XI).GT.0) THEN
793 PC = XI / (XW+XI)
794 ELSE
795 PC = 1.
796 ENDIF
799 ! SATURATION VAPOUR DENSITY DIFFERENCE BETWTEEN STONE AND CLOUD
800 CALL VAPORCLOSE(DELRW,PC,TS,TC,ITYPE)
803 !!!!!!!!!!!!!!!!!! 6. STONE'S MASS GROWTH !!!!!!!!!!!!!!!!!!!!
804 CALL MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
805 TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,secdel,ITYPE,DELRW)
808 !!!!!!!!!!!!!!!!!! 7. HEAT BUDGET OF HAILSTONE !!!!!!!!!!!!!!!
809 CALL HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, &
810 DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,secdel,ITYPE,P)
813 !!!!! 8. TEST DIAMETER OF STONE AND HEIGHT ABOVE GROUND !!!!!!!
814 !!! TEST IF DIAMETER OF STONE IS GREATER THAN CRITICAL LIMIT, IF SO
815 !!! BREAK UP
816 WATER=FW*GM !KG
817 ! CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S SURFACE
818 CRIT = 2.0E-4
819 IF (WATER.GT.CRIT)THEN
820 CALL BREAKUP(DENSE,D,GM,FW,CRIT)
821 ENDIF
823 !!! Has stone reached below cloud base?
824 !IF (Z .LE. 0) GOTO 200
825 IF (Z .LE. ZBAS) GOTO 200
827 !calculate ice-only diameter size
828 D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333
830 !Has the stone entirely melted and it's below the freezing level?
831 IF ((D_ICE .LT. 1.E-8) .AND. (TC.GT.273.155)) GOTO 300
833 !move values to previous timestep value
834 TSm1 = TS
835 TSm2 = TSm1
837 ENDDO !end cloud lifetime loop
839 100 CONTINUE !outside pressure levels in model
840 200 CONTINUE !stone reached surface
841 300 CONTINUE !stone has entirely melted and is below freezing level
843 !!!!!!!!!!!!!!!!!! 9. MELT STONE BELOW CLOUD !!!!!!!!!!!!!!!!!!!!
844 !Did the stone shoot out the top of the storm?
845 !Then let's assume it's lost in the murky "outside storm" world.
846 IF (P.lt.PA(nz)) THEN
847 D=0.0
848 !Is the stone entirely water? Then set D=0 and exit.
849 ELSE IF(ABS(FW - 1.0).LT.0.001) THEN
850 D=0.0
851 ELSE IF (Z.GT.0) THEN
852 !If still frozen, then use melt routine to melt below cloud
853 ! based on mean below-cloud conditions.
855 !Calculate mean sub-cloud layer conditions
856 TSUM = 0.
857 RSUM = 0.
858 PSUM = 0.
859 DO k=1,KBAS
860 TSUM = TSUM + TCA(k)
861 PSUM = PSUM + PA(k)
862 RSUM = RSUM + RA(k)
863 ENDDO
864 TLAYER = TSUM / KBAS
865 PLAYER = PSUM / KBAS
866 RLAYER = RSUM / KBAS
868 !MELT is expecting a hailstone of only ice. At the surface
869 !we're only interested in the actual ice diameter of the hailstone,
870 !so let's shed any excess water now.
871 D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333
872 D = D_ICE
873 CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT,TS,DENSE)
875 ENDIF !end check for melting call
877 !RAS - 20190423 - !
878 400 CONTINUE !Check for if embryo insertion point is below cloud base
880 !Require embryo to have stayed aloft for at least some time
881 IF (sec.LT.60) D = 0.
883 !Check to make sure something hasn't gone horribly wrong
884 IF (D.GT.0.254) D = 0. !just consider missing for now if > 10 in
886 !assign hail size in mm for output
887 dhails(i) = D * 1000
889 ENDDO !end embryo size loop
892 dhail1 = dhails(1)
893 dhail2 = dhails(2)
894 dhail3 = dhails(3)
895 dhail4 = dhails(4)
896 dhail5 = dhails(5)
898 END SUBROUTINE hailstone_driver
901 SUBROUTINE INTERPP(PA,PVAL,TA,TVAL,IFOUT,ITEL)
902 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
903 !!!!
904 !!!! INTERP: to linearly interpolate value of pval at temperature tval
905 !!!! between two levels of pressure array pa and temperatures ta
906 !!!!
907 !!!! INPUT: PA 1D array of pressure, to be interpolated
908 !!!! TA 1D array of temperature
909 !!!! TVAL temperature value at which we want to calculate pressure
910 !!!! IFOUT set to 0 if TVAL outside range of TA
911 !!!! ITEL number of vertical levels
912 !!!! OUTPUT: PVAL interpolated pressure variable at temperature tval
913 !!!!
914 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
915 IMPLICIT NONE
917 REAL PVAL, TVAL
918 REAL, DIMENSION( ITEL) :: TA, PA
919 INTEGER ITEL, IFOUT
920 !local variables
921 INTEGER I
922 REAL FRACT
924 IFOUT=1
926 DO I=1,ITEL-1
927 IF ( (TVAL .LT. TA(I) .AND. TVAL .GE. TA(I+1)) .or. & ! dT/dz < 0
928 (TVAL .GT. TA(I) .AND. TVAL .LE. TA(I+1)) ) THEN ! dT/dz > 0
930 FRACT = (TA(I) - TVAL) / (TA(I) - TA(I+1))
931 !.... compute the pressure value pval at temperature tval
932 PVAL = ((1.0 - FRACT) * PA(I)) + (FRACT * PA(I+1))
934 !End loop
935 IFOUT=0
936 EXIT
937 ENDIF
938 ENDDO
940 END SUBROUTINE INTERPP
944 SUBROUTINE INTERP(AA,A,P,IFOUT,PA,ITEL)
945 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
946 !!!!
947 !!!! INTERP: to linearly interpolate values of A at level P
948 !!!! between two levels of AA (at levels PA)
949 !!!!
950 !!!! INPUT: AA 1D array of variable
951 !!!! PA 1D array of pressure
952 !!!! P new pressure level we want to calculate A at
953 !!!! IFOUT set to 0 if P outside range of PA
954 !!!! ITEL number of vertical levels
955 !!!! OUTPUT: A variable at pressure level P
956 !!!!
957 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
958 IMPLICIT NONE
960 REAL A, P
961 REAL, DIMENSION( ITEL) :: AA, PA
962 INTEGER ITEL, IFOUT
963 !local variables
964 INTEGER I
965 REAL PDIFF, VDIFF, RDIFF, VERH, ADIFF
967 IFOUT=1
969 DO I=1,ITEL-1
970 IF (P.LE.PA(I) .AND. P.GT.PA(I+1)) THEN
971 !Calculate ratio between vdiff and pdiff
972 PDIFF = PA(I)-PA(I+1)
973 VDIFF = PA(I)-P
974 VERH = VDIFF/PDIFF
976 !Calculate the difference between the 2 A values
977 RDIFF = AA(I+1) - AA(I)
979 !Calculate new value
980 A = AA(I) + RDIFF*VERH
982 !End loop
983 IFOUT=0
984 EXIT
985 ENDIF
986 ENDDO
988 END SUBROUTINE INTERP
991 SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC)
992 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
993 !!!!
994 !!!! INTERP: Calculate terminal velocity of the hailstone
995 !!!!
996 !!!! INPUT: DENSA density of updraft air (kg/m3)
997 !!!! DENSE density of hailstone
998 !!!! D diameter of hailstone (m)
999 !!!! TC updraft temperature (K)
1000 !!!! OUTPUT:VT hailstone terminal velocity (m/s)
1001 !!!!
1002 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1003 IMPLICIT NONE
1005 REAL*8 D
1006 REAL DENSA, DENSE, TC, VT
1007 REAL GMASS, GX, RE, W, Y
1008 REAL, PARAMETER :: PI = 3.141592654, G = 9.78956
1009 REAL ANU
1011 !Mass of stone in kg
1012 GMASS = (DENSE * PI * (D**3.)) / 6.
1014 !Dynamic viscosity
1015 ANU = (0.00001718)*(273.155+120.)/(TC+120.)*(TC/273.155)**(1.5)
1017 !CALC THE BEST NUMBER, X AND REYNOLDS NUMBER, RE
1018 GX=(8.0*GMASS*G*DENSA)/(PI*(ANU*ANU))
1019 RE=(GX/0.6)**0.5
1021 !SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON
1022 !THE BEST NUMBER. Follows RH87 pg. 2762, but fixes an error in their (B1).
1023 IF (GX.LT.550) THEN
1024 W=LOG10(GX)
1025 Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0)
1026 RE=10**Y
1027 VT=ANU*RE/(D*DENSA)
1028 ELSE IF (GX.GE.550.AND.GX.LT.1800) THEN
1029 W=LOG10(GX)
1030 Y= -1.81391 + 1.34671*W - 0.12427*(W**2.0) + 0.0063*(W**3.0)
1031 RE=10**Y
1032 VT=ANU*RE/(D*DENSA)
1033 ELSE IF (GX.GE.1800.AND.GX.LT.3.45E08) THEN
1034 RE=0.4487*(GX**0.5536)
1035 VT=ANU*RE/(D*DENSA)
1036 ELSE
1037 RE=(GX/0.6)**0.5
1038 VT=ANU*RE/(D*DENSA)
1039 ENDIF
1041 END SUBROUTINE TERMINL
1045 SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE)
1046 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1047 !!! VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY
1048 !!! BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD
1049 !!! AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT,
1050 !!! AND IF THE STONE IS IN WET OR DRY GROWTH REGIME
1051 !!!
1052 !!! INPUT: PC fraction of updraft water that is frozen
1053 !!! TS temperature of hailstone (K)
1054 !!! TC temperature of updraft air (K)
1055 !!! ITYPE wet (2) or dry (1) growth regime
1056 !!! OUTPUT: DELRW difference in sat vap. dens. between hail and air
1057 !!! (kg/m3)
1058 !!!
1059 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1061 IMPLICIT NONE
1062 REAL DELRW, PC, TS, TC
1063 INTEGER ITYPE
1064 !local variables
1065 REAL RV, ALV, ALS, RATIO
1066 DATA RV/461.48/,ALV/2500000./,ALS/2836050./
1067 REAL ESAT, RHOKOR, ESATW, RHOOMGW, ESATI, RHOOMGI, RHOOMG
1069 !!! FOR HAILSTONE: FIRST TEST IF STONE IS IN WET OR DRY GROWTH
1070 RATIO = 1./273.155
1071 IF(ITYPE.EQ.2) THEN !!WET GROWTH
1072 ESAT=611.*EXP(ALV/RV*(RATIO-1./TS))
1073 ELSE !!DRY GROWTH
1074 ESAT=611.*EXP(ALS/RV*(RATIO-1./TS))
1075 ENDIF
1076 RHOKOR=ESAT/(RV*TS)
1078 !!! NOW FOR THE AMBIENT/IN-CLOUD CONDITIONS
1079 ESATW=611.*EXP(ALV/RV*(RATIO-1./TC))
1080 RHOOMGW=ESATW/(RV*TC)
1081 ESATI=611.*EXP(ALS/RV*(RATIO-1./TC))
1082 RHOOMGI=ESATI/(RV*TC)
1083 !RHOOMG=PC*(RHOOMGI-RHOOMGW)+RHOOMGW
1084 RHOOMG = RHOOMGI !done as in hailtraj.f
1086 !!! CALC THE DIFFERENCE(KG/M3): <0 FOR CONDENSATION,
1087 !!! >0 FOR EVAPORATION
1088 DELRW=(RHOKOR-RHOOMG)
1090 END SUBROUTINE VAPORCLOSE
1094 SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
1095 TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,ITYPE,DELRW)
1096 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1097 !!! CALC THE STONE'S INCREASE IN MASS
1098 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1100 IMPLICIT NONE
1101 REAL*8 D
1102 REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DI,ANU,RE,AE, &
1103 TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,DELRW
1104 INTEGER ITYPE
1105 !local variables
1106 REAL PI, D0, GMW2, GMI2, EW, EI,DGMV
1107 REAL DENSEL, DENSELI, DENSELW
1108 REAL DC !MEAN CLOUD DROPLET DIAMETER (MICRONS, 1E-6M)
1109 REAL VOLL, VOLT !VOLUME OF NEW LAYER, TOTAL (M3)
1110 REAL VOL1, DGMW_NOSOAK, SOAK, SOAKM
1111 REAL DENSAC, E, AFACTOR, NS, TSCELSIUS, VIMP, W
1112 PI=3.141592654
1114 !!! CALCULATE THE DIFFUSIVITY DI (m2/s)
1115 D0=0.226*1.E-4 ! change to m2/s, not cm2/s
1116 DI=D0*(TC/273.155)**1.81*(100000./P)
1118 !!! COLLECTION EFFICIENCY FOR WATER AND ICE
1119 !EW=1.0
1120 !!!! IF TS WARMER THAN -5C THEN ACCRETE ALL THE ICE (EI=1.0)
1121 !!!! OTHERWISE EI=0.21
1122 !IF(TS.GE.268.15)THEN
1123 ! EI=1.00
1124 !ELSE
1125 ! EI=0.21
1126 !ENDIF
1128 !!! COLLECTION EFFICIENCY FOR WATER AND ICE
1129 EW=1.0
1130 !!! Linear function for ice accretion efficiency
1131 IF (TC .GE. 273.155) THEN
1132 EI=1.00
1133 ELSE IF (TC.GE.233.155) THEN
1134 EI= 1.0 - ( (273.155 - TS) / 40. )
1135 ELSE !cooler than -40C
1136 EI = 0.0
1137 ENDIF
1139 !!! CALCULATE THE VENTILATION COEFFICIENT - NEEDED FOR GROWTH FROM VAPOR
1140 !The coefficients in the ventilation coefficient equations have been
1141 !experimentally derived, and are expecting cal-C-g units. Do some conversions.
1142 DENSAC = DENSA * (1.E3) * (1.E-6)
1143 !dynamic viscosity
1144 ANU=1.717E-5*(393.0/(TC+120.0))*(TC/273.155)**1.5 !RAS units fix to agree with Roger and Yau ps. 102
1145 !!! CALCULATE THE REYNOLDS NUMBER - unitless
1146 RE=D*VT*DENSAC/ANU
1147 E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv)
1148 !!! SELECT APPROPRIATE VALUES OF AE ACCORDING TO RE
1149 IF(RE.LT.6000.0)THEN
1150 AE=0.78+0.308*E
1151 ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN
1152 !RAS 190713
1153 AE=0.5*0.76*E !bug fix to match RasmussenHeymsfield1987
1154 ELSEIF(RE.GE.20000.0) THEN
1155 AE = 0.5*(0.57+9.0E-6*RE)*E !bug fix to match RasmussenHeymsfield1987
1156 ENDIF
1159 !!! CALC HAILSTONE'S MASS (GM), MASS OF WATER (GMW) AND THE
1160 !!! MASS OF ICE IN THE STONE (GMI)
1161 GM=PI/6.*(D**3.)*DENSE
1162 GMW=FW*GM
1163 GMI=GM-GMW
1165 !!! STORE THE MASS
1166 GM1=GM
1168 !!! NEW MASS GROWTH CALCULATIONS WITH VARIABLE RIME
1169 !!! LAYER DENSITY BASED ON ZIEGLER ET AL. (1983)
1171 !!! CALCULATE INCREASE IN MASS DUE INTERCEPTED CLD WATER, USE
1172 !!! ORIGINAL DIAMETER
1173 GMW2=GMW+SEKDEL*(PI/4.*D**2.*VT*XW*EW)
1174 DGMW=GMW2-GMW
1175 GMW=GMW2
1177 !!! CALCULATE THE INCREASE IN MASS DUE INTERCEPTED CLOUD ICE
1178 GMI2=GMI+SEKDEL*(PI/4.*D**2.*VT*XI*EI)
1179 DGMI=GMI2-GMI
1180 GMI=GMI2
1182 !!! CALCULATE INCREASE IN MASS DUE TO SUBLIMATION/CONDENSATION OF
1183 !!! WATER VAPOR
1184 DGMV = SEKDEL*2*PI*D*AE*DI*DELRW
1185 IF (DGMV .LT. 0) DGMV=0
1187 !!! CALCULATE THE TOTAL MASS CHANGE
1188 DGM=DGMW+DGMI+DGMV
1189 !Dummy arguments in the event of no water, vapor, or ice growth
1190 DENSELW = 900.
1191 DENSELI = 700.
1192 !!! CALCULATE DENSITY OF NEW LAYER, DEPENDS ON FW AND ITYPE
1193 IF (ITYPE.EQ.1) THEN !DRY GROWTH
1194 !If hailstone encountered supercooled water, calculate new layer density
1195 ! using Macklin form
1196 IF ((DGMW.GT.0).OR.(DGMV.GT.0)) THEN
1197 !MEAN CLOUD DROPLET RADIUS, ASSUME CLOUD DROPLET CONC OF 3E8 M-3 (300 CM-3)
1198 DC = (0.74*XW / (3.14159*1000.*3.E8))**0.33333333 * 1.E6 !MICRONS
1199 !!! FIND THE STOKES NUMBER (rasmussen heymsfield 1985)
1200 NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50) !need hail radius in cm
1201 !!! FIND IMPACT VELOCITY (rasmussen heymsfield 1985)
1202 !W = LOG10(NS) - bug fix 20221216 - RAS !
1203 IF (NS.LT.0.1)THEN
1204 W=-1.
1205 ELSE
1206 W = LOG10(NS)
1207 ENDIF
1208 IF (RE.GT.200) THEN
1209 IF (NS.LT.0.1) THEN
1210 VIMP = 0.0
1211 ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN
1212 VIMP = (0.356 + 0.4738*W - 0.1233*W**2. &
1213 -0.1618*W**3. + 0.0807*W**4.)*VT
1214 ELSEIF (NS.GT.10) THEN
1215 VIMP = 0.63*VT
1216 ENDIF
1217 ELSEIF ((RE.GT.65).AND.(RE.LE.200)) THEN
1218 IF (NS.LT.0.1) THEN
1219 VIMP = 0.0
1220 ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN
1221 VIMP = (0.3272 + 0.4907*W - 0.09452*W**2. &
1222 -0.1906*W**3. + 0.07105*W**4.)*VT
1223 ELSEIF (NS.GT.10) THEN
1224 VIMP = 0.61*VT
1225 ENDIF
1226 ELSEIF ((RE.GT.20).AND.(RE.LE.65)) THEN
1227 IF (NS.LT.0.1) THEN
1228 VIMP = 0.0
1229 ELSEIF ((NS.GE.0.1).AND.(NS.LE.10)) THEN
1230 VIMP = (0.2927 + 0.5085*W - 0.03453*W**2. &
1231 -0.2184*W**3. + 0.03595*W**4.)*VT
1232 ELSEIF (NS.GT.10) THEN
1233 VIMP = 0.59*VT
1234 ENDIF
1235 ELSEIF (RE.LE.20) THEN
1236 IF (NS.LT.0.4) THEN
1237 VIMP = 0.0
1238 ELSEIF ((NS.GE.0.4).AND.(NS.LE.10)) THEN
1239 VIMP = (0.1701 + 0.7246*W + 0.2257*W**2. &
1240 -1.13*W**3. + 0.5756*W**4.)*VT
1241 ELSEIF (NS.GT.10) THEN
1242 VIMP = 0.57*VT
1243 ENDIF
1244 ENDIF
1246 !RIME LAYER DENSITY, HEYMSFIELD AND PFLAUM 1985 FORM
1247 TSCELSIUS = TS - 273.16
1248 AFACTOR = -DC*VIMP/TSCELSIUS
1249 IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN
1250 DENSELW = 0.30*(AFACTOR)**0.44
1251 ELSE
1252 DENSELW = EXP(-0.03115 - 1.7030*AFACTOR + &
1253 0.9116*AFACTOR**2. - 0.1224*AFACTOR**3.)
1254 ENDIF
1256 DENSELW = DENSELW * 1000. !KG M-3
1257 !BOUND POSSIBLE DENSITIES
1258 IF (DENSELW.LT.100) DENSELW=100
1259 IF (DENSELW.GT.900) DENSELW=900
1260 !WRITE(12,*) 'MASSAGR, PFLAUM, MACKLIN: ', DENSELW, &
1261 ! MACKLIN_DENSELW
1262 ENDIF
1263 IF (DGMI.GT.0) THEN
1264 !Ice collection main source of growth, so set new density layer
1265 DENSELI = 700.
1266 ENDIF
1268 !All liquid water contributes to growth, none is soaked into center.
1269 DGMW_NOSOAK = DGMW !All liquid water contributes to growth,
1270 ! none of it is soaked into center.
1272 ELSE !WET GROWTH
1273 !Collected liquid water can soak into the stone before freezing,
1274 ! increasing mass and density but leaving volume constant.
1275 !Volume of current drop, before growth
1276 VOL1 = GM/DENSE
1277 !Difference b/w mass of stone if density is 900 kg/m3, and
1278 ! current mass
1279 SOAK = 900*VOL1 - GM
1280 !Liquid mass available
1281 SOAKM = DGMW
1282 !Soak up as much liquid as we can, up to a density of 900 kg/m3
1283 IF (SOAKM.GT.SOAK) SOAKM=SOAK
1284 GM = GM+SOAKM !Mass of current drop, plus soaking
1285 !New density of current drop, including soaking but before growth
1286 DENSE = GM/VOL1
1287 !Mass increment of liquid water growth that doesn't
1288 ! include the liquid water we just soaked into the stone.
1289 DGMW_NOSOAK = DGMW - SOAKM
1291 !Whatever growth does occur has high density
1292 DENSELW = 900. !KG M-3
1293 DENSELI = 900.
1295 ENDIF
1297 !!!VOLUME OF NEW LAYER
1298 !VOLL = (DGM) / DENSEL
1299 !VOLL = (DGMI+DGMV+DGMW_NOSOAK) / DENSEL
1300 !VOLL = (DGMI) / DENSELI + (DGMW_NOSOAK+DGMV) / DENSELW
1301 IF (DGMI.LE.0) THEN
1302 VOLL = (DGMW_NOSOAK+DGMV) / DENSELW
1303 ELSE IF (DGMW.LE.0) THEN
1304 VOLL = (DGMI) / DENSELI
1305 ELSE
1306 VOLL = (DGMI) / DENSELI + (DGMW_NOSOAK+DGMV) / DENSELW
1307 ENDIF
1309 !!!NEW TOTAL VOLUME, DENSITY, DIAMETER
1310 VOLT = VOLL + GM/DENSE
1311 !DENSE = (GM+DGM) / VOLT
1312 DENSE = (GM+DGMI+DGMV+DGMW_NOSOAK) / VOLT
1313 !D=D+SEKDEL*0.5*VT/DENSE*(XW*EW+XI*EI)
1314 GM = GM+DGMI+DGMW_NOSOAK+DGMV
1315 D = ( (6*GM) / (PI*DENSE) )**0.33333333
1317 END SUBROUTINE MASSAGR
1321 SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, &
1322 DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,ITYPE,P)
1323 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1324 !!! CALCULATE HAILSTONE'S HEAT BUDGET
1325 !!! See Rasmussen and Heymsfield 1987; JAS
1326 !!! Original Hailcast's variable units
1327 !!! TS - Celsius
1328 !!! FW - unitless, between 0 and 1
1329 !!! TC - Celsius
1330 !!! VT - m/s
1331 !!! D - m
1332 !!! DELRW - g/cm3 (per comment)
1333 !!! DENSA - g/cm3 (per comment)
1334 !!! GM1, DMG, DGMW, DGMV, DGMI, GMW, GMI - should all be kg
1335 !!! DI - cm2 / sec
1336 !!! P - hPa
1337 !!! Original HAILCAST HEATBUD subroutine uses c-g-s units, so do some conversions
1338 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1339 IMPLICIT NONE
1340 REAL*8 D
1341 REAL TS,TSm1,TSm2,FW,TC,VT,DELRW,DENSA,GM1,GM,DGM,DGMW,DGMV, &
1342 DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,P
1343 INTEGER ITYPE
1345 REAL RV, RD, G, PI, ALF, ALV, ALS, CI, CW, AK
1346 REAL H, AH, TCC, TSC, DELRWC, DENSAC, TDIFF
1347 REAL DCM, GMC, DGMWC, DGMIC, GM1C, DGMC
1348 REAL DMLT
1349 REAL TSCm1, TSCm2
1350 DATA RV/461.48/,RD/287.04/,G/9.78956/
1351 DATA PI/3.141592654/,ALF/79.7/,ALV/597.3/
1352 DATA ALS/677.0/,CI/0.5/,CW/1./
1354 !Convert values to non-SI units here to match all the constants
1355 TSC = TS - 273.155
1356 TSCm1 = TSm1 - 273.155
1357 TSCm2 = TSm2 - 273.155
1358 TCC = TC - 273.155
1359 DELRWC = DELRW * (1.E3) * (1.E-6)
1360 DENSAC = DENSA * (1.E3) * (1.E-6)
1361 !DI still in cm2/sec
1362 DCM = D * 100. !m to cm
1363 GMC = GM * 1000. !kg to g
1364 GM1C = GM1 * 1000. !kg to g
1365 DGMWC = DGMW * 1000. !kg to g
1366 DGMIC = DGMI * 1000. !kg to g
1367 DGMC = DGM * 1000. !kg to g
1370 !!! CALCULATE THE CONSTANTS
1371 AK=(5.8+0.0184*TCC)*1.E-5 !thermal conductivity - cal/(cm*sec*K)
1372 !dynamic viscosity - calculated in MASSAGR
1373 !ANU=1.717E-5*(393.0/(TC+120.0))*(TC/273.155)**1.5
1375 !!! CALCULATE THE REYNOLDS NUMBER - unitless
1376 !RE=D*VT*DENSAC/ANU - calculated in MASSAGR
1378 H=(0.71)**(0.333333333)*(RE**0.50) !ventilation coefficient heat (fh)
1379 !E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv)
1381 !!! SELECT APPROPRIATE VALUES OF AH AND AE ACCORDING TO RE
1382 IF(RE.LT.6000.0)THEN
1383 AH=0.78+0.308*H
1384 !AE=0.78+0.308*E
1385 ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN
1386 !RAS 190713
1387 AH=0.5*0.76*H !bug fix to match RasmussenHeymsfield1987
1388 ELSEIF(RE.GE.20000.0) THEN
1389 !RAS 190713
1390 AH=0.5*(0.57+9.0E-6*RE)*H !bug fix to match RasmussenHeymsfield1987
1391 !AE=(0.57+9.0E-6*RE)*E calculated in MASSAGR
1392 ENDIF
1394 !!! FOR DRY GROWTH FW=0, CALCULATE NEW TS, ITIPE=1
1395 !!! FOR WET GROWTH TS=0, CALCULATE NEW FW, ITIPE=2
1398 IF(ITYPE.EQ.1) THEN
1399 !!! DRY GROWTH; CALC NEW TEMP OF THE STONE
1400 !Original Hailcast algorithm (no time differencing)
1401 !TSC=TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* &
1402 ! (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ &
1403 ! DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)
1404 !units fix
1405 TSC=0.6*(TSC-TSC*DGMC/GM1C+SEKDEL/(GM1C*CI)* &
1406 (2.*PI*DCM*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ &
1407 DGMWC/SEKDEL*(ALF+CW*TCC)+DGMIC/SEKDEL*CI*TCC)) + &
1408 0.2*TSCm1 + 0.2*TSCm2
1410 TS = TSC+273.155
1411 IF (TS.GE.273.155) THEN
1412 TS=273.155
1413 ENDIF
1414 TDIFF = ABS(TS-273.155)
1415 IF (TDIFF.LE.1.E-6) ITYPE=2 !NOW IN WET GROWTH
1417 ELSE IF (ITYPE.EQ.2) THEN
1418 !!! WET GROWTH; CALC NEW FW
1420 IF (TCC.LT.0.) THEN
1421 !Original Hailcast algorithm
1422 !FW=FW-FW*DGM/GM1+SEKDEL/(GM1*ALF)* &
1423 ! (2.*PI*D*(AH*AK*TCC-AE*ALV*DI*DELRWC)+ &
1424 ! DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)
1425 !units fixed
1426 FW = FW-FW*DGMC/GMC+SEKDEL/(GMC*ALF)* &
1427 (2.*PI*DCM*(AH*AK*TCC-AE*ALV*DI*DELRWC)+ &
1428 DGMWC/SEKDEL*(ALF+CW*TCC)+DGMIC/SEKDEL*CI*TCC)
1429 ELSE
1430 !Calculate decrease in ice mass due to melting
1431 !!Bug fix 28Jul2021 RAS - non-SI units
1432 !More non-SI units - fixed 20220307 RAS
1433 ! Extra SEKDEL in the calculation - 20221102 RAS
1434 DMLT = SEKDEL / ALF * (2.*PI*DCM*AH*AK*TCC + 2.*PI*DCM*AE*ALV*DI*DELRWC + &
1435 DGMWC/SEKDEL*CW*TCC)
1436 FW = (FW*GM1C + DMLT + DGMWC) / GMC
1437 ENDIF
1439 IF(FW.GT.1.)FW=1.
1440 IF(FW.LT.0.)FW=0.
1442 !IF ALL OUR ACCRETED WATER WAS FROZEN, WE ARE BACK IN DRY GROWTH
1443 IF(FW.LE.1.E-6) THEN
1444 ITYPE=1
1445 ENDIF
1447 ENDIF
1449 END SUBROUTINE HEATBUD
1453 SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT)
1454 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1455 !!! TEST IF AMOUNT OF WATER ON SURFACE EXCEEDS CRTICAL LIMIT-
1456 !!! IF SO INVOKE SHEDDING SCHEME
1457 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1459 IMPLICIT NONE
1460 REAL*8 D
1461 REAL DENSE, GM, FW, CRIT
1462 !local variables
1463 REAL WATER, GMI, WAT, PI
1464 DATA PI/3.141592654/
1466 WATER=FW*GM !KG
1467 !GMI=(GM-WATER) !KG
1468 !REMAIN = CRIT*GM
1470 ! CALC CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S
1471 ! SURFACE
1472 !CRIT=0.268+0.1389*GMI
1473 !CRIT=0.268*1.E-3 + 0.1389*1.E-3*GMI !mass now in kg instead of g
1474 !CRIT = 1.0E-10
1475 !CRIT - now passed from main subroutine
1477 WAT=WATER-CRIT
1478 GM=GM-WAT
1479 FW=(CRIT)/GM
1481 IF(FW.GT.1.0) FW=1.0
1482 IF(FW.LT.0.0) FW=0.0
1484 ! RECALCULATE DIAMETER AFTER SHEDDING
1485 ! Assume density remains the same
1486 D=(6.*GM/(PI*DENSE))**(0.333333333)
1487 END SUBROUTINE BREAKUP
1490 SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT,TS,DENSE)
1491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1492 !!! This is a spherical hail melting estimate based on the Goyer
1493 !!! et al. (1969) eqn (3). The depth of the warm layer, estimated
1494 !!! terminal velocity, and mean temperature of the warm layer are
1495 !!! used. DRB. 11/17/2003.
1496 !!!
1497 !!! Incorporated some known variables into the subroutine
1498 !!! to use instead of constants (VT, TS, DENSE). RAS 10/2/2019
1499 !!!
1500 !!! Note - this subroutine assumes melted water is immediately
1501 !!! shed. Possibly not the most accurate assumption.
1502 !!!
1503 !!! INPUT: TLAYER mean sub-cloud layer temperature (K)
1504 !!! PLAYER mean sub-cloud layer pressure (Pa)
1505 !!! RLAYER mean sub-cloud layer mixing ratio (kg/kg)
1506 !!! VT terminal velocity of stone (m/s)
1507 !!! OUTPUT: D diameter (m)
1508 !!!
1509 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1510 IMPLICIT NONE
1512 REAL*8 D
1513 REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT, TS, DENSE
1514 REAL eenv, delta, ewet, de, der, wetold, wetbulb, wetbulbk
1515 REAL tdclayer, tclayer, eps, b, hplayer
1516 REAL*8 a
1517 REAL sd, lt, ka, lf, lv, t0, dv, pi, rv, rhoice, &
1518 tres, re, delt, esenv, rhosenv, essfc, rhosfc, dsig, &
1519 dmdt, mass, massorg, newmass, gamma, r, rho
1520 INTEGER wcnt
1522 !Convert temp to Celsius, calculate dewpoint in celsius
1523 eps = 0.622
1524 tclayer = TLAYER - 273.155
1525 a = 2.53E11
1526 b = 5.42E3
1527 tdclayer = b / LOG(a*eps / (rlayer*player))
1528 hplayer = player / 100.
1530 !Calculate partial vapor pressure
1531 eenv = (player*rlayer) / (rlayer+eps)
1532 eenv = eenv / 100. !convert to mb
1534 !Estimate wet bulb temperature (C)
1535 gamma = 6.6E-4*player
1536 delta = (4098.0*eenv)/((tdclayer+237.7)*(tdclayer+237.7))
1537 wetbulb = ((gamma*tclayer)+(delta*tdclayer))/(gamma+delta)
1539 !Iterate to get exact wet bulb
1540 wcnt = 0
1541 DO WHILE (wcnt .lt. 11)
1542 ewet = 6.108*(exp((17.27*wetbulb)/(237.3 + wetbulb)))
1543 de = (0.0006355*hplayer*(tclayer-wetbulb))-(ewet-eenv)
1544 der= (ewet*(.0091379024 - (6106.396/(273.155+wetbulb)**2))) &
1545 - (0.0006355*hplayer)
1546 wetold = wetbulb
1547 wetbulb = wetbulb - de/der
1548 wcnt = wcnt + 1
1549 IF ((abs(wetbulb-wetold)/wetbulb.gt.0.0001)) THEN
1550 EXIT
1551 ENDIF
1552 ENDDO
1554 wetbulbk = wetbulb + 273.155 !convert to K
1555 ka = .02 ! thermal conductivity of air
1556 lf = 3.34e5 ! latent heat of melting/fusion
1557 lv = 2.5e6 ! latent heat of vaporization
1558 t0 = 273.155 ! temp of ice/water melting interface
1559 dv = 0.25e-4 ! diffusivity of water vapor (m2/s)
1560 pi = 3.1415927
1561 rv = 1004. - 287. ! gas constant for water vapor
1562 !rhoice = 917.0 ! density of ice (kg/m**3)
1563 rhoice = DENSE ! we know the stone density, let's use it
1564 r = D/2. ! radius of stone (m)
1566 !Compute residence time in warm layer
1567 tres = LDEPTH / VT
1569 !Calculate dmdt based on eqn (3) of Goyer et al. (1969)
1570 !Reynolds number...from pg 317 of Atmo Physics (Salby 1996)
1571 !Just use the density of air at 850 mb...close enough.
1572 !rho = 85000./(287.*TLAYER)
1573 rho = player/(287.*TLAYER) !we have the layer pressure, just use it
1574 !units fix - r is now in meters, not mm. Plus we need D, not r. RAS 20220308
1575 re = rho*r*2*VT/1.7e-5
1577 !Temperature difference between environment and hailstone surface
1578 !We know the stone surface temperature, let's use it
1579 delt = wetbulb - (TS - 273.155) !again, wet bulb is in C
1581 !Difference in vapor density of air stream and equil vapor
1582 !density at the sfc of the hailstone
1583 esenv = 610.8*(exp((17.27*wetbulb)/ &
1584 (237.3 + wetbulb))) ! es environment in Pa
1585 rhosenv = esenv/(rv*wetbulbk)
1586 essfc = 610.8*(exp((17.27*(t0-273.155))/ &
1587 (237.3 + (t0-273.155)))) ! es environment in Pa
1588 rhosfc = essfc/(rv*t0)
1589 dsig = rhosenv - rhosfc
1591 !Calculate new mass growth
1592 dmdt = (-1.7*pi*r*(re**0.5)/lf)*((ka*delt)+((lv-lf)*dv*dsig))
1593 IF (dmdt.gt.0.) dmdt = 0
1594 mass = dmdt*tres ! < 0
1596 !Find the new hailstone diameter
1597 massorg = 1.33333333*pi*r*r*r*rhoice
1598 newmass = massorg + mass
1599 if (newmass.lt.0.0) newmass = 0.0
1600 D = 2.*(0.75*newmass/(pi*rhoice))**0.333333333
1601 END SUBROUTINE MELT
1604 END MODULE module_diag_hailcast
1606 #endif