!*******************************************************************************
!$RCSfile$
!$Revision$
!$Date$
!*******************************************************************************
SUBROUTINE set_evap_step()

USE scipuff_fi
USE surface_fi
USE srfevap_fi

!  Sets surface evaporation time step information

IMPLICIT NONE

INTEGER ilev
REAL    umax, qmax, slmax
REAL    xx, yy

!-----  Find maximum BL velocity for the large timestep

CALL SWIMgetUmax( umax,qmax,slmax )

xx = xmin + 0.5*FLOAT(nx-1)*dxg
yy = ymin + 0.5*FLOAT(ny-1)*dyg

CALL mapfac( xx,yy,xmap_evap,ymap_evap )

!-----  Set time levels

DO ilev = 0,MAXLEV_EVAP
  lev_evap(ilev) = -1
END DO

mxlev_evap  = -1
mxtlev_evap = -1

dx_evap = MIN(Psrfdep%dx/xmap_evap,Psrfdep%dy/ymap_evap)

CALL set_del_evap( umax,qmax,slmax,del_evap,facu_evap )
CALL get_lim_evap( del_evap,facu_evap )
IF( nError /= NO_ERROR )GOTO 9999

9999 CONTINUE

RETURN
END

!===============================================================================

SUBROUTINE set_del_evap( umax,qmax,slmax,del,facu )
USE scipuff_fi
USE surface_fi
USE srfevap_fi

IMPLICIT NONE

REAL, INTENT( IN  ) :: umax, qmax, slmax
REAL, INTENT( OUT ) :: del, facu

del = MAX(umax*xmap_evap/Psrfdep%dx,umax*ymap_evap/Psrfdep%dy)
IF( del > SMALL )THEN
  del = 1./del
ELSE
  del = delt
END IF

IF( qmax > SMALL )THEN
  facu = slmax**(0.333333)/qmax
ELSE
  facu = 1.E+20
END IF

RETURN
END

!===============================================================================

SUBROUTINE get_lim_evap( del,facu )

USE scipuff_fi
USE surface_fi
USE srfevap_fi
USE srfparam_fd

!  Gets maximum surface evaporation grid depth

IMPLICIT NONE

REAL, INTENT( IN ) :: del, facu

INTEGER i, ix, iy, mlevi, nxs, nys, ivm, ivv, ivs
INTEGER icell0

INTERFACE
  RECURSIVE SUBROUTINE limevap( srf,icell0,ivm,ivv,ivs,mlev0,dx0,dt0,facu )
    USE sagstr_fd
    TYPE( SAGgrid_str ), POINTER      :: srf
    INTEGER,             INTENT( IN ) :: ivm, ivv, ivs
    INTEGER,             INTENT( IN ) :: icell0, mlev0
    REAL,                INTENT( IN ) :: dx0, dt0, facu
  END SUBROUTINE limevap
END INTERFACE

nxs = Psrfdep%nx
nys = Psrfdep%ny

IF( ndep_blocks == 0 )GOTO 9999

DO i = 1,ndep_blocks

  IF( srf_block(i)%type == SBLK_EVAP )THEN
    ivm = (srf_block(i)%field-1)*Psrfdep%mxgrd
    ivv = ivm + Psrfdep%mxgrd
    ivs = ivv + Psrfdep%mxgrd

    DO ix = 1,nxs
      DO iy = 1,nys
        icell0 = (iy-1)*nxs + ix
        mlevi  = 0
        CALL limevap( Psrfdep,icell0,ivm,ivv,ivs,mlevi,dx_evap,del,facu )
        IF( nError /= NO_ERROR )GOTO 9999
      END DO
    END DO

  END IF

END DO

9999 CONTINUE

RETURN
END

!===============================================================================

RECURSIVE SUBROUTINE limevap( srf,icell0,ivm,ivv,ivs,mlev0,dx0,dt0,facu )

USE scipuff_fi
USE files_fi
USE srfevap_fi
USE sagstr_fd

IMPLICIT NONE

TYPE( SAGgrid_str ), POINTER      :: srf
INTEGER,             INTENT( IN ) :: ivm, ivv, ivs
INTEGER,             INTENT( IN ) :: icell0, mlev0
REAL,                INTENT( IN ) :: dx0, dt0, facu

INTEGER mlev, i, icell, ilev, jlev, ios
REAL    si, fac, dti, dx

INTEGER, EXTERNAL :: time_level

IF( srf%ipdat(icell0+ivm) > SMALL .AND. srf%ipdat(icell0+ivv) > SMALL )THEN
  ilev = MAX(mlev0-2,0)
  si   = srf%ipdat(icell0+ivs)/srf%ipdat(icell0+ivv)
  fac  = 0.5**ilev
  dx   = dx0*fac
  dti  = dt0*fac
  dti  = MIN(si**TwoThirds*facu,dti)

  jlev = time_level( dti )
  IF( nError /= NO_ERROR )THEN
    eRoutine = 'limevap'
    WRITE(lun_log,*,IOSTAT=ios)'******* TIME LEVEL ERROR ********'
    WRITE(lun_log,*,IOSTAT=ios)TRIM(eRoutine)
    WRITE(lun_log,*,IOSTAT=ios)TRIM(eInform)
    WRITE(lun_log,*,IOSTAT=ios)'DT=',dti
    WRITE(lun_log,*,IOSTAT=ios)'DELT=',delt
    WRITE(lun_log,*,IOSTAT=ios)'LEVEL=',jlev,MAXTLV,mlev0
    WRITE(lun_log,*,IOSTAT=ios)'dx,si,facu=',dx,si,facu
    WRITE(lun_log,*,IOSTAT=ios)'ivar,icell0=',ivm,icell0
    GOTO 9999
  END IF

  lev_evap(ilev) = MAX(lev_evap(ilev),jlev)
  mxlev_evap     = MAX(mxlev_evap,ilev)
  mxtlev_evap    = MAX(mxtlev_evap,lev_evap(ilev))
