REM ****                 PROGRAM TRW                 ****
REM ****           STANLEY DAVID GEDZELMAN           ****
REM ****  STARTED 12 FEB 1998 COMPLETED 15 FEB 1988  ****
REM ****    THIS PROGRAM SIMULATES A THUNDERSTORM    ****

REM *** GRID SIZE, NUMBERS OF PARTICLES, AND DIMENSIONS
XMAX = 20: ZMAX = .75 * XMAX: ZTROP = ZMAX * .7: REM KM
NHM = 500: NT = 500: NDT = 1000
DIM M(NHM), X(NHM), Z(NHM), W(NHM): REM HYDROMETEOR PROPERTIES
DIM XT(NT), ZT(NT), WT(NT), TT(NT), QT(NT): REM TRACER PROPERTIES
DIM NTG(XMAX, ZMAX), NHG(XMAX, ZMAX): REM BOX (FLUID) PROPERTIES
DIM WG(XMAX, ZMAX), QG(XMAX, ZMAX), QHG(XMAX, ZMAX), DQH(XMAX, ZMAX)
DIM KZTMAX(XMAX), PPT(XMAX)

REM ** CONSTANTS AND GRAPHICS SCREEN
G = .00981: U = .005: RHY = .04: REM KM
GAMMA = -7: GAMMAS = -6
T0 = 293: T2 = T0 + 2 * GAMMA * ZTROP
DT = 8: REM TIME STEP
KFRZ = INT(263 - T0) / GAMMA + 1: KMLT = INT(273 - T0) / GAMMA
QCRIT = .002: KRATE = .1: NGRO = 5: DQDZ = .0025: REM INIT CONDENSATION RATE
SCREEN 12
WINDOW (0, 0)-(XMAX, ZMAX)
LOCATE 4, 1: PRINT "ACCUMULATED": LOCATE 5, 1: PRINT "RAINFALL "
LOCATE 23, 1: PRINT "INITIAL": LOCATE 24, 1: PRINT "UPDRAFT"
RAIN = 10: SNOW = 13: REM COLORS OF RAIN AND SNOW

REM *** INITIAL POSITION AND TEMPERATURE OF TRACERS
RANDOMIZE TIMER
FOR I = 1 TO NT
XT(I) = RND * XMAX
ZT(I) = RND: REM MOIST AIR IN LOWEST KM
TT(I) = T0 + GAMMA * ZT(I)
PSET (XT(I), ZT(I))
NEXT I

REM *** START THE STORM
FOR J = 1 TO NDT

REM *** MOVE TRACERS AND CONDENSE VAPOR
FOR I = 1 TO NT
PSET (XT(I), ZT(I)), 0
XT(I) = XT(I) + U * DT
DZ = WT(I) * DT
TT(I) = TT(I) + DZ * GAMMAS
ZT(I) = ZT(I) + DZ
IF ZT(I) < ZTROP THEN
T = T0 + ZT(I) * GAMMA
ELSE
T = T2 - ZT(I) * GAMMA
END IF
QT(I) = QT(I) + DQDZ * DZ * (ZMAX - ZT(I)) / ZMAX
IF QT(I) < 0 THEN QT(I) = 0
BUOY = G * ((TT(I) - T) / TT(I) - QT(I))
WT(I) = WT(I) + DT * BUOY
IF XT(I) > XMAX OR ZT(I) < 0 THEN
XT(I) = RND
ZT(I) = RND
TT(I) = T0 + GAMMA * ZT(I)
WT(I) = .004 + .001 * RND
QT(I) = 0
END IF
PSET (XT(I), ZT(I))
IXT = INT(XT(I)) + 1: KZT = INT(ZT(I)) + 1
NTG(IXT, KZT) = NTG(IXT, KZT) + 1
QG(IXT, KZT) = QG(IXT, KZT) + QT(I): REM THESE MUST BE DIVIDED BY NTG LATER
WG(IXT, KZT) = WG(IXT, KZT) + WT(I)
IF KZTMAX(IXT) < KZT THEN KZTMAX(IXT) = KZT
NEXT I

