Skip to content

Commit 0eeab62

Browse files
adams-selinadams-selin
authored
HAILCAST bug fixes: tiling, non-SI units (wrf-model#1877)
TYPE: bug fix KEYWORDS: HAILCAST, tiling issues, undefined variables SOURCE: Becky Adams-Selin DESCRIPTION OF CHANGES: Problem: When running using tiling/quilted processes, the haildtacttime variable was being reset by every tile, thereby only activating HAILCAST to run on a single tile. Also, the RWA_new calculation can sometimes be undefined if the cloud base index is not set before the loop, and some SI/non-SI unit conversions were not handled correctly in the HEATBUD subroutine. Solution: The haildtacttime variable was changed to an gridded variable so could be set individually by gridpoint. Cloud base index definition added before the loop. SI/non-SI unit conversions fixed. LIST OF MODIFIED FILES: M Registry/Registry.EM_COMMON M phys/module_diag_hailcast.F TESTS CONDUCTED: Modifications fixed problem, and are currently used in operational HRRR v4.0. The Jenkins tests are all passing. RELEASE NOTE: HAILCAST quilted processes/tiling and a few other bug fixes.
1 parent e9b2182 commit 0eeab62

File tree

2 files changed

+113
-53
lines changed

2 files changed

+113
-53
lines changed

‎Registry/Registry.EM_COMMON

+1-1
Original file line numberDiff line numberDiff line change
@@ -3275,7 +3275,7 @@ state real HAILCAST_DIAM_MEAN ij misc 1 - -
32753275
state real HAILCAST_DIAM_STD ij misc 1 - - "HAILCAST_DIAM_STD" "WRF-HAILCAST Stand. Dev. Hail Diameter" "mm"
32763276
state real HAILCAST_WUP_MASK ij misc 1 - r "HAILCAST_WUP_MASK" "Updraft mask, 1 if > 10m/s" ""
32773277
state real HAILCAST_WDUR ij misc 1 - r "HAILCAST_WDUR" "Updraft duration" "sec"
3278-
state real haildtacttime - - - - r "haildtacttime" "HAILDTACTTIME" "HAILCAST ACTIVATION TIME in s"
3278+
state real haildtacttime ij misc 1 - r "haildtacttime" "HAILDTACTTIME" "HAILCAST ACTIVATION TIME in s"
32793279

32803280

32813281
package hailcast hailcast_opt==1 - state:hailcast_diam_max,hailcast_diam_mean,hailcast_diam_std,hailcast_dhail1,hailcast_dhail2,hailcast_dhail3,hailcast_dhail4,hailcast_dhail5

‎phys/module_diag_hailcast.F

+112-52
Original file line numberDiff line numberDiff line change
@@ -50,7 +50,8 @@ SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags &
5050

5151
REAL, INTENT(IN ),OPTIONAL :: haildt
5252
REAL, INTENT(IN ),OPTIONAL :: curr_secs
53-
REAL, INTENT(INOUT),OPTIONAL :: haildtacttime
53+
REAL, DIMENSION( ims:ime, jms:jme ), &
54+
INTENT(INOUT),OPTIONAL :: haildtacttime
5455
INTEGER, INTENT(IN ) :: itimestep
5556
REAL, INTENT(IN ) :: dt
5657

@@ -209,10 +210,18 @@ SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags &
209210
( (.not. grid%nested) .or. &
210211
( (grid%nested) .and. (abs(grid%dtbc) < 0.0001) ) ) )THEN
211212
doing_adapt_dt = .TRUE.
212-
IF ( haildtacttime .eq. 0. ) THEN
213-
haildtacttime = CURR_SECS + haildt
214-
END IF
215-
END IF
213+
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
216225
217226
! Calculate STEPHAIL
218227
stephail = NINT(haildt / dt)
@@ -270,12 +279,19 @@ SUBROUTINE hailcast_diagnostic_driver ( grid , config_flags &
270279
decided = .TRUE.
271280
END IF
272281
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
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
279295
END IF
280296
281297
write ( message, * ) 'timestep to run HAILCAST? ', run_param
@@ -442,6 +458,13 @@ END SUBROUTINE hailcast_diagnostic_driver
442458
!!!! 06 May 2017...................................Becky Adams-Selin AER
443459
!!!! Bug fixes for some memory issues in the adiabatic liquid
444460
!!!! 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.
445468
!!!!
446469
!!!! See Adams-Selin and Ziegler 2016, MWR for further documentation.
447470
!!!!
@@ -609,7 +632,8 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,&
609632
!the difference between the saturation mixing ratio at each level
610633
!and the cloud base water vapor mixing ratio
611634
DO k=1,nz
612-
IF (k.LT.KBAS) THEN
635+
RWA_new(k) = RWA(k)
636+
IF (k.LT.KBAS) THEN
613637
RWA_adiabat(k) = 0.
614638
CYCLE
615639
ENDIF
@@ -640,19 +664,10 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,&
640664
ELSE IF ((k.ge.KBAS+1).AND.(RWA_adiabat(k).ge.1.E-12)) THEN
641665
RWA_new(k) = RWA_adiabat(k)*(h1d(k)-h1d(k-1)) - RWA_new(k-1)
642666
IF (RWA_new(k).LT.0) RWA_new(k) = 0.
643-
ELSE
644-
RWA_new(k) = RWA(k)
645667
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
653668
ENDDO
654669
!remove the height factor from RWA_new
655-
DO k=KBAS,nz
670+
DO k=KBAS+1,nz !bug fix, ensure index start 1 higher
656671
RWA_new(k) = RWA_new(k) / (h1d(k)-h1d(k-1))
657672
ENDDO
658673
@@ -696,6 +711,11 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,&
696711
CALL INTERP(VUU_pert, VUFZL, WFZLP, IFOUT, PA, nz)
697712
CALL INTERP(rho1d, DENSAFZL, WFZLP, IFOUT, PA, nz)
698713
714+
!Check if height of cloud base is above embryo insertion point
715+
!RAS-20190423-!
716+
IF (ZFZL < ZBAS-1000) GOTO 400
717+
718+
699719
!Begin hail simulation time (seconds)
700720
sec = 0.
701721
@@ -850,9 +870,15 @@ SUBROUTINE hailstone_driver ( TCA, h1d, ht, PA, rho1d,&
850870
!so let's shed any excess water now.
851871
D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333
852872
D = D_ICE
853-
CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
873+
CALL MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT,TS,DENSE)
854874

855875
ENDIF !end check for melting call
876+
877+
!RAS - 20190423 - !
878+
400 CONTINUE !Check for if embryo insertion point is below cloud base
879+
880+
!Require embryo to have stayed aloft for at least some time
881+
IF (sec.LT.60) D = 0.
856882

857883
!Check to make sure something hasn't gone horribly wrong
858884
IF (D.GT.0.254) D = 0. !just consider missing for now if > 10 in
@@ -993,7 +1019,7 @@ SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC)
9931019
RE=(GX/0.6)**0.5
9941020
9951021
!SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON
996-
!THE BEST NUMBER
1022+
!THE BEST NUMBER. Follows RH87 pg. 2762, but fixes an error in their (B1).
9971023
IF (GX.LT.550) THEN
9981024
W=LOG10(GX)
9991025
Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0)
@@ -1115,17 +1141,18 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
11151141
!experimentally derived, and are expecting cal-C-g units. Do some conversions.
11161142
DENSAC = DENSA * (1.E3) * (1.E-6)
11171143
!dynamic viscosity
1118-
ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5
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
11191145
!!! CALCULATE THE REYNOLDS NUMBER - unitless
11201146
RE=D*VT*DENSAC/ANU
11211147
E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv)
11221148
!!! SELECT APPROPRIATE VALUES OF AE ACCORDING TO RE
11231149
IF(RE.LT.6000.0)THEN
11241150
AE=0.78+0.308*E
11251151
ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN
1126-
AE=0.76*E
1152+
!RAS 190713
1153+
AE=0.5*0.76*E !bug fix to match RasmussenHeymsfield1987
11271154
ELSEIF(RE.GE.20000.0) THEN
1128-
AE=(0.57+9.0E-6*RE)*E
1155+
AE = 0.5*(0.57+9.0E-6*RE)*E !bug fix to match RasmussenHeymsfield1987
11291156
ENDIF
11301157
11311158
@@ -1159,6 +1186,9 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
11591186