END IF

IF( srf%ipgrd(icell0) /= 0 )THEN
  DO i = 0,3
    icell = srf%ipgrd(icell0) + i
    mlev  = mlev0 + 1
    CALL limevap( srf,icell,ivm,ivv,ivs,mlev,dx0,dt0,facu )
  END DO
END IF

9999 CONTINUE

RETURN
END


!===============================================================================

SUBROUTINE step_srf_evap( lev,lev2 )

USE scipuff_fi
USE surface_fi
USE srfevap_fi
USE srfparam_fd

!  Steps surface deposition variables for secondary evaporation

IMPLICIT NONE

INTEGER, INTENT( IN    ) :: lev
INTEGER, INTENT( INOUT ) :: lev2

INTEGER imat, ilev, nxs, nys
INTEGER i, ix, iy, icell0, ivar
REAL    dxs, dys, xx, yy, dx, dy, xmins, ymins

TYPE( puff_material ) pmatlIn

INTERFACE
  RECURSIVE SUBROUTINE srfgrd_evap( srf,ivar,icell0,ilev,xx,yy,dx,dy,lev,lev2 )
    USE sagstr_fd
    TYPE( SAGgrid_str ), POINTER         :: srf    !SAG grid structure
    INTEGER,             INTENT( IN    ) :: ivar   !Variable offset
    INTEGER,             INTENT( IN    ) :: icell0 !Grid cell
    INTEGER,             INTENT( IN    ) :: ilev   !Current grid level
    INTEGER,             INTENT( IN    ) :: lev    !Current time level
    INTEGER,             INTENT( INOUT ) :: lev2   !Current max time level
    REAL,                INTENT( IN    ) :: xx,yy  !Cell centroid location
    REAL,                INTENT( IN    ) :: dx,dy  !Cell size (project coord)
  END SUBROUTINE srfgrd_evap
END INTERFACE

INTEGER, EXTERNAL :: num_puff_types

!------ Check time level for this time step

IF( lev > mxtlev_evap )RETURN

!------ Loop over grid

nxs   = Psrfdep%nx
nys   = Psrfdep%ny
xmins = Psrfdep%xmin
ymins = Psrfdep%ymin
dxs   = Psrfdep%dx
dys   = Psrfdep%dy

IF( nsrf_blocks == 0 )GOTO 9999

DO i = 1,nsrf_blocks

  IF( srf_block(i)%type == SBLK_EVAP )THEN

    imat = MOD(srf_block(i)%id,65536)
    ivar = (srf_block(i)%field-1)*Psrfdep%mxgrd

!---- Get liquid puff properties

    ityp_evap = material(imat)%ioffp + 2
    CALL get_puff_material( ityp_evap,pmatlIn )
    pmatl = TRANSFER(pmatlIn,pmatl)
    evap_min = material(imat)%prop(3)

!---- Reset type for vapor puff

    ityp_evap = ityp_evap - 1

    DO ix = 1,nxs
      DO iy = 1,nys
        icell0 = (iy-1)*nxs + ix
        ilev   = 0
        xx = xmins + (FLOAT(ix)-0.5)*dxs
        yy = ymins + (FLOAT(iy)-0.5)*dys
        CALL mapfac( xx,yy,xmap_evap,ymap_evap )
        dx = dxs
        dy = dys
        srf_met = .FALSE.
        rate    = 0.
        CALL srfgrd_evap( Psrfdep,ivar,icell0,ilev,xx,yy,dx,dy,lev,lev2 )
        IF( nError /= NO_ERROR )GOTO 9999
      END DO
    END DO

  END IF

END DO

9999 CONTINUE

RETURN
END

!===============================================================================

SUBROUTINE get_srf_met( xx,yy )

USE SWIMpuff_fd
USE SWIMparam_fd
USE srfevap_fi

!  Calculates surface meteorology needed for evaporation

IMPLICIT NONE

REAL, INTENT( IN ) :: xx, yy

TYPE( PuffMet ) :: Met

INTEGER irv

INTEGER, EXTERNAL :: SWIMcurrentSrfMet

irv = SWIMcurrentSrfMet( xx,yy,Met )
IF( irv /= SWIMsuccess )THEN
  CALL setSWIMerror( 'SWIMcurrentSrfMet' )
  GOTO 9999
END IF

zisrf   = Met%BL%MixingHt
hflxsrf = Met%BL%HeatFlux
tsrf    = Met%Mean%T
ubsrf   = Met%Mean%U
vbsrf   = Met%Mean%V
ustr    = Met%BL%UstarDep
z0      = Met%BL%zruf
srf_met = .TRUE.

9999 CONTINUE

RETURN
END

!===============================================================================

RECURSIVE SUBROUTINE srfgrd_evap( srf,ivar,icell0,ilev,xx,yy,dx,dy,lev,lev2 )

USE scipuff_fi
USE srfevap_fi
USE sagstr_fd

!  Calculates surface deposition from the adaptive grid