FOR IG = 1 TO XMAX
FOR JG = 2 TO ZMAX
IF NTG(IG, JG) > 0 THEN
QG(IG, JG) = QG(IG, JG) / NTG(IG, JG)
WG(IG, JG) = WG(IG, JG) / NTG(IG, JG)
NTG(IG, JG) = 0
END IF
NEXT JG
NEXT IG

REM *** MOVE THE HYDROMETEORS
FOR I = 1 TO NH
JJ = I - NRES
X(JJ) = X(I): Z(JJ) = Z(I): M(JJ) = M(I)
CIRCLE (X(JJ), Z(JJ)), RHY, 0
X(JJ) = X(JJ) + U * DT
IF X(JJ) < XMAX THEN
IX = INT(X(JJ)) + 1: KZ = INT(Z(JJ)) + 1
        IF KZ > KMLT THEN
        WTERM = .003: HUE = SNOW
        ELSE
        WTERM = .006: HUE = RAIN
        END IF
Z(JJ) = Z(JJ) + (WG(IX, KZ) - WTERM) * DT
        IF Z(JJ) > 0 THEN
        CIRCLE (X(JJ), Z(JJ)), RHY, HUE
        KZ = INT(Z(JJ)) + 1
        QHG(IX, KZ) = QHG(IX, KZ) + M(JJ): REM ADD MASS OF EACH DROP IN BOX
        NHG(IX, KZ) = NHG(IX, KZ) + 1
        END IF
END IF
IF Z(JJ) < 0 THEN PPT(IX) = PPT(IX) + M(JJ)
IF X(JJ) > XMAX OR Z(JJ) < 0 THEN NRES = NRES + 1
NEXT I
NH = NH - NRES: NRES = 0

REM *** GROW HYDROMETEORS
FOR I = 1 TO NH
IX = INT(X(I)) + 1: KZ = INT(Z(I)) + 1
IF QG(IX, KZ) > 0 THEN
DQH(IX, KZ) = QG(IX, KZ) * (1 - EXP(-NHG(IX, KZ) / NGRO))
M(I) = (QHG(IX, KZ) + DQH(IX, KZ)) / NHG(IX, KZ)
END IF
NEXT I

REM *** ADJUST VERTICAL VELOCITY AND PRODUCE NEW PARTICLES
FOR IG = 1 TO XMAX
FOR JG = 2 TO ZMAX
WG(IG, JG) = WG(IG, JG) - G * QHG(IG, JG) * DT
IF QG(IG, JG) > 0 THEN
QG(IG, JG) = QG(IG, JG) - DQH(IG, JG)
        IF JG < KZTMAX(IG) AND QG(IG, JG) > QCRIT AND JG >= KFRZ THEN
        CTR2 = CTR2 + 1
        DM = (QG(IG, JG) - QCRIT) * KRATE: REM KRATE IS CREATION RATE
        Z(NH + CTR2) = RND + JG - 1
        X(NH + CTR2) = RND + IG - 1
        M(NH + CTR2) = DM
        IF JG > KMLT THEN HUE = SNOW ELSE HUE = RAIN
        CIRCLE (X(NH + CTR2), Z(NH + CTR2)), RHY, HUE
        QG(IG, JG) = QG(IG, JG) - DM
        END IF
END IF
NHG(IG, JG) = 0
QHG(IG, JG) = 0
DQH(IG, JG) = 0
NEXT JG
NEXT IG
NH = NH + CTR2
CTR2 = 0

REM *** AVERAGE VERTICAL VELOCITY OF AIR IN EACH BOX
FOR I = 1 TO NT
IG = INT(XT(I)) + 1: KG = INT(ZT(I)) + 1
IF KG > 1 THEN WT(I) = WG(IG, KG): QT(I) = QG(IG, KG)
NEXT I

FOR IG = 1 TO XMAX
FOR JG = 2 TO ZMAX: REM STARTING AT JG = 2 ALLOWS RANDOM VERTICAL MOTIONS IN THE LOWEST LAYER
WG(IG, JG) = 0: QG(IG, JG) = 0
NEXT JG
KZTMAX(IG) = 0
LINE (IG - .75, ZMAX - 1.5)-(IG - .25, ZMAX - 1.5 + PPT(IG)), 12, BF
NEXT IG

NEXT J