!Copyright (c) 2016-2017 by Atmospheric and Environmental Research, Inc. (AER). 
!Software designed and tested in collaboration with the National Severe Storms 
!Laboratory, Norman, Oklahoma. 
!Documentation provided in Adams-Selin, R. and C. Ziegler, 2016: Forecasting hail 
!using a one-dimensional hail growth model within WRF. Mon. Wea. Rev., 144, 4919-4939.

PROGRAM hailcast_trajectory
    
    !to compile on Mac OS X:
    ! gfortran -o hailcast_trajectory hailcast_trajectory.f90 -I/opt/local/include -L/opt/local/lib -lnetcdf -lnetcdff
    
    !Given an embyro size and location, use the trajectory model
    ! to move it around in the storm until it exits the domain or
    ! falls below the surface.
    !
    !Set initial embryo size and X, Y, Z location at line 155-158.  X, Y, Z
    ! origin is in lower left corner.
    !
    !Prints out to a file a range of hailstone/in-cloud characteristics at 
    ! every trajectory timestep.  Output file name set at line 159.
    ! 
    !Currently outputs:
    ! Code variables:
    !   sec, X, Y, Z, D*1000., D_ICE*1000., DENSE, TS, ITYPE, FW, VT, RI*1000., RW*1000., TC
    ! 
    ! Which translates to:
        ! time, X Y Z location (m), diameter (mm), ice-only diameter (mm),
        ! density (kg/m3), hailstone temperature (K), dry(1)/wet(2) growth mode,
        ! fraction of hailstone that is water, terminal velocity (m/s), 
        ! in-cloud ice mixing ratio at hailstone point (g/kg),
        ! in-cloud liquid water mixing ratio at hailstone point (g/kg),
        ! in-cloud temperature at hailstone point (K)
    ! 
    
    USE netcdf
     
    IMPLICIT NONE
     
    !variables to use when reading in netcdf file
    character(200) :: filename
    integer :: ncid
     
    !dimensions
    integer, dimension(nf90_max_var_dims) :: dimIDs
    character(20), allocatable :: dimnames(:)
    integer, allocatable       :: dimvals(:)
    integer :: ndims, dims, varid
    integer :: nx, ny, nz
    integer :: ISTAT
    
    !variables to read in 
    REAL, ALLOCATABLE :: u(:,:,:), v(:,:,:), w(:,:,:), &
                         qc(:,:,:), qr(:,:,:), qv(:,:,:), qs(:,:,:), qi(:,:,:), &
                         temp(:,:,:), pres(:,:,:), rho(:,:,:)
    REAL, ALLOCATABLE :: x(:), y(:), z(:)
    
    !variables to pass to hailstone subroutine
    REAL, ALLOCATABLE :: RIA(:,:,:), RWA(:,:,:)
    REAL size  !hail size (mm)
    INTEGER :: e_x, e_y, e_z  !embryo location (m)
    REAL :: hail_x, hail_y, hail_z
    
    !unit number for writing out data
    INTEGER :: funit
    CHARACTER(5) :: e_xstr, e_ystr, e_zstr
    
    !loopers
    INTEGER :: i, j, k, n
    
    
    !file unit to write out trajectory info
    funit=10
    
    
    !open netcdf file
    filename = '3D_composites_WKsndumax31.nc'
    call check( nf90_open(filename, NF90_NOWRITE, ncid) )
    
    !get the dimensions
    call check( nf90_inquire(ncid, nDimensions = dims))
    allocate(dimnames(dims))
    allocate(dimvals(dims))
    do n=1,dims
      call check( nf90_inquire_dimension(ncid, n, name = dimnames(n), len = dimvals(n)))
      if (trim(dimnames(n)) .eq. 'nx') then
          nx=dimvals(n)
      endif
      if (trim(dimnames(n)) .eq. 'ny') then
          ny=dimvals(n)
      endif
      if (trim(dimnames(n)) .eq. 'nz') then
          nz=dimvals(n)
      endif
    enddo
    
    !allocate arrays
    ALLOCATE(u(nx,ny,nz))
    ALLOCATE(v(nx,ny,nz))
    ALLOCATE(w(nx,ny,nz))
    ALLOCATE(qc(nx,ny,nz))
    ALLOCATE(qr(nx,ny,nz))
    ALLOCATE(qv(nx,ny,nz))
    ALLOCATE(qs(nx,ny,nz))
    ALLOCATE(qi(nx,ny,nz))
    ALLOCATE(temp(nx,ny,nz))
    ALLOCATE(pres(nx,ny,nz))
    ALLOCATE(rho(nx,ny,nz))
    ALLOCATE(x(nx))
    ALLOCATE(y(ny))
    ALLOCATE(z(nz))
    
    !read in data
    call check( nf90_inq_varid(ncid, 'comp_u', varid))
    call check( nf90_get_var(ncid, varid, u))
    call check( nf90_inq_varid(ncid, 'comp_v', varid))
    call check( nf90_get_var(ncid, varid, v))
    call check( nf90_inq_varid(ncid, 'comp_w', varid))
    call check( nf90_get_var(ncid, varid, w))
    call check( nf90_inq_varid(ncid, 'comp_qc', varid))
    call check( nf90_get_var(ncid, varid, qc))
    call check( nf90_inq_varid(ncid, 'comp_qr', varid))
    call check( nf90_get_var(ncid, varid, qr))
    call check( nf90_inq_varid(ncid, 'comp_qs', varid))
    call check( nf90_get_var(ncid, varid, qs))
    call check( nf90_inq_varid(ncid, 'comp_qi', varid))
    call check( nf90_get_var(ncid, varid, qi))
    call check( nf90_inq_varid(ncid, 'comp_qv', varid))
    call check( nf90_get_var(ncid, varid, qv))
    call check( nf90_inq_varid(ncid, 'comp_t', varid))
    call check( nf90_get_var(ncid, varid, temp))
    call check( nf90_inq_varid(ncid, 'comp_p', varid))
    call check( nf90_get_var(ncid, varid, pres))
    call check( nf90_inq_varid(ncid, 'comp_rho', varid))
    call check( nf90_get_var(ncid, varid, rho))
    call check( nf90_inq_varid(ncid, 'x', varid))
    call check( nf90_get_var(ncid, varid, x))
    call check( nf90_inq_varid(ncid, 'y', varid))
    call check( nf90_get_var(ncid, varid, y))
    call check( nf90_inq_varid(ncid, 'z', varid))
    call check( nf90_get_var(ncid, varid, z))
    ISTAT = NF90_CLOSE(ncid)
    
    
    !!Calculate total ice, liquid water content
    ALLOCATE(RIA(nx,ny,nz))
    ALLOCATE(RWA(nx,ny,nz))
    DO k=1,nz
    DO j=1,ny
    DO i=1,nx
       RIA(i,j,k) = qi(i,j,k) + qs(i,j,k)
       RWA(i,j,k) = qc(i,j,k)  ! didn't include rain
    ENDDO
    ENDDO
    ENDDO
        
    
    
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!!!   CHANGE YOUR EMBRYO LOCATION, SIZE HERE
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    size=7.5E-3  ! m  (recommend b/w 2-10 mm)
    hail_x = 16000 !m   e_x;  0 western edge
    hail_y = 19000 !m   e_y;  0 southern edge
    hail_z = 11000 !m   e_z;  0 surface
    OPEN (unit=funit, file='traj_16_19_11.0.output')
    CALL hailstone_trajectory(temp, pres, rho, qv, RIA, RWA, &
                              u, v, w, z, y, x, nz, ny, nx, &
                              hail_z, hail_y, hail_x, funit, size)

    CLOSE(funit)