IMPLICIT NONE

TYPE( SAGgrid_str ), POINTER         :: srf    !SAG grid structure
INTEGER,             INTENT( IN    ) :: ivar   !Variable offset
INTEGER,             INTENT( IN    ) :: icell0 !Grid cell
INTEGER,             INTENT( IN    ) :: ilev   !Current grid level
INTEGER,             INTENT( IN    ) :: lev    !Current time level
INTEGER,             INTENT( INOUT ) :: lev2   !Current max time level
REAL,                INTENT( IN    ) :: xx,yy  !Cell centroid location
REAL,                INTENT( IN    ) :: dx,dy  !Cell size (project coord)

INTEGER ilevx, i, iy, icell
REAL    dtx, dxx, dyy, xxx, yyy

INTERFACE
  SUBROUTINE srf_evap( srf,ivar,icell0,ilev,xx,yy,dx,dy,dts,lev2 )
    USE sagstr_fd
    TYPE( SAGgrid_str ), POINTER         :: srf    !SAG grid structure
    INTEGER,             INTENT( IN    ) :: ivar   !Variable offset
    INTEGER,             INTENT( IN    ) :: icell0 !Grid cell
    INTEGER,             INTENT( IN    ) :: ilev   !Current grid level
    INTEGER,             INTENT( INOUT ) :: lev2   !Current max time level
    REAL,                INTENT( IN    ) :: xx,yy  !Cell centroid location
    REAL,                INTENT( IN    ) :: dx,dy  !Cell size (project coord)
    REAL,                INTENT( IN    ) :: dts    !Local time step
  END SUBROUTINE srf_evap
END INTERFACE

IF( lev_evap(ilev) < lev )THEN

!-------- Level below minimum, so continue to refine grid

  IF( srf%ipgrd(icell0) > 0 )THEN
    ilevx = ilev + 1
    dxx = 0.5*dx
    dyy = 0.5*dy
    DO i = 0,3
      icell = srf%ipgrd(icell0) + i
      iy = INT(i/2)
      xxx = xx + (FLOAT(i-iy-iy)-0.5)*dxx
      yyy = yy + (FLOAT(iy)-0.5)*dyy
      CALL srfgrd_evap( srf,ivar,icell,ilevx,xxx,yyy,dxx,dyy,lev,lev2 )
      IF( nError /= NO_ERROR )GOTO 9999
    END DO
  ELSE
    RETURN
  END IF

ELSE

  IF( srf%ipgrd(icell0) > 0 )THEN

    dtx = delt*(0.5**lev_evap(ilev))
    CALL srf_evap( srf,ivar,icell0,ilev,xx,yy,dx,dy,dtx,lev2 )
    IF( nError /= NO_ERROR )GOTO 9999

!---------- Check levels below if necessary

    IF( ilev < mxlev_evap .AND. found_subcell )THEN
      ilevx = ilev + 1
      dxx   = 0.5*dx
      dyy   = 0.5*dy
      DO i = 0,3
        icell = srf%ipgrd(icell0) + i
        iy = INT(i/2)
        xxx = xx + (FLOAT(i-iy-iy)-0.5)*dxx
        yyy = yy + (FLOAT(iy)-0.5)*dyy
        CALL srfgrd_evap( srf,ivar,icell,ilevx,xxx,yyy,dxx,dyy,lev,lev2 )
        IF( nError /= NO_ERROR )GOTO 9999
      END DO
    ELSE
      RETURN
    END IF
  ELSE
    RETURN
  END IF

END IF

9999 CONTINUE

RETURN
END

!===============================================================================

SUBROUTINE srf_evap( srf,ivar,icell0,ilev,xx,yy,dx,dy,dts,lev2 )

USE scipuff_fi
USE met_fi
USE srfevap_fi
USE sagstr_fd

!  Calculates surface deposition from the adaptive grid

IMPLICIT NONE

REAL, PARAMETER :: RT12 = 2.*SQRT3

TYPE( SAGgrid_str ), POINTER         :: srf    !SAG grid structure
INTEGER,             INTENT( IN    ) :: ivar   !Variable offset
INTEGER,             INTENT( IN    ) :: icell0 !Grid cell
INTEGER,             INTENT( IN    ) :: ilev   !Current grid level
INTEGER,             INTENT( INOUT ) :: lev2   !Current max time level
REAL,                INTENT( IN    ) :: xx,yy  !Cell centroid location
REAL,                INTENT( IN    ) :: dx,dy  !Cell size (project coord)
REAL,                INTENT( IN    ) :: dts    !Local time step

REAL, DIMENSION(7) :: msrf
REAL, DIMENSION(5) :: asrf

REAL    smass, dsrf, depth
REAL    mass, tau_evap
INTEGER iv, i, j, icell, jcell, i0, nv

INTEGER iy,jy
REAL    xxx, yyy, dxx, dyy
REAL    xxxx,yyyy,dxxx,dyyy

INTERFACE
  SUBROUTINE srf_cell_evap( srf,ivar,icell0,ilev )
    USE sagstr_fd
    TYPE( SAGgrid_str ), POINTER      :: srf    !Dep grid structure
    INTEGER,             INTENT( IN ) :: ivar   !Start of evap block
    INTEGER,             INTENT( IN ) :: ilev   !Grid refinement level
    INTEGER,             INTENT( IN ) :: icell0 !Upper cell location
  END SUBROUTINE srf_cell_evap
END INTERFACE