11601187
!!! CALCULATE THE TOTAL MASS CHANGE
11611188
DGM=DGMW+DGMI+DGMV
1189+
!Dummy arguments in the event of no water, vapor, or ice growth
1190+
DENSELW = 900.
1191+
DENSELI = 700.
11621192
!!! CALCULATE DENSITY OF NEW LAYER, DEPENDS ON FW AND ITYPE
11631193
IF (ITYPE.EQ.1) THEN !DRY GROWTH
11641194
!If hailstone encountered supercooled water, calculate new layer density
@@ -1169,7 +1199,12 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
11691199
!!! FIND THE STOKES NUMBER (rasmussen heymsfield 1985)
11701200
NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50) !need hail radius in cm
11711201
!!! FIND IMPACT VELOCITY (rasmussen heymsfield 1985)
1172-
W = LOG10(NS)
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
11731208
IF (RE.GT.200) THEN
11741209
IF (NS.LT.0.1) THEN
11751210
VIMP = 0.0
@@ -1213,7 +1248,7 @@ SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
12131248
AFACTOR = -DC*VIMP/TSCELSIUS
12141249
IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN
12151250
DENSELW = 0.30*(AFACTOR)**0.44
1216-
ELSEIF (TSCELSIUS.GT.-5.) THEN
1251+
ELSE
12171252
DENSELW = EXP(-0.03115 - 1.7030*AFACTOR + &
12181253
0.9116*AFACTOR**2. - 0.1224*AFACTOR**3.)
12191254
ENDIF
@@ -1309,26 +1344,33 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, &
13091344
13101345
REAL RV, RD, G, PI, ALF, ALV, ALS, CI, CW, AK
13111346
REAL H, AH, TCC, TSC, DELRWC, DENSAC, TDIFF
1347+
REAL DCM, GMC, DGMWC, DGMIC, GM1C, DGMC
13121348
REAL DMLT
13131349
REAL TSCm1, TSCm2
13141350
DATA RV/461.48/,RD/287.04/,G/9.78956/
13151351
DATA PI/3.141592654/,ALF/79.7/,ALV/597.3/
13161352
DATA ALS/677.0/,CI/0.5/,CW/1./
13171353
1318-
!Convert values to non-SI units here
1354+
!Convert values to non-SI units here to match all the constants
13191355
TSC = TS - 273.155
13201356
TSCm1 = TSm1 - 273.155
13211357
TSCm2 = TSm2 - 273.155
13221358
TCC = TC - 273.155
13231359
DELRWC = DELRW * (1.E3) * (1.E-6)
13241360
DENSAC = DENSA * (1.E3) * (1.E-6)
13251361
!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
13261368
13271369
13281370
!!! CALCULATE THE CONSTANTS
13291371
AK=(5.8+0.0184*TCC)*1.E-5 !thermal conductivity - cal/(cm*sec*K)
13301372
!dynamic viscosity - calculated in MASSAGR
1331-
!ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5
1373+
!ANU=1.717E-5*(393.0/(TC+120.0))*(TC/273.155)**1.5
13321374
13331375
!!! CALCULATE THE REYNOLDS NUMBER - unitless
13341376
!RE=D*VT*DENSAC/ANU - calculated in MASSAGR
@@ -1341,10 +1383,11 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, &
13411383
AH=0.78+0.308*H
13421384
!AE=0.78+0.308*E
13431385
ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN
1344-
AH=0.76*H
1345-
!AE=0.76*E
1386+
!RAS 190713
1387+
AH=0.5*0.76*H !bug fix to match RasmussenHeymsfield1987
13461388
ELSEIF(RE.GE.20000.0) THEN
1347-
AH=(0.57+9.0E-6*RE)*H
1389+
!RAS 190713
1390+
AH=0.5*(0.57+9.0E-6*RE)*H !bug fix to match RasmussenHeymsfield1987
13481391
!AE=(0.57+9.0E-6*RE)*E calculated in MASSAGR
13491392
ENDIF
13501393
@@ -1358,11 +1401,12 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, &
13581401
!TSC=TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)* &
13591402
! (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+ &
13601403
! 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)) + &
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)) + &
13641408
0.2*TSCm1 + 0.2*TSCm2
1365-
1409+
13661410
TS = TSC+273.155
13671411
IF (TS.GE.273.155) THEN
13681412
TS=273.155
@@ -1375,14 +1419,21 @@ SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW, &
13751419
13761420
IF (TCC.LT.0.) THEN
13771421
!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)
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)
13811429
ELSE
13821430
!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
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
13861437
ENDIF
13871438
13881439
IF(FW.GT.1.)FW=1.
@@ -1436,13 +1487,19 @@ SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT)
14361487
END SUBROUTINE BREAKUP
14371488