CONTAINS

  SUBROUTINE hailstone_trajectory ( TCA, PA, rho, RA, RIA, RWA,  &
                                UA, VA, VUU, ZA, YA, XA, nz, ny, nx, &
                                Z, Y, X, funit, dhail )
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!! Hailstone driver, adapted from hailstone subroutine in HAILCAST
!!!!  Inputs:
!!!!    3-d (nx, ny, nz)
!!!!     TCA          temperature (K) 
!!!!     PA           total pressure (Pa)
!!!!     rho          density (kg/m3)
!!!!     RA           vapor mixing ratio (kg/kg)
!!!!     RWA          cloud water mixing ratio (kg/kg)
!!!!     RIA          cloud ice and snow mixing ratio (kg/kg)
!!!!     VUU          updraft speed at each level (m/s)
!!!!     UA           u wind (m/s)
!!!!     VA           v wind (m/s)
!!!!    1-d
!!!!     ZA          height above ground (m)
!!!!     XA, YA      distance in x, y direction, starts at 0 in ll corner (m)
!!!!    Float 
!!!!     dhail      embryo size (m)
!!!!     X, Y, Z    embryo location (m)
!!!!    Integer
!!!!     funit      output file unit
!!!!     nz, ny, nz size of XA, YA, ZA
!!!!
!!!!  Output:
!!!!     dhail      hail diameter in mm one hailstone either:
!!!!                1) falls below cloud base
!!!!                2) reaches surface
!!!!                3) is advected outside the domain
!!!!     X, Y, Z    hailstone location (m)
!!!!
!!!!  13 Aug 2013 .................................Becky Adams-Selin AER
!!!!     adapted from hailstone subroutine in SPC's HAILCAST
!!!!  18 Mar 2014 .................................Becky Adams-Selin AER
!!!!     added variable rime layer density, per Ziegler et al. (1983)
!!!!  4 Jun 2014 ..................................Becky Adams-Selin AER
!!!!     removed initial embryo size dependency on microphysic scheme
!!!!  5 Jun 2014 ..................................Becky Adams-Selin AER
!!!!     used smaller initial embryo sizes
!!!!  25 Jun 2015..................................Becky Adams-Selin AER
!!!!     Significant revamping.  Fixed units bug in HEATBUD that caused
!!!!     hailstone temperature instabilities.  Similar issue fixed in BREAKUP
!!!!     subroutine.  Removed graupel from ice content.  Changed initial
!!!!     embryo size and location to better match literature.  Added
!!!!     enhanced melting when hailstone collides with liquid water
!!!!     in regions above freezing.  Final diameter returned is ice diameter
!!!!     only. Added hailstone temperature averaging over previous timesteps 
!!!!     to decrease initial temperature instability at small embyro diameters.  
!!!!  3 Sep 2015...................................Becky Adams-Selin AER
!!!!    Insert embryos at -13C; interpret pressure and other variables to
!!!!    that exact temperature level.
!!!! 16 Nov 2015...................................Becky Adams-Selin AER
!!!!     Hailstone travels horizontally through updraft instead of being
!!!!     locked in the center.
!!!!  9 Jul 2016...................................Becky Adams-Selin AER
!!!!     Uses an adiabatic liquid cloud water profile generated from
!!!!     the saturation vapor pressure using the model temperature
!!!!     profile.
!!!! 15 Apr 2017...................................Becky Adams-Selin AER
!!!!     Converted to a trajectory model for use in 3D steady-state simulations.
!!!!     Uses cloud water directly from netcdf file, no adiabatic profile.
!!!!
!!!!    
!!!! See Adams-Selin and Ziegler 2016, MWR for further documentation.
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    IMPLICIT NONE
    INTEGER, INTENT(IN ) :: funit, nz, ny, nx

    REAL, DIMENSION( nx, ny, nz ),             &
         INTENT(IN   ) ::                                  TCA  & ! temperature (K)
                                              ,            rho  &
                                              ,             PA  & ! pressure (Pa)
                                              ,             RA  & ! vapor mixing ratio (kg/kg)
                                              ,            VUU  & ! updraft speed (m/s)
                                              ,            RIA  & ! kg/kg
                                              ,            RWA  & ! kg/kg
                                              ,             UA  & ! m/s
                                              ,             VA    ! m/s

    REAL, DIMENSION( nx ),             &
         INTENT(IN   ) ::                                    XA ! m
    REAL, DIMENSION( ny ),             &
         INTENT(IN   ) ::                                    YA ! m
    REAL, DIMENSION( nz ),             &
         INTENT(IN   ) ::                                    ZA ! m
    
    !Output: hail diameter, location returned
    REAL, INTENT(INOUT) ::                              dhail  ! hail diameter (mm);
    REAL, INTENT(INOUT) ::                            X, Y, Z  ! hail location (m)


    !Local variables
    REAL, DIMENSION( nx,ny ) ::  ZBASA  ! height of cloud base across domain
    REAL ZBAS                  ! cloud base in this grid column
    REAL g                     ! gravity (m/s)
    REAL r_d                   ! constant
    !hailstone parameters
    REAL*8 D, D_ICE            ! hail diameter (m)
    REAL VT                    ! terminal velocity (m/s)
    REAL V                     ! actual stone velocity (m/s)
    REAL TS                    ! hailstone temperature (K)
    !RAS15.10 HAILSTONE temperature differencing
    REAL TSm1, TSm2            ! hailstone temperature at previous 3 timesteps
    REAL FW                    ! fraction of stone that is liquid
    REAL WATER                 ! mass of stone that is liquid
    REAL CRIT                  ! critical water fraction before shedding
    REAL DENSE                 ! hailstone density (kg/m3)
    INTEGER ITYPE              ! wet (2) or dry (1) growth regime
    !in-cloud updraft parameters at location of hailstone
    REAL P                     ! in-cloud pressure (Pa)
    REAL RS                    ! in-cloud saturation mixing ratio 
    REAL RI, RW                ! ice, liquid water mix. ratio (kg/kg)
    REAL XI, XW                ! ice, liquid water content (kg/m3 air)
    REAL PC                    ! in-cloud fraction of frozen water
    REAL TC                    ! in-cloud temperature (K)
    REAL VU                    ! in-cloud updraft speed (m/s)
    REAL UC                    ! in-cloud u wind speed (m/s)
    REAL VC                    ! in-cloud v wind speed (m/s)
    REAL DENSA                 ! in-cloud updraft density (kg/m3)
    REAL DELRW                 ! diff in sat vap. dens. between hail and air (kg/m3)
                                     !used as initial hail embryo size
    !mean sub-cloud layer variables
    REAL TLAYER,RLAYER,PLAYER  ! mean sub-cloud temp, mix ratio, pres
    REAL TSUM,RSUM,PSUM        ! sub-cloud layer T, R, P sums
    REAL LDEPTH                ! layer depth
    !internal function variables
    REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE
    REAL dum, DIFF, icefactor, adiabatic_factor
      
    REAL sec, secdel           ! time step, increment in seconds
    INTEGER i, j, k, IFOUT, ind(1), inject, t
    CHARACTER*256 :: message

    !set some constants
    secdel = 1.0
    g=9.81
    r_d = 287.



    !given x, y, z position of embryo, interpolate state variables
    CALL INTERP3D(PA, P, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
    CALL INTERP3D(RA, RS, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
    !CALL INTERP3D(VUU, VU, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
    CALL INTERP3D(rho, DENSA, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
    CALL INTERP3D(TCA, TC, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)

    
    !Set initial hailstone parameters
    TS = TC
    TSm1 = TS  !temperature differencing to account for instabilities
    TSm2 = TS
    D = dhail
    FW = 0.0
    DENSE = 500.  !kg/m3
    ITYPE = 1  !Assume starts in dry growth
    
    !Begin hail simulation time (seconds)
    sec = 0.


    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !write out some interesting data to the output file
    WRITE(funit, '(8(a12,1x),a2,5(a12,1x))') 'sec', 'X', 'Y', 'Z', 'D(mm)', &
                    'D_ICE(mm)', 'DENSE', 'TS', &
                    'ITYPE', 'FW', 'VT', 'RI(g/kg)', 'RW(g/kg)', 'TC'
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    

    
    !Start advection loop
    DO WHILE ( (X.gt.XA(1)) .AND. (X.lt.XA(nx)) .AND. &
               (Y.gt.YA(1)) .AND. (Z.lt.YA(ny)) .AND. &
               (Z.gt.ZA(1)) .AND. (Z.lt.ZA(nz)) )
               

        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        !write out some interesting data to the output file
        WRITE(funit, '(8(f12.5,1x),i2,5(f12.5,1x))') sec, X, Y, Z, D*1000., &
                        D_ICE*1000., DENSE, TS, &
                        ITYPE, FW, VT, RI*1000., RW*1000., TC
        !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


        sec = sec + secdel

        !Intepolate vertical, horizontal velocity for our position
        CALL INTERP3D(VUU, VU, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        CALL INTERP3D(UA, UC, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        CALL INTERP3D(VA, VC, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        
        !Calculate terminal velocity of the hailstone 
        ! (use previous values)
        CALL TERMINL(DENSA,DENSE,D,VT,TC)

        !Actual velocity of hailstone (upwards positive)
        V = VU - VT

        !Calc new height, location through advection
        Z = Z + V*secdel
        X = X + UC*secdel
        Y = Y + VC*secdel

        !Interpolate cloud temp, qvapor, pressure at new location
        CALL INTERP3D(TCA, TC, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        CALL INTERP3D(RA, RS, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        CALL INTERP3D(PA, P, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)

        !Outside model domain?  If so, exit loop
        IF (IFOUT.EQ.1) GOTO 100

        !New density of in-cloud air
        DENSA=P/(r_d*(1.+0.609*RS/(1.+RS))*TC)

        !Interpolate liquid, frozen water mix ratio at new level
        CALL INTERP3D(RIA, RI, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        CALL INTERP3D(RWA, RW, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
        !Are we inside any cloud yet?  If not, skip to next timestep
        IF ( (RI+RW) .LT. 1.E-12) THEN
            CYCLE
        ENDIF
        
        !Calculate ice/liquid water content
        XI = RI * DENSA
        XW = RW * DENSA 
        IF( (XW+XI).GT.0) THEN
        PC = XI / (XW+XI)
        ELSE
        PC = 1.
        ENDIF
        
        
        ! SATURATION VAPOUR DENSITY DIFFERENCE BETWTEEN STONE AND CLOUD
        CALL VAPORCLOSE(DELRW,PC,TS,TC,ITYPE)
        WRITE(*,*), sec, D*1000.
      

        !!!!!!!!!!!!!!!!!!!!   STONE'S MASS GROWTH !!!!!!!!!!!!!!!!!!!!
        CALL MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,&
                 TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,secdel,ITYPE,DELRW) 


        !!!!!!!!!!!!!!!!!!!!   HEAT BUDGET OF HAILSTONE !!!!!!!!!!!!!!!
        CALL HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW,  & 
                     DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,secdel,ITYPE,P)


 
        !!!!!!!! TEST DIAMETER OF STONE AND HEIGHT ABOVE GROUND !!!!!!!
        !!!  TEST IF DIAMETER OF STONE IS GREATER THAN CRITICAL LIMIT, IF SO  
        !!!  BREAK UP 
        WATER=FW*GM  !KG
        ! CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S SURFACE 
        CRIT = 2.0E-4
        IF (WATER.GT.CRIT)THEN
           CALL BREAKUP(DENSE,D,GM,FW,CRIT)
        ENDIF
        
                
        !calculate ice-only diameter size
        D_ICE = ( (6*GM*(1.-FW)) / (3.141592654*DENSE) )**0.33333333 


        !Has the stone entirely melted and it's below the freezing level?  
        IF ((D_ICE .LT. 1.E-8) .AND. (TC.GT.273.155)) GOTO 300

        !move values to previous timestep value
        TSm1 = TS
        TSm2 = TSm1
      ENDDO  ! end domain loop

100   CONTINUE !outside model domain (second in-loop check)
200   CONTINUE !stone reached below cloud base
300   CONTINUE !stone has entirely melted and is below freezing level


      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !write out some interesting data to the output file
      WRITE(funit, '(8(f12.5,1x),i2,5(f12.5,1x))') sec, X, Y, Z, D*1000., &
                    D_ICE*1000., DENSE, TS, &
                    ITYPE, FW, VT, RI*1000., RW*1000., TC
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


      !Right now this does not include melting below the cloud.  
      ! HAILCAST does melting
      ! by estimating all melting in subcloud at once using mean
      ! subcloud conditions to estimate wet bulb temperature. 
      ! Could include this option later.

  END SUBROUTINE hailstone_trajectory


  SUBROUTINE INTERP3D(DUMA, DUM, XA, YA, ZA, X, Y, Z, IFOUT, nx, ny, nz)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!! INTERP: to trilinearly interpolate value of dum at point x, y, z
!!!!   within 3d array of DUMA with dimension arrays XA, YA, ZA
!!!!   Borrowed from the ever-so-helpful wikipedia article on 
!!!!   trilinear interpolation
!!!!
!!!! INPUT: DUMA        3D state variable array, to be interpolated
!!!!        XA, YA, ZA  1D array of X, Y, Z distance (m)
!!!!        X, Y, Z     point at which we want to calculate value of DUMA
!!!!        IFOUT       set to 1 if X, Y, Z outside range of XA, YA, ZA
!!!!        nx, ny, nz  number of hor/vertical gridpoints
!!!! OUTPUT: DUM interpolated DUMA variable at point x, y, z
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    IMPLICIT NONE

    INTEGER :: nx, ny, nz, IFOUT
    REAL, DIMENSION( nx, ny, nz) :: DUMA
    REAL :: DUM, X, Y, Z
    REAL, DIMENSION( nx ) :: XA
    REAL, DIMENSION( nx ) :: YA
    REAL, DIMENSION( nx ) :: ZA
    !local variables
    INTEGER I, J, K  !loopers
    REAL xd, yd, zd
    REAL c00, c01, c10, c11, c0, c1


    IFOUT = 1

    DO I=1,nx-1
    IF ( X .GT. XA(I) .AND. X .LE. XA(I+1) ) THEN
        DO J=1,ny-1
        IF ( Y .GT. YA(J) .AND. Y .LE. YA(J+1) ) THEN
            DO K=1,nz-1
            IF ( Z .GT. ZA(K) .AND. Z .LE. ZA(K+1) ) THEN
                xd = (X - XA(I)) / (XA(I+1) - XA(I))
                yd = (Y - YA(J)) / (YA(J+1) - YA(J))
                zd = (Z - ZA(K)) / (ZA(K+1) - ZA(K))
                c00 = DUMA(  I,  J,  K)*(1.-xd) + DUMA(I+1,  J,  K)*xd
                c01 = DUMA(  I,  J,K+1)*(1.-xd) + DUMA(I+1,  J,K+1)*xd
                c10 = DUMA(  I,J+1,  K)*(1.-xd) + DUMA(I+1,J+1,  K)*xd
                c11 = DUMA(  I,I+1,K+1)*(1.-xd) + DUMA(I+1,J+1,K+1)*xd
                c0 = c00*(1.-yd) + c10*yd
                c1 = c01*(1.-yd) + c11*yd
                DUM = c0*(1.-zd) + c1*zd
                
                IFOUT = 0
                EXIT  !leave loop
            ENDIF
            ENDDO
        ENDIF
        ENDDO
    ENDIF
    ENDDO
    
  END SUBROUTINE INTERP3D
  
  
  SUBROUTINE INTERP2D(DUMA, DUM, XA, YA, X, Y, IFOUT, nx, ny)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!
!!!! INTERP: to bilinearly interpolate value of dum at point x, y
!!!!   within 2d array of DUMA with dimension arrays XA, YA
!!!!   Borrowed from the ever-so-helpful wikipedia article on 
!!!!   bilinear interpolation
!!!!
!!!! INPUT: DUMA        2D state variable array, to be interpolated
!!!!        XA, YA      1D array of X, Y distance (m)
!!!!        X, Y        point at which we want to calculate value of DUMA
!!!!        IFOUT       set to 1 if X, Y outside range of XA, YA 
!!!!        nx, ny      number of hor/vertical gridpoints
!!!! OUTPUT: DUM interpolated DUMA variable at point x, y
!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    IMPLICIT NONE

    INTEGER :: nx, ny, IFOUT
    REAL, DIMENSION( nx, ny) :: DUMA
    REAL :: DUM, X, Y
    REAL, DIMENSION( nx ) :: XA
    REAL, DIMENSION( nx ) :: YA
    !local variables
    INTEGER I, J, K  !loopers
    REAL fxy1, fxy2


    IFOUT = 1

    DO I=1,nx-1
    IF ( X .GT. XA(I) .AND. X .LE. XA(I+1) ) THEN
        DO J=1,ny-1
        IF ( Y .GT. YA(J) .AND. Y .LE. YA(J+1) ) THEN
            fxy1 = (XA(I+1) - X) / (XA(I+1) - XA(I)) * DUMA(I,J) + &
                   (X - XA(I))   / (XA(I+1) - XA(I)) * DUMA(I+1,J)
            fxy2 = (XA(I+1) - X) / (XA(I+1) - XA(I)) * DUMA(I,J+1) + &
                   (X - XA(I))   / (XA(I+1) - XA(I)) * DUMA(I+1,J+1)
            
            DUM = (YA(J+1) - Y) / (YA(J+1) - YA(J)) * fxy1 + &
                  (Y - YA(J))   / (YA(J+1) - YA(J)) * fxy2
            
            IFOUT = 0
            EXIT  !leave loop
        ENDIF
        ENDDO
    ENDIF
    ENDDO
    
  END SUBROUTINE INTERP2D
                

  SUBROUTINE TERMINL(DENSA,DENSE,D,VT,TC)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!
  !!!! INTERP: Calculate terminal velocity of the hailstone
  !!!!
  !!!! INPUT: DENSA  density of updraft air (kg/m3)
  !!!!        DENSE  density of hailstone
  !!!!        D      diameter of hailstone (m)
  !!!!        TC     updraft temperature (K)
  !!!! OUTPUT:VT     hailstone terminal velocity (m/s)
  !!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IMPLICIT NONE
      
      REAL*8 D
      REAL DENSA, DENSE, TC, VT
      REAL GMASS, GX, RE, W, Y
      REAL, PARAMETER :: PI = 3.141592654, G = 9.78956
      REAL ANU
      
      !Mass of stone in kg
      GMASS = (DENSE * PI * (D**3.)) / 6.
      
      !Dynamic viscosity
      ANU = (0.00001718)*(273.155+120.)/(TC+120.)*(TC/273.155)**(1.5)
      
      !CALC THE BEST NUMBER, X AND REYNOLDS NUMBER, RE 
      GX=(8.0*GMASS*G*DENSA)/(PI*(ANU*ANU))
      RE=(GX/0.6)**0.5

      !SELECT APPROPRIATE EQUATIONS FOR TERMINAL VELOCITY DEPENDING ON 
      !THE BEST NUMBER
      IF (GX.LT.550) THEN
        W=LOG10(GX)
        Y= -1.7095 + 1.33438*W - 0.11591*(W**2.0)      
        RE=10**Y
        VT=ANU*RE/(D*DENSA)
      ELSE IF (GX.GE.550.AND.GX.LT.1800) THEN
        W=LOG10(GX)
        Y= -1.81391 + 1.34671*W - 0.12427*(W**2.0) + 0.0063*(W**3.0)
        RE=10**Y
        VT=ANU*RE/(D*DENSA)
      ELSE IF (GX.GE.1800.AND.GX.LT.3.45E08) THEN
        RE=0.4487*(GX**0.5536)
        VT=ANU*RE/(D*DENSA)
      ELSE 
        RE=(GX/0.6)**0.5
        VT=ANU*RE/(D*DENSA)
      ENDIF
      
  END SUBROUTINE TERMINL   

   
   
  SUBROUTINE VAPORCLOSE(DELRW,PC,TS,TC,ITYPE)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!  VAPORCLOSE: CALC THE DIFFERENCE IN SATURATION VAPOUR DENSITY 
  !!!  BETWEEN THAT OVER THE HAILSTONE'S SURFACE AND THE IN-CLOUD 
  !!!  AIR, DEPENDS ON THE WATER/ICE RATIO OF THE UPDRAFT, 
  !!!  AND IF THE STONE IS IN WET OR DRY GROWTH REGIME
  !!!
  !!!  INPUT:  PC    fraction of updraft water that is frozen
  !!!          TS    temperature of hailstone (K)
  !!!          TC    temperature of updraft air (K)
  !!!          ITYPE wet (2) or dry (1) growth regime
  !!!  OUTPUT: DELRW difference in sat vap. dens. between hail and air
  !!!          (kg/m3)
  !!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      IMPLICIT NONE
      REAL DELRW, PC, TS, TC
      INTEGER ITYPE
      !local variables
      REAL RV, ALV, ALS, RATIO
      DATA RV/461.48/,ALV/2500000./,ALS/2836050./ 
      REAL ESAT, RHOKOR, ESATW, RHOOMGW, ESATI, RHOOMGI, RHOOMG

      !!!  FOR HAILSTONE:  FIRST TEST IF STONE IS IN WET OR DRY GROWTH
      RATIO = 1./273.155
      IF(ITYPE.EQ.2) THEN !!WET GROWTH
        ESAT=611.*EXP(ALV/RV*(RATIO-1./TS))
      ELSE  !!DRY GROWTH
        ESAT=611.*EXP(ALS/RV*(RATIO-1./TS))
      ENDIF
      RHOKOR=ESAT/(RV*TS)
      
      !!!  NOW FOR THE AMBIENT/IN-CLOUD CONDITIONS 
      ESATW=611.*EXP(ALV/RV*(RATIO-1./TC))
      RHOOMGW=ESATW/(RV*TC)
      ESATI=611.*EXP(ALS/RV*(RATIO-1./TC))
      RHOOMGI=ESATI/(RV*TC)
      !RHOOMG=PC*(RHOOMGI-RHOOMGW)+RHOOMGW
      RHOOMG = RHOOMGI  !done as in hailtraj.f

      !!!  CALC THE DIFFERENCE(KG/M3): <0 FOR CONDENSATION, 
      !!!  >0 FOR EVAPORATION
      DELRW=(RHOKOR-RHOOMG) 

  END SUBROUTINE VAPORCLOSE
     
      

  SUBROUTINE MASSAGR(D,GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DGMV,DI,ANU,RE,AE,& 
                 TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,ITYPE,DELRW) 
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!! CALC THE STONE'S INCREASE IN MASS 
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            
      IMPLICIT NONE
      REAL*8 D
      REAL GM,GM1,GMW,GMI,DGM,DGMW,DGMI,DI,ANU,RE,AE,  &
                 TC,TS,P,DENSE,DENSA,FW,VT,XW,XI,SEKDEL,DELRW
      INTEGER ITYPE 
      !local variables
      REAL PI, D0, GMW2, GMI2, EW, EI,DGMV
      REAL DENSEL, DENSELI, DENSELW 
      REAL DC !MEAN CLOUD DROPLET DIAMETER (MICRONS, 1E-6M)
      REAL VOLL, VOLT !VOLUME OF NEW LAYER, TOTAL (M3)
      REAL VOL1, DGMW_NOSOAK, SOAK, SOAKM
      REAL DENSAC, E, AFACTOR, NS, TSCELSIUS, VIMP, W
      PI=3.141592654

      !!!  CALCULATE THE DIFFUSIVITY DI (m2/s)
      D0=0.226*1.E-4  ! change to m2/s, not cm2/s
      DI=D0*(TC/273.155)**1.81*(100000./P)
  
      !!!  COLLECTION EFFICIENCY FOR WATER AND ICE 
      !EW=1.0      
      !!!!   IF TS WARMER THAN -5C THEN ACCRETE ALL THE ICE (EI=1.0) 
      !!!!   OTHERWISE EI=0.21      
      !IF(TS.GE.268.15)THEN
      !  EI=1.00
      !ELSE
      !  EI=0.21
      !ENDIF

      !!!  COLLECTION EFFICIENCY FOR WATER AND ICE 
      EW=1.0      
      !!!  Linear function for ice accretion efficiency
      IF (TC .GE. 273.155) THEN
         EI=1.00
      ELSE IF (TC.GE.233.155) THEN
         EI= 1.0 - ( (273.155 - TS) / 40. )
      ELSE  !cooler than -40C
         EI = 0.0
      ENDIF

      !!!  CALCULATE THE VENTILATION COEFFICIENT - NEEDED FOR GROWTH FROM VAPOR
      !The coefficients in the ventilation coefficient equations have been 
      !experimentally derived, and are expecting cal-C-g units.  Do some conversions.
      DENSAC = DENSA * (1.E3) * (1.E-6)
      !dynamic viscosity 
      ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5
      !!!  CALCULATE THE REYNOLDS NUMBER - unitless
      RE=D*VT*DENSAC/ANU   
      E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv)
      !!!   SELECT APPROPRIATE VALUES OF AE ACCORDING TO RE
      IF(RE.LT.6000.0)THEN
         AE=0.78+0.308*E
      ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN
         AE=0.76*E
      ELSEIF(RE.GE.20000.0) THEN
         AE=(0.57+9.0E-6*RE)*E
      ENDIF


      !!!  CALC HAILSTONE'S MASS (GM), MASS OF WATER (GMW) AND THE  
      !!!  MASS OF ICE IN THE STONE (GMI)
      GM=PI/6.*(D**3.)*DENSE
      GMW=FW*GM
      GMI=GM-GMW
  
      !!!  STORE THE MASS
      GM1=GM
      
      !!! NEW MASS GROWTH CALCULATIONS WITH VARIABLE RIME 
      !!! LAYER DENSITY BASED ON ZIEGLER ET AL. (1983)
      
      !!! CALCULATE INCREASE IN MASS DUE INTERCEPTED CLD WATER, USE
      !!! ORIGINAL DIAMETER
      GMW2=GMW+SEKDEL*(PI/4.*D**2.*VT*XW*EW)
      DGMW=GMW2-GMW 
      GMW=GMW2

      !!!  CALCULATE THE INCREASE IN MASS DUE INTERCEPTED CLOUD ICE
      GMI2=GMI+SEKDEL*(PI/4.*D**2.*VT*XI*EI)
      DGMI=GMI2-GMI 
      GMI=GMI2
  
      !!! CALCULATE INCREASE IN MASS DUE TO SUBLIMATION/CONDENSATION OF 
      !!! WATER VAPOR
      DGMV = SEKDEL*2*PI*D*AE*DI*DELRW
      IF (DGMV .LT. 0) DGMV=0

      !!!  CALCULATE THE TOTAL MASS CHANGE 
      DGM=DGMW+DGMI+DGMV
      !!! CALCULATE DENSITY OF NEW LAYER, DEPENDS ON FW AND ITYPE
      IF (ITYPE.EQ.1) THEN !DRY GROWTH
          !If hailstone encountered supercooled water, calculate new layer density 
          ! using Macklin form
          IF ((DGMW.GT.0).OR.(DGMV.GT.0)) THEN
             !MEAN CLOUD DROPLET RADIUS, ASSUME CLOUD DROPLET CONC OF 3E8 M-3 (300 CM-3)
             DC = (0.74*XW / (3.14159*1000.*3.E8))**0.33333333 * 1.E6 !MICRONS
             !!! FIND THE STOKES NUMBER  (rasmussen heymsfield 1985)
             NS = 2*VT*100.*(DC*1.E-4)**2. / (9*ANU*D*50)  !need hail radius in cm
             !!! FIND IMPACT VELOCITY (rasmussen heymsfield 1985)
             W = LOG10(NS)
             IF (RE.GT.200) THEN
                IF (NS.LT.0.1) THEN
                   VIMP = 0.0
                ELSEIF ((NS.GT.0.1).AND.(NS.LE.10)) THEN
                   VIMP = (0.356 + 0.4738*W - 0.1233*W**2. &
                           -0.1618*W**3. + 0.0807*W**4.)*VT
                ELSEIF (NS.GT.10) THEN
                   VIMP = 0.63*VT
                ENDIF
             ELSEIF ((RE.GT.65).AND.(RE.LE.200)) THEN
                IF (NS.LT.0.1) THEN
                   VIMP = 0.0
                ELSEIF ((NS.GT.0.1).AND.(NS.LE.10)) THEN
                   VIMP = (0.3272 + 0.4907*W - 0.09452*W**2. &
                           -0.1906*W**3. + 0.07105*W**4.)*VT
                ELSEIF (NS.GT.10) THEN
                   VIMP = 0.61*VT
                ENDIF
             ELSEIF ((RE.GT.20).AND.(RE.LE.65)) THEN
                IF (NS.LT.0.1) THEN
                   VIMP = 0.0
                ELSEIF ((NS.GT.0.1).AND.(NS.LE.10)) THEN
                   VIMP = (0.2927 + 0.5085*W - 0.03453*W**2. &
                           -0.2184*W**3. + 0.03595*W**4.)*VT
                ELSEIF (NS.GT.10) THEN
                   VIMP = 0.59*VT
                ENDIF
             ELSEIF (RE.LE.20) THEN
                IF (NS.LT.0.4) THEN
                   VIMP = 0.0
                ELSEIF ((NS.GT.0.4).AND.(NS.LE.10)) THEN
                   VIMP = (0.1701 + 0.7246*W + 0.2257*W**2. &
                           -1.13*W**3. + 0.5756*W**4.)*VT
                ELSEIF (NS.GT.10) THEN
                   VIMP = 0.57*VT
                ENDIF
             ENDIF
              
             !RIME LAYER DENSITY, HEYMSFIELD AND PFLAUM 1985 FORM
             TSCELSIUS = TS - 273.16 
             AFACTOR = -DC*VIMP/TSCELSIUS
             IF ((TSCELSIUS.LE.-5.).OR.(AFACTOR.GE.-1.60)) THEN
                 DENSELW = 0.30*(AFACTOR)**0.44
             ELSEIF (TSCELSIUS.GT.-5.) THEN
                 DENSELW = EXP(-0.03115 - 1.7030*AFACTOR + &
                               0.9116*AFACTOR**2. - 0.1224*AFACTOR**3.)
             ENDIF

             DENSELW = DENSELW * 1000. !KG M-3
             !BOUND POSSIBLE DENSITIES
             IF (DENSELW.LT.100) DENSELW=100
             IF (DENSELW.GT.900) DENSELW=900
             !WRITE(12,*) 'MASSAGR, PFLAUM, MACKLIN: ', DENSELW, &
             !             MACKLIN_DENSELW
          ENDIF
          IF (DGMI.GT.0) THEN
             !Ice collection main source of growth, so set new density layer
             DENSELI = 700.
          ENDIF
          
          !All liquid water contributes to growth, none is soaked into center.
          DGMW_NOSOAK = DGMW  !All liquid water contributes to growth,
                              ! none of it is soaked into center.

      ELSE !WET GROWTH
          !Collected liquid water can soak into the stone before freezing,
          ! increasing mass and density but leaving volume constant.
          !Volume of current drop, before growth 
          VOL1 = GM/DENSE
          !Difference b/w mass of stone if density is 900 kg/m3, and
          ! current mass
          SOAK = 900*VOL1 - GM
          !Liquid mass available
          SOAKM = DGMW
          !Soak up as much liquid as we can, up to a density of 900 kg/m3
          IF (SOAKM.GT.SOAK) SOAKM=SOAK
          GM = GM+SOAKM  !Mass of current drop, plus soaking
          !New density of current drop, including soaking but before growth
          DENSE = GM/VOL1 
          !Mass increment of liquid water growth that doesn't
          ! include the liquid water we just soaked into the stone.
          DGMW_NOSOAK = DGMW - SOAKM
          
          !Whatever growth does occur has high density
          DENSELW = 900.  !KG M-3
          DENSELI = 900.
         
      ENDIF

      !!!VOLUME OF NEW LAYER
      !VOLL = (DGM) / DENSEL
      !VOLL = (DGMI+DGMV+DGMW_NOSOAK) / DENSEL
      !VOLL = (DGMI) / DENSELI + (DGMW_NOSOAK+DGMV) / DENSELW
      IF (DGMI.LE.0) THEN
         VOLL = (DGMW_NOSOAK+DGMV) / DENSELW
      ELSE IF (DGMW.LE.0) THEN
         VOLL = (DGMI) / DENSELI
      ELSE
         VOLL = (DGMI) / DENSELI + (DGMW_NOSOAK+DGMV) / DENSELW
      ENDIF

      !!!NEW TOTAL VOLUME, DENSITY, DIAMETER
      VOLT = VOLL + GM/DENSE
      !DENSE = (GM+DGM) / VOLT
      DENSE = (GM+DGMI+DGMV+DGMW_NOSOAK) / VOLT
      !D=D+SEKDEL*0.5*VT/DENSE*(XW*EW+XI*EI)      
      GM = GM+DGMI+DGMW_NOSOAK+DGMV
      D = ( (6*GM) / (PI*DENSE) )**0.33333333 

  END SUBROUTINE MASSAGR



  SUBROUTINE HEATBUD(TS,TSm1,TSm2,FW,TC,VT,DELRW,D,DENSA,GM1,GM,DGM,DGMW,     &
                     DGMV,DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,ITYPE,P)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!! CALCULATE HAILSTONE'S HEAT BUDGET 
  !!! See Rasmussen and Heymsfield 1987; JAS
  !!! Original Hailcast's variable units
  !!! TS - Celsius
  !!! FW - unitless, between 0 and 1
  !!! TC - Celsius
  !!! VT - m/s
  !!! D  - m
  !!! DELRW - g/cm3 (per comment)
  !!! DENSA - g/cm3 (per comment)
  !!! GM1, DMG, DGMW, DGMV, DGMI, GMW, GMI - should all be kg
  !!! DI - cm2 / sec
  !!! P  - hPa
  !!! Original HAILCAST HEATBUD subroutine uses c-g-s units, so do some conversions
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IMPLICIT NONE
      REAL*8 D
      REAL TS,TSm1,TSm2,FW,TC,VT,DELRW,DENSA,GM1,GM,DGM,DGMW,DGMV,  &
                    DGMI,GMW,GMI,DI,ANU,RE,AE,SEKDEL,P
      INTEGER ITYPE
      
      REAL RV, RD, G, PI, ALF, ALV, ALS, CI, CW, AK
      REAL H, AH, TCC, TSC, DELRWC, DENSAC, TDIFF
      REAL DMLT
      REAL TSCm1, TSCm2
      DATA RV/461.48/,RD/287.04/,G/9.78956/
      DATA PI/3.141592654/,ALF/79.7/,ALV/597.3/
      DATA ALS/677.0/,CI/0.5/,CW/1./
      
      !Convert values to non-SI units here
      TSC = TS - 273.155
      TSCm1 = TSm1 - 273.155
      TSCm2 = TSm2 - 273.155
      TCC = TC - 273.155
      DELRWC = DELRW * (1.E3) * (1.E-6)
      DENSAC = DENSA * (1.E3) * (1.E-6)
      !DI still in cm2/sec


      !!!  CALCULATE THE CONSTANTS 
      AK=(5.8+0.0184*TCC)*1.E-5  !thermal conductivity - cal/(cm*sec*K)
      !dynamic viscosity - calculated in MASSAGR
      !ANU=1.717E-4*(393.0/(TC+120.0))*(TC/273.155)**1.5

      !!!  CALCULATE THE REYNOLDS NUMBER - unitless
      !RE=D*VT*DENSAC/ANU  - calculated in MASSAGR  
      
      H=(0.71)**(0.333333333)*(RE**0.50) !ventilation coefficient heat (fh)
      !E=(0.60)**(0.333333333)*(RE**0.50) !ventilation coefficient vapor (fv)

      !!!   SELECT APPROPRIATE VALUES OF AH AND AE ACCORDING TO RE
      IF(RE.LT.6000.0)THEN
         AH=0.78+0.308*H
         !AE=0.78+0.308*E
      ELSEIF(RE.GE.6000.0.AND.RE.LT.20000.0)THEN
         AH=0.76*H
         !AE=0.76*E
      ELSEIF(RE.GE.20000.0) THEN
         AH=(0.57+9.0E-6*RE)*H
         !AE=(0.57+9.0E-6*RE)*E  calculated in MASSAGR
      ENDIF

      !!!  FOR DRY GROWTH FW=0, CALCULATE NEW TS, ITIPE=1 
      !!!  FOR WET GROWTH TS=0, CALCULATE NEW FW, ITIPE=2


      IF(ITYPE.EQ.1) THEN
      !!!  DRY GROWTH; CALC NEW TEMP OF THE STONE 
         !Original Hailcast algorithm (no time differencing)
         !TSC=TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)*                &
         !   (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+     &
         !   DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)
         TSC=0.6*(TSC-TSC*DGM/GM1+SEKDEL/(GM1*CI)*                &
            (2.*PI*D*(AH*AK*(TCC-TSC)-AE*ALS*DI*DELRWC)+     &
            DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)) + &
            0.2*TSCm1 + 0.2*TSCm2
         
         TS = TSC+273.155
         IF (TS.GE.273.155) THEN 
            TS=273.155
         ENDIF
         TDIFF = ABS(TS-273.155)         
         IF (TDIFF.LE.1.E-6) ITYPE=2  !NOW IN WET GROWTH
     
      ELSE IF (ITYPE.EQ.2) THEN
      !!!  WET GROWTH; CALC NEW FW          
         
         IF (TCC.LT.0.) THEN
            !Original Hailcast algorithm
            FW=FW-FW*DGM/GM1+SEKDEL/(GM1*ALF)*               &
                (2.*PI*D*(AH*AK*TCC-AE*ALV*DI*DELRWC)+          &
                DGMW/SEKDEL*(ALF+CW*TCC)+DGMI/SEKDEL*CI*TCC)
         ELSE
            !Calculate decrease in ice mass due to melting
            DMLT = (2.*PI*D*AH*AK*TCC + 2.*PI*D*AE*ALV*DI*DELRWC + &
                    DGMW/SEKDEL*CW*TCC) / ALF
            FW = (FW*GM + DMLT) / GM
         ENDIF
         
         IF(FW.GT.1.)FW=1.
         IF(FW.LT.0.)FW=0.

         !IF ALL OUR ACCRETED WATER WAS FROZEN, WE ARE BACK IN DRY GROWTH
         IF(FW.LE.1.E-6) THEN
            ITYPE=1  
         ENDIF
         
      ENDIF

  END SUBROUTINE HEATBUD


  
  SUBROUTINE BREAKUP(DENSE,D,GM,FW,CRIT)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!  TEST IF AMOUNT OF WATER ON SURFACE EXCEEDS CRTICAL LIMIT- 
  !!!  IF SO INVOKE SHEDDING SCHEME 
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

      IMPLICIT NONE
      REAL*8 D
      REAL DENSE, GM, FW, CRIT
      !local variables
      REAL WATER, GMI, WAT, PI
      DATA PI/3.141592654/

      WATER=FW*GM  !KG
      !GMI=(GM-WATER) !KG
      !REMAIN = CRIT*GM

      ! CALC CRTICAL MASS CAPABLE OF BEING "SUPPORTED" ON THE STONE'S 
      ! SURFACE 
      !CRIT=0.268+0.1389*GMI 
      !CRIT=0.268*1.E-3 + 0.1389*1.E-3*GMI  !mass now in kg instead of g
      !CRIT = 1.0E-10
      !CRIT - now passed from main subroutine

      WAT=WATER-CRIT
      GM=GM-WAT
      FW=(CRIT)/GM
    
      IF(FW.GT.1.0) FW=1.0
      IF(FW.LT.0.0) FW=0.0

      ! RECALCULATE DIAMETER AFTER SHEDDING 
      ! Assume density remains the same
      D=(6.*GM/(PI*DENSE))**(0.333333333)
  END SUBROUTINE BREAKUP
  
  
  
  SUBROUTINE MELT(D,TLAYER,PLAYER,RLAYER,LDEPTH,VT)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!  This is a spherical hail melting estimate based on the Goyer 
  !!!  et al. (1969) eqn (3).  The depth of the warm layer, estimated 
  !!!  terminal velocity, and mean temperature of the warm layer are 
  !!!  used.  DRB.  11/17/2003.
  !!!
  !!!  INPUT:  TLAYER   mean sub-cloud layer temperature (K)
  !!!          PLAYER   mean sub-cloud layer pressure (Pa)
  !!!          RLAYER   mean sub-cloud layer mixing ratio (kg/kg)
  !!!          VT       terminal velocity of stone (m/s)
  !!!  OUTPUT: D        diameter (m)
  !!!          
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      IMPLICIT NONE

      REAL*8 D
      REAL TLAYER, PLAYER, RLAYER, LDEPTH, VT
      REAL eenv, delta, ewet, de, der, wetold, wetbulb, wetbulbk
      REAL tdclayer, tclayer, eps, b, hplayer
      REAL*8 a
      REAL sd, lt, ka, lf, lv, t0, dv, pi, rv, rhoice, &
           tres, re, delt, esenv, rhosenv, essfc, rhosfc, dsig, &
           dmdt, mass, massorg, newmass, gamma, r, rho
      INTEGER wcnt
      
      !Convert temp to Celsius, calculate dewpoint in celsius
      tclayer = TLAYER - 273.155
      a = 2.53E11
      b = 5.42E3
      tdclayer = b / LOG(a*eps / (rlayer*player))
      hplayer = player / 100.
      
      !Calculate partial vapor pressure
      eps = 0.622
      eenv = (player*rlayer) / (rlayer+eps)
      eenv = eenv / 100.  !convert to mb
      
      !Estimate wet bulb temperature (C)
      gamma = 6.6E-4*player
      delta = (4098.0*eenv)/((tdclayer+237.7)*(tdclayer+237.7))
      wetbulb = ((gamma*tclayer)+(delta*tdclayer))/(gamma+delta)
      
      !Iterate to get exact wet bulb
      wcnt = 0
      DO WHILE (wcnt .lt. 11)
        ewet = 6.108*(exp((17.27*wetbulb)/(237.3 + wetbulb))) 
        de = (0.0006355*hplayer*(tclayer-wetbulb))-(ewet-eenv)
        der= (ewet*(.0091379024 - (6106.396/(273.155+wetbulb)**2))) &
             - (0.0006355*hplayer)
        wetold = wetbulb
        wetbulb = wetbulb - de/der
        wcnt = wcnt + 1
        IF ((abs(wetbulb-wetold)/wetbulb.gt.0.0001)) THEN
           EXIT
        ENDIF
      ENDDO
      
      wetbulbk = wetbulb + 273.155  !convert to K
      ka = .02 ! thermal conductivity of air
      lf = 3.34e5 ! latent heat of melting/fusion
      lv = 2.5e6  ! latent heat of vaporization
      t0 = 273.155 ! temp of ice/water melting interface
      dv = 0.25e-4 ! diffusivity of water vapor (m2/s)
      pi = 3.1415927
      rv = 1004. - 287. ! gas constant for water vapor
      rhoice = 917.0 ! density of ice (kg/m**3)
      r = D/2. ! radius of stone (m)
      
      !Compute residence time in warm layer
      tres = LDEPTH / VT
        
      !Calculate dmdt based on eqn (3) of Goyer et al. (1969)
      !Reynolds number...from pg 317 of Atmo Physics (Salby 1996)
      !Just use the density of air at 850 mb...close enough.
      rho = 85000./(287.*TLAYER)
      re = rho*r*VT*.01/1.7e-5
      
      !Temperature difference between environment and hailstone surface
      delt = wetbulb !- 0.0 !assume stone surface is at 0C
                            !wetbulb is in Celsius

      !Difference in vapor density of air stream and equil vapor
      !density at the sfc of the hailstone
      esenv = 610.8*(exp((17.27*wetbulb)/  &
               (237.3 + wetbulb))) ! es environment in Pa
      rhosenv = esenv/(rv*wetbulbk)
      essfc = 610.8*(exp((17.27*(t0-273.155))/  &
               (237.3 + (t0-273.155)))) ! es environment in Pa
      rhosfc = essfc/(rv*t0)
      dsig = rhosenv - rhosfc

      !Calculate new mass growth
      dmdt = (-1.7*pi*r*(re**0.5)/lf)*((ka*delt)+((lv-lf)*dv*dsig))
      IF (dmdt.gt.0.) dmdt = 0
      mass = dmdt*tres
      
      !Find the new hailstone diameter
      massorg = 1.33333333*pi*r*r*r*rhoice
      newmass = massorg + mass
      if (newmass.lt.0.0) newmass = 0.0
      D = 2.*(0.75*newmass/(pi*rhoice))**0.333333333
  END SUBROUTINE MELT



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! SUBROUTINE: CHECK
!
! PURPOSE
! =======
! Checks NETCDF error status and stops program if there is an error.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE check(status)

  integer, intent (in) :: status

  if(status /= nf90_noerr) then
     print *, trim(nf90_strerror(status))
     stop 2
  end if

  END SUBROUTINE check

END PROGRAM hailcast_trajectory