!------ Accumulate subcells for this level
!       msrf(1) = sum mean C
!       msrf(2) = sum variance C
!       msrf(3) = sum variance * length Scale
!       msrf(4) = sum mean * drop diameter
!       msrf(5) = total wetted area
!       msrf(6) = sum sfc mean C
!       msrf(7) = sum mean * front depth (m)
!
!       asrf(1) = sum mean C**2
!       asrf(2) = sum mean * xbar
!       asrf(3) = sum mean * xbar**2
!       asrf(4) = sum mean * ybar
!       asrf(5) = sum mean * ybar**2
!
!       sum total variance = msrf(2) + msrf(6) - msrf(1)**2/16

found_subcell = .FALSE.
IF( dx >= MetGrid(1)%dx )srf_met = .FALSE.

IF( substrate_type == 0 )THEN
  nv = 5
ELSE
  nv = 7
END IF

msrf = 0.0
asrf = 0.0

dxx  = 0.5*dx
dyy  = 0.5*dy
dxxx = 0.5*dxx
dyyy = 0.5*dyy
DO i = 0,3
  iy = INT(i/2)
  xxx = (FLOAT(i-iy-iy)-0.5)*dxx
  yyy = (FLOAT(iy)-0.5)*dyy
  icell = srf%ipgrd(icell0) + i
  IF( ilev == 0 )THEN
    DO iv = 1,nv
      i0 = (iv-1)*srf%mxgrd
      msrf(iv) = msrf(iv) + 4.*srf%ipdat(icell+ivar+i0)
    END DO
    smass   = srf%ipdat(icell+ivar)
    asrf(1) = asrf(1) + 4.*smass*smass
    asrf(2) = asrf(2) + 4.*smass*xxx
    asrf(3) = asrf(3) + 4.*smass*(xxx*xxx + (dxx/RT12)**2)
    asrf(4) = asrf(4) + 4.*smass*yyy
    asrf(5) = asrf(5) + 4.*smass*(yyy*yyy + (dyy/RT12)**2)
  END IF
  IF( srf%ipgrd(icell) > 0 )THEN
    found_subcell = .TRUE.
    DO j = 0,3
      jcell = srf%ipgrd(icell) + j
      DO iv = 1,nv
        i0 = (iv-1)*srf%mxgrd
        msrf(iv) = msrf(iv) + srf%ipdat(jcell+ivar+i0)
      END DO
      jy = INT(j/2)
      xxxx = xxx + (FLOAT(j-jy-jy)-0.5)*dxxx
      yyyy = yyy + (FLOAT(jy)-0.5)*dyyy
      smass   = srf%ipdat(jcell+ivar)
      asrf(1) = asrf(1) + smass*smass
      asrf(2) = asrf(2) + smass*xxxx
      asrf(3) = asrf(3) + smass*(xxxx*xxxx + (dxxx/RT12)**2)
      asrf(4) = asrf(4) + smass*yyyy
      asrf(5) = asrf(5) + smass*(yyyy*yyyy + (dyyy/RT12)**2)
    END DO
  END IF
END DO

IF( ilev == 0 )THEN
  DO iv = 1,nv
    i0 = (iv-1)*srf%mxgrd
    msrf(iv) = msrf(iv) + 16.*srf%ipdat(icell0+ivar+i0)
  END DO
  smass   = srf%ipdat(icell0+ivar)
  asrf(1) = asrf(1) + 16.*smass*smass
  asrf(3) = asrf(3) + 16.*smass*((dx/RT12)**2)
  asrf(5) = asrf(5) + 16.*smass*((dy/RT12)**2)
END IF

IF( asrf(1) > SMALL .AND. msrf(2) > SMALL )THEN !check for underflow

  xxx = asrf(2)/msrf(1)
  yyy = asrf(4)/msrf(1)
  dxx = SQRT(asrf(3)/msrf(1) - xxx*xxx)*RT12
  dyy = SQRT(asrf(5)/msrf(1) - yyy*yyy)*RT12
  xxx = xx + xxx
  yyy = yy + yyy

  IF( .NOT.srf_met )THEN
    CALL get_srf_met( xx,yy )
    IF( nError /= NO_ERROR )GOTO 9999
  END IF

!--- Set cell properties

  dsrf = MAX(msrf(4)/msrf(1),1.E-6)  !Mean droplet diameter (m)

  IF( substrate_type == 0 )THEN   !Impermeable surface
    smass = msrf(1)
    depth = 0.0
  ELSE
    smass = msrf(6)
    depth = msrf(7)/msrf(1)       !Front depth below sfc
  END IF

  CALL srf_evap_rate( msrf(1),dsrf,msrf(5),depth,smass,dts )

  IF( 1.0-rate <= 0.0 )GOTO 9999   ! Skip out if no evaporated mass

!---- Total evaporated vapor mass

  mass = msrf(1)*(1.0-rate)*dx*dy/(xmap_evap*ymap_evap*16.)

  tau_evap = MIN(2.*dts/MAX((1.0-rate),1.E-6),t+dts)

  CALL srf_puff_create( mass,msrf,asrf,xxx,yyy,dxx,dyy,dts,tau_evap,lev2 )
  IF( nError /= NO_ERROR )GOTO 9999

ELSE    !Zero any underflows

  rate   = 0.0
  rate_a = 0.0
  rate_d = 0.0

END IF

CALL srf_cell_evap( srf,ivar,icell0,ilev )

9999 CONTINUE

RETURN
END

!===============================================================================

SUBROUTINE srf_evap_rate( dmass,diam,awet,depth,dsurf,dts )

USE constants_fd
USE default_fd
USE srfevap_fi
USE substrate_fi

IMPLICIT NONE

REAL, INTENT( IN ) :: dmass    !Deposition mass per unit area
REAL, INTENT( IN ) :: diam     !Surface droplet diameter (m)
REAL, INTENT( IN ) :: awet     !Wetted area per unit area
REAL, INTENT( IN ) :: depth    !Front depth (m)
REAL, INTENT( IN ) :: dsurf    !Deposition mass per unit area on surface
REAL, INTENT( IN ) :: dts      !Time step (s)

REAL tx, mtot, msfc, hvd, narea

!--- Set agent properties

CALL lqd( pmatl,tsrf,rhod,cs,hvd,difd )

darea = 0.25*PI*diam**2             !Droplet surface area

IF( porosity == NOT_SET_R )THEN   !Impermeable surface

  CALL srf_drop_evap( diam )
  rate_s = rate
  rate_a = 0.0
  rate_d = 0.0
  dfac   = 0.0

ELSE

  narea = awet/darea                !Droplet number density
  tx    = 0.0

  IF( dsurf > SMALL )THEN            !Surface mass is present
    mtot = dmass/narea               !Single droplet mass total
    msfc = dsurf/narea               !Single droplet surface mass
    CALL srf_drop_evap( diam )
    CALL srf_drop_absorb( mtot,msfc,dts,tx )
    IF( depth > 0.0 )tx = 0.0
  ELSE
    rate   = 0.0
    rate_s = 0.0
    rate_a = 0.0
  END IF

  IF( depth > 0.0 .OR. msfc == 0.0 )THEN !Subsurface mass is exposed
    dfac = rhod*darea*porosity
    CALL srf_drop_desorb( diam,depth,dts-tx )
    rate_d = rate_d*(1.0-tx/dts)
  ELSE
    dfac   = 0.0
    rate_d = 0.0
  END IF

END IF

!---- Set rate factors for surface grid - surface mass (absorption)

IF( porosity /= NOT_SET_R )THEN
  IF( dsurf > SMALL )THEN
    IF( depth > SMALL )THEN
      rate_d = rate_d*(1.0-dsurf/dmass)  ! Adjust rates for mixed evaporation
    END IF
    rate_a = 1.0 - (rate_s+rate_a)*awet*dts/(darea*dsurf)   !Fractional reduction
    rate_a = MAX(0.0,rate_a)
  ELSE
    rate_a = 1.0
  END IF
END IF

!---- Rate factor for total surface mass

rate = 1.0 - (rate_s+rate_d*dfac)*awet*dts/(darea*dmass)      !Fractional reduction

IF( rate < SMALL )rate = 0.0

!---- Rate factor for frontal depth

rate_d = rate_d*dts        !Delta-h frontal depth (m)

RETURN
END

!===============================================================================

SUBROUTINE srf_drop_evap( diam )

USE constants_fd
USE basic_fi
USE srfevap_fi

IMPLICIT NONE

REAL, INTENT( IN ) :: diam   !Surface droplet diameter

!--- Calculate evaporation rate for single surface droplet
!    based on Baines and James model

REAL pe

REAL, EXTERNAL :: sherwd_sec

!------ Compute the Peclet number, Pe

pe = rhoair*(ustr*diam)**2 / (rmuair*difd)

!------ Compute the droplet mass flux

rate = PI2*difd*diam*cs * sherwd_sec( pe )

RETURN
END

!===============================================================================

SUBROUTINE srf_drop_absorb( mtot,msfc,dt,tx )

USE srfevap_fi
USE substrate_fi

IMPLICIT NONE

REAL, INTENT( IN    ) :: mtot        !Total mass (kg)
REAL, INTENT( INOUT ) :: msfc        !Mass on surface (kg)
REAL, INTENT( IN    ) :: dt          !Time step (s)
REAL, INTENT( OUT   ) :: tx          !Evaporation time (s)

REAL kperm, mnew, msub, mtem, tem

!------ Calculate absorption rate

msub = mtot - msfc           !Mass in substrate (kg)

kperm  = k_substrate*pmatl%st*(rhod*darea)**2/pmatl%viscosity
mnew   = SQRT( msub*msub + 2.0*kperm*dt )
rate_a = (mnew - msub)/dt

!----- Limit evaporation/absorption if all surface mass is lost

IF( (rate+rate_a)*dt > msfc )THEN !All sfc mass is absorbed/evaporated

  mtem = mtot*mtot-msub*msub

!---- Calculate time (tx) when sfc mass becomes zero

  IF( rate*(mtot+msub) > 1.0E-3*kperm )THEN  !Check evap time vs absorb time
    tem = mtot + kperm/rate
    tx  = (tem - SQRT(tem*tem-mtem))/rate
  ELSE
    tx  = 0.5*mtem/kperm
  END IF

!---- Adjust sfc mass loss rates to reflect shorter duration

  rate_a = 1.1*(msfc-rate*tx)/dt  !Factor of 1.1 ensures new sfc mass = 0
  rate_s = rate*tx/dt
  msfc   = 0.0                    !Set so desorption occurs for remainder

ELSE

  tx     = 0.0
  rate_s = rate

END IF

RETURN
END

!===============================================================================

SUBROUTINE srf_drop_desorb( diam,depth,dt )

USE constants_fd
USE srfevap_fi
USE substrate_fi

IMPLICIT NONE

REAL, INTENT( IN )  :: diam        !Droplet diameter (m)
REAL, INTENT( IN )  :: depth       !Frontal depth (m)
REAL, INTENT( IN )  :: dt          !Time step (s)

REAL deff, vr, tem, dnew

!------ Calculate desorption rate

deff = difd*porosity/tortuosity

IF( rate == 0.0 )CALL srf_drop_evap( diam )

rate_d = rate

vr = 4.0*rate_d/(PI*diam*diam*cs)

tem  = rate_d*deff*dt/dfac + deff*depth + 0.5*vr*depth*depth
dnew = (-deff + SQRT(deff*deff + 2.0*vr*tem))/vr

!-- Integral of dh/dt = rate_d*deff/dfac/(deff + vr*depth)

rate_d = (dnew - depth)/dt

RETURN
END

!-----------------------------------------------------------------------

REAL FUNCTION sherwd_sec( pe )

IMPLICIT NONE

REAL, INTENT( IN ) :: pe

sherwd_sec = 0.105*pe**0.33333333

RETURN
END

!===============================================================================

SUBROUTINE srf_cell_evap( srf,ivar,icell0,ilev )

USE scipuff_fi
USE srfevap_fi
USE sagstr_fd

IMPLICIT NONE


TYPE( SAGgrid_str ), POINTER      :: srf    !Dep grid structure
INTEGER,             INTENT( IN ) :: ivar   !Start of evap block
INTEGER,             INTENT( IN ) :: ilev   !Grid refinement level
INTEGER,             INTENT( IN ) :: icell0 !Upper cell location

INTEGER i, j, icell, jcell, i0, i1, i2, i3, i4
REAL    rate2, rate_awet, rate_diam

rate2 = rate*rate

rate_awet = rate2**0.3333333
rate_diam = rate_awet*rate_awet

i0 = ivar

i1 = i0 + srf%mxgrd
i2 = i1 + srf%mxgrd
i3 = i2 + srf%mxgrd
i4 = i3 + srf%mxgrd

DO i = 0,3
  icell = srf%ipgrd(icell0) + i
  IF( ilev == 0 )THEN
    srf%ipdat(icell+i0) = srf%ipdat(icell+i0)*rate
    srf%ipdat(icell+i1) = srf%ipdat(icell+i1)*rate2
    srf%ipdat(icell+i2) = srf%ipdat(icell+i2)*rate2
    srf%ipdat(icell+i3) = srf%ipdat(icell+i3)*rate_diam
    srf%ipdat(icell+i4) = srf%ipdat(icell+i4)*rate_awet
  END IF
  IF( srf%ipgrd(icell) > 0 )THEN
    DO j = 0,3
      jcell = srf%ipgrd(icell) + j
      srf%ipdat(jcell+i0) = srf%ipdat(jcell+i0)*rate
      srf%ipdat(jcell+i1) = srf%ipdat(jcell+i1)*rate2
      srf%ipdat(jcell+i2) = srf%ipdat(jcell+i2)*rate2
      srf%ipdat(jcell+i3) = srf%ipdat(jcell+i3)*rate_diam
      srf%ipdat(jcell+i4) = srf%ipdat(jcell+i4)*rate_awet
    END DO
  END IF
END DO

IF( ilev == 0 )THEN
  srf%ipdat(icell0+i0) = srf%ipdat(icell0+i0)*rate
  srf%ipdat(icell0+i1) = srf%ipdat(icell0+i1)*rate2
  srf%ipdat(icell0+i2) = srf%ipdat(icell0+i2)*rate2
  srf%ipdat(icell0+i3) = srf%ipdat(icell0+i3)*rate_diam
  srf%ipdat(icell0+i4) = srf%ipdat(icell0+i4)*rate_awet
END IF

IF( substrate_type /= 0 )THEN

  i1 = i4 + srf%mxgrd
  i2 = i1 + srf%mxgrd
  DO i = 0,3
    icell = srf%ipgrd(icell0) + i
    IF( ilev == 0 )THEN
      srf%ipdat(icell+i1) = srf%ipdat(icell+i1)*rate_a
      srf%ipdat(icell+i2) = srf%ipdat(icell+i2)*rate &
                          + srf%ipdat(icell+i0)*rate_d
    END IF
    IF( srf%ipgrd(icell) > 0 )THEN
      DO j = 0,3
        jcell = srf%ipgrd(icell) + j
        srf%ipdat(jcell+i1) = srf%ipdat(jcell+i1)*rate_a
        srf%ipdat(jcell+i2) = srf%ipdat(jcell+i2)*rate &
                            + srf%ipdat(jcell+i0)*rate_d
      END DO
    END IF
  END DO

  IF( ilev == 0 )THEN
    srf%ipdat(icell0+i1) = srf%ipdat(icell0+i1)*rate_a
    srf%ipdat(icell0+i2) = srf%ipdat(icell0+i2)*rate &
                         + srf%ipdat(icell0+i0)*rate_d
  END IF

END IF

RETURN
END

!===============================================================================

SUBROUTINE SetStaticEvap( p,pq,dt,t_srf,fac_evap )

USE scipuff_fi
USE met_fi
USE step_p_fi
USE srfevap_fi

IMPLICIT NONE

TYPE( puff_str ),    INTENT( IN  ) :: p
TYPE( puff_liquid ), INTENT( IN  ) :: pq
REAL,                INTENT( IN  ) :: dt
REAL,                INTENT( IN  ) :: t_srf
REAL,                INTENT( OUT ) :: fac_evap  !Retained deposition fraction

REAL diam, depth, dmass, awet

TYPE( puff_material ) pmatlIn

!---- Get liquid puff properties

ityp_evap = material(imat)%ioffp + 2
CALL get_puff_material( ityp_evap,pmatlIn )
pmatl = TRANSFER(pmatlIn,pmatl)

dmass = 1.0
depth = 0.0
diam  = pq%d*pmatl%sf
awet  = 3.0*dmass*pmatl%sf*pmatl%sf/(2.0*pq%d*pmatl%rho)

CALL get_srf_met( p%xbar,p%ybar )

CALL srf_evap_rate( dmass,diam,awet,depth,dmass,dt )

IF( rate > 0.0 )THEN
  fac_evap = MIN(1.0,-dt/t_srf/LOG(rate))
ELSE
  fac_evap = 0.0
END IF

RETURN
END

!===============================================================================

SUBROUTINE srf_puff_create( mass,msrf,asrf,xx,yy,dx,dy,dts,tau_evap,lev2 )

USE scipuff_fi
USE met_fi
USE files_fi
USE srfevap_fi

IMPLICIT NONE

REAL,               INTENT( IN    ) :: mass         !Puff vapor mass
REAL, DIMENSION(*), INTENT( IN    ) :: msrf, asrf   !Grid-averaged data
REAL,               INTENT( IN    ) :: xx, yy       !Puff location
REAL,               INTENT( IN    ) :: dx, dy       !Puff spread
REAL,               INTENT( IN    ) :: dts          !Time step
REAL,               INTENT( IN    ) :: tau_evap     !Estimated evaporation time
INTEGER,            INTENT( INOUT ) :: lev2         !Max time level

TYPE ( puff_dynamics ) pd

INTEGER ilev, ios, ifld
REAL    xrele, yrele, delxe, delye, sigxe, sigye, sigze
REAL    hp, hx, hy, uubt, vvbt, uvbt, sly, rat
REAL    sxx, syy, szz, sxy, vol, si, velc, dtx
REAL    rat1, rat2, rat3

REAL,    EXTERNAL :: rlimit
REAL,    EXTERNAL :: SkewFraction
INTEGER, EXTERNAL :: time_level
LOGICAL, EXTERNAL :: check_slope
LOGICAL, EXTERNAL :: IsMulti
INTEGER, EXTERNAL :: allocatePuffAux

sigxe = 0.5*dx/xmap_evap
sigye = 0.5*dy/ymap_evap
sigze = rlimit( VONK*ustr*dts/A,5.*zruf,0.8*zisrf )

CALL get_met( xx,yy,0.0,sigze,0.0,1,Shh=0.5*(sigxe**2+sigye**2),outField=ifld )

dtx   = 0.5*dts
delxe = ub*dtx
delye = vb*dtx

xrele = xx + delxe*xmap_evap
yrele = yy + delye*ymap_evap

si = msrf(3)/msrf(2) + ustr*dts

IF( lsv_oper )CALL reset_lsv( si )

!------ scale turbulence for time-averaging

IF( t_avg /= DEF_VAL_R )THEN
  velc = ub*ub + vb*vb
  CALL turb_rescale( si,sigze,velc )
END IF

CALL set_turb( uubt,vvbt,uvbt )

!---- Set composite horizontal scale for limiting lengthwise puff scale

rat1 = 2.*uubl/(uubt+vvbt+SMALL)
rat2 = 2.*vvbl/(uubt+vvbt+SMALL)
rat3 = (uub+vvb)/(uubt+vvbt+SMALL)
IF( rat1+rat2+rat3 > SMALL )THEN
  sly = rat1*sbls + rat2*sbl + rat3*sby
ELSE ! for bl_type = none
  sly = MAX(sbls,sbl,sby)
END IF

velc = SQRT(ub**2 + vb**2 + qqsh + qqb + qql)

sxx  = sigxe*sigxe + delxe*delxe + uubt*dtx**2
syy  = sigye*sigye + delye*delye + vvbt*dtx**2
sxy  = delxe*delye + uvbt*dtx**2
szz  = sigze*sigze

vol = PI3*SQRT(szz*(sxx*syy-sxy**2))

rat = msrf(2)/asrf(1)

IF( SQRT(rat+1.)*mass < evap_min*vol )GOTO 9999

CALL check_newpuff()
IF( nError /= NO_ERROR )THEN
  eRoutine = 'srf_puff_create'
  GOTO 9999
END IF

npuf = npuf + 1

CALL zero_puff( puff(npuf) )

puff(npuf)%naux = typeID(ityp_evap)%npaux
ios = allocatePuffAux( puff(npuf) )
IF( ios /= 0 )THEN
  nError   = UK_ERROR
  eRoutine = 'srf_puff_create'
  eMessage = 'Error allocating scipuff puff auxiliary array'
  WRITE(eInform,'(A,I0,A,I0)')'Request=',puff(npuf)%naux,'  : Error = ',ios
  GOTO 9999
END IF

CALL setPuffifld( puff(npuf),ifld )

puff(npuf)%c    = mass
puff(npuf)%xbar = xrele
puff(npuf)%ybar = yrele
puff(npuf)%zbar = 0.0
puff(npuf)%sxx  = sxx
puff(npuf)%syy  = syy
puff(npuf)%sxy  = sxy
puff(npuf)%szz  = szz
puff(npuf)%si   = si
puff(npuf)%si2  = rlimit( sly,puff(npuf)%si,0.25*tau_evap*velc )
puff(npuf)%sv   = sigze
puff(npuf)%cfo  = 1.