14381489

1439-
SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
1490+
SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT,TS,DENSE)
14401491
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
14411492
!!! This is a spherical hail melting estimate based on the Goyer
14421493
!!! et al. (1969) eqn (3). The depth of the warm layer, estimated
14431494
!!! terminal velocity, and mean temperature of the warm layer are
14441495
!!! used. DRB. 11/17/2003.
14451496
!!!
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+
!!!
14461503
!!! INPUT: TLAYER mean sub-cloud layer temperature (K)
14471504
!!! PLAYER mean sub-cloud layer pressure (Pa)
14481505
!!! RLAYER mean sub-cloud layer mixing ratio (kg/kg)
@@ -1453,7 +1510,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
14531510
IMPLICIT NONE
14541511

14551512
REAL*8 D
1456-
REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT
1513+
REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT, TS, DENSE
14571514
REAL eenv, delta, ewet, de, der, wetold, wetbulb, wetbulbk
14581515
REAL tdclayer, tclayer, eps, b, hplayer
14591516
REAL*8 a
@@ -1502,7 +1559,8 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
15021559
dv = 0.25e-4 ! diffusivity of water vapor (m2/s)
15031560
pi = 3.1415927
15041561
rv = 1004. - 287. ! gas constant for water vapor
1505-
rhoice = 917.0 ! density of ice (kg/m**3)
1562+
!rhoice = 917.0 ! density of ice (kg/m**3)
1563+
rhoice = DENSE ! we know the stone density, let's use it
15061564
r = D/2. ! radius of stone (m)
15071565
15081566
!Compute residence time in warm layer
@@ -1511,12 +1569,14 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
15111569
!Calculate dmdt based on eqn (3) of Goyer et al. (1969)
15121570
!Reynolds number...from pg 317 of Atmo Physics (Salby 1996)
15131571
!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
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
15161576
15171577
!Temperature difference between environment and hailstone surface
1518-
delt = wetbulb !- 0.0 !assume stone surface is at 0C
1519-
!wetbulb is in Celsius
1578+
!We know the stone surface temperature, let's use it
1579+
delt = wetbulb - (TS - 273.155) !again, wet bulb is in C
15201580

15211581
!Difference in vapor density of air stream and equil vapor
15221582
!density at the sfc of the hailstone
@@ -1531,7 +1591,7 @@ SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
15311591
!Calculate new mass growth
15321592
dmdt = (-1.7*pi*r*(re**0.5)/lf)*((ka*delt)+((lv-lf)*dv*dsig))
15331593
IF (dmdt.gt.0.) dmdt = 0
1534-
mass = dmdt*tres
1594+
mass = dmdt*tres ! < 0
15351595

15361596
!Find the new hailstone diameter
15371597
massorg = 1.33333333*pi*r*r*r*rhoice

0 commit comments

Comments
 (0)