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