IF( lsmp )THEN
  CALL set_tscale( puff(npuf)%si,puff(npuf)%si2,velc,puff(npuf)%sr )
END IF

IF( lter )THEN
  CALL get_topogIn( xrele,yrele,hp,hx,hy,ifld )

  puff(npuf)%zbar = hp
  puff(npuf)%zi   = hp + zisrf

  IF( check_slope(hx,hy) )CALL set_puff_rot( puff(npuf),hx,hy )
ELSE
  puff(npuf)%zi = zisrf
END IF

IF( hflxsrf >= 0.0 )THEN
  puff(npuf)%zc = puff(npuf)%zi
ELSE
  puff(npuf)%zc = 0.0
END IF

puff(npuf)%cc   = -mass*MAX(rat,0.01)

puff(npuf)%ityp = ityp_evap
puff(npuf)%inxt = 0
puff(npuf)%iprv = 0
CALL setPuffipgrd( puff(npuf),0 )
puff(npuf)%idtn = 0

puff(npuf)%idtl = time_level( dts )
IF( nError /= NO_ERROR )THEN
  eRoutine = 'srf_puff_create'
  WRITE(lun_log,*,IOSTAT=ios)'******* TIME LEVEL ERROR ********'
  WRITE(lun_log,*,IOSTAT=ios)TRIM(eRoutine)
  WRITE(lun_log,*,IOSTAT=ios)TRIM(eInform)
  WRITE(lun_log,*,IOSTAT=ios)'DT=',dts
  WRITE(lun_log,*,IOSTAT=ios)'DELT=',delt
  WRITE(lun_log,*,IOSTAT=ios)'LEVEL=',puff(npuf)%idtl,MAXTLV
  CALL dump_puff( npuf,puff(npuf) )
  GOTO 9999
END IF

IF( dynamic )CALL get_dynamics( puff(npuf),pd )

CALL siginv( puff(npuf) )
CALL init_tlev( puff(npuf),pd,ilev,xmap_evap,ymap_evap )
puff(npuf)%idtl = MAX(puff(npuf)%idtl,ilev)

lev2 = MAX(lev2,puff(npuf)%idtl)

IF( SkewFraction(0.0) < 1.0 )CALL setPuffiskew( puff(npuf),SKEW_UP )

IF( IsMulti(typeID(ityp_evap)%icls) )THEN
  CALL InitMCvapor( puff(npuf) )
  IF( nError /= NO_ERROR )GOTO 9999
END IF

9999 CONTINUE

RETURN
END

!===============================================================================

SUBROUTINE set_puff_rot( p,hx,hy )

USE puffstruct_fd

IMPLICIT NONE

TYPE( puff_str ),     INTENT( INOUT ) :: p
REAL,                 INTENT( IN    ) :: hx, hy

REAL(8), DIMENSION(3,3) :: arot, atrot, smat, stem
REAL(8), DIMENSION(3)   :: e1, e2, e3

REAL(8) sxx, syy, sxy, s11, s22, s33, hxd, hyd
REAL(8) c, s, ang

!------ Puff horizontal spread

sxx = DBLE(p%sxx); syy = DBLE(p%syy); sxy = DBLE(p%sxy)

!------ Find horizontal principle axes

ang = 0.5*ATAN2( 2.*sxy,sxx-syy )

c = COS(ang); s = SIN(ang)

s11 = c*c*sxx + s*s*syy + 2.*s*c*sxy
s22 = s*s*sxx + c*c*syy - 2.*s*c*sxy
s33 = DBLE(p%szz)

!------ Set tensor with principal axes (assume sxz=syz=0)

smat(:,1) = (/ s11 , 0.D0, 0.D0 /)
smat(:,2) = (/ 0.D0, s22 , 0.D0 /)
smat(:,3) = (/ 0.D0, 0.D0,  s33 /)

!------ Define rotation matrix [ e1
!                                e2
!                                e3 ]

!------ Normal to terrain

hxd = DBLE(hx)
hyd = DBLE(hy)

e3 = (/ -hxd,-hyd,1.D0 /) / SQRT(1.D0+hxd*hxd+hyd*hyd)

!------ Align rotated "x" coordinate along principle axis; set e1.e2 = 0

e1 = (/ c,s,hxd*c+hyd*s /) / SQRT(1.D0+(hxd*c+hyd*s)**2)

!------ Other axis defined by cross-product

e2 = (/ -e1(2)*e3(3)+e3(2)*e1(3),e1(1)*e3(3)-e3(1)*e1(3),-e1(1)*e3(2)+e3(1)*e1(2) /)

!------ Rotation matrix

arot(:,1) = e1
arot(:,2) = e2
arot(:,3) = e3

!------ Rotate principal axes into global coordinate system

atrot = TRANSPOSE( arot )
stem = MATMUL( arot,smat )
smat = MATMUL( stem,atrot )

!------ Copy out into puff structure

p%sxx = SNGL(smat(1,1)); p%sxy = SNGL(smat(2,1)); p%sxz = SNGL(smat(3,1))
                         p%syy = SNGL(smat(2,2)); p%syz = SNGL(smat(3,2))
                                                  p%szz = SNGL(smat(3,3))

RETURN
END
