!*******************************************************************************
!$RCSfile$
!$Revision$
!$Date$
!*******************************************************************************
SUBROUTINE splitz( p,lev1,lev2,nchg )

USE scipuff_fi
USE step_p_fi

!----- ***** this routine may be hazardous to your puff's health *****

IMPLICIT NONE

TYPE( puff_str ), INTENT( INOUT ) :: p
INTEGER,          INTENT( IN    ) :: lev1
INTEGER,          INTENT( IN    ) :: lev2
INTEGER,          INTENT( INOUT ) :: nchg

INTEGER iprv, inxt, idtn, il, ios

INTEGER, EXTERNAL :: deallocatePuffAux

!------ Save variables to move into original puff

iprv = p%iprv
inxt = p%inxt
idtn = p%idtn
il   = p%idtl


!------ Vertical split depends on whether puff is well-mixed

IF( sz >= WMX*(p%zc-hp) .AND. p%zc > 0.0 )THEN
  CALL splitz_wm( p,lev1,lev2,nchg )
ELSE
  CALL splitz_nwm( p,lev1,lev2,nchg )
END IF
IF( nError /= NO_ERROR )GOTO 9999

!------ Move last split puff into original puff location

ios = deallocatePuffAux(p)
IF( ios /= 0 )THEN
  nError   = UK_ERROR
  eRoutine = 'splitz'
  eMessage = 'Error deallocating puff aux array'
  WRITE(eInform,'(A,I0)')'IOS=',ios
  GOTO 9999
END IF
CALL move_puff( puff(npuf),p )

p%iprv = iprv
p%inxt = inxt
p%idtn = idtn

npuf = npuf - 1

!---- Check for change in time level

IF( il /= p%idtl )nchg = nchg + 1

9999 CONTINUE

RETURN
END

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

SUBROUTINE splitz_nwm( p,lev1,lev2,nchg )

USE scipuff_fi
USE step_p_fi
USE mpi_fi, ONLY: pufNo
USE StepMC_mpi_fi, ONLY: puffNos

!------ split part of puff above zi if it previously was below zi and zi has
!       decreased - assumes puff is not well-mixed

IMPLICIT NONE

REAL, PARAMETER :: ZFAC = 3.0000000 ! set to match value in step_p
REAL, PARAMETER :: R2   = 0.7071068
REAL, PARAMETER :: EPS  = 1.E-2     ! set to match value in gauss_int

TYPE( puff_str ), INTENT( INOUT ) :: p
INTEGER,          INTENT( IN    ) :: lev1
INTEGER,          INTENT( IN    ) :: lev2
INTEGER,          INTENT( INOUT ) :: nchg

INTEGER npufo, jpuf, i, nzp
REAL    szz, dzz, tfac, zbot, zrfl
REAL    zlow, ctot, frac, dzp, delz
LOGICAL lrfl_top, lrfl_bot

REAL,    EXTERNAL :: erfc

npufo = npuf

!------ determine vertical range for newly created puffs

zbar = p%zbar
dzz  = ZFAC*sz
tfac = R2/sz
zrfl = p%zc

!------ top-most extent

IF( zrfl >= zbar + dzz .OR. zrfl == 0. )THEN
  ztop = zbar + dzz
  lrfl_top = .FALSE.
ELSE
  ztop = zrfl
  lrfl_top = .TRUE.
END IF

!------ bottom-most extent

zbot = zbar-dzz
zlow = MAX(p%zi,zbot)

lrfl_bot = (zbar - hp) < dzz

!------ check if splitting will create puffs of significant (relative) mass

ctot = 0.5 * erfc( (zlow-zbar)*tfac )

IF( ctot < EPS )THEN

  CALL check_newpuff()
  IF( nError /= NO_ERROR) GOTO 9999

  npuf = npuf + 1

  CALL copy_puff( p,puff(npuf) )
  IF( nError /= NO_ERROR) GOTO 9999

  RETURN

END IF

!------ create puffs

IF( lsplitz )THEN
  nzp  = MAX(INT((ztop-zlow)/dzg+0.5),1)
ELSE
  nzp = 1
END IF

dzp  = (ztop-zlow)/FLOAT(nzp)
ctot = 0.

DO i = 1,nzp

  ztop = zlow + dzp

  CALL check_newpuff()
  IF( nError /= NO_ERROR) GOTO 9999

  npuf = npuf + 1

  CALL chop_puff( p,sz,tfac,zlow,ztop,zrfl,hp,lrfl_top,lrfl_bot, &
                  frac,delz,szz )

  IF( frac == 0.0 )THEN
    npuf = npuf - 1
  ELSE
    CALL set_puff( p,szz,frac,delz,xmap,ymap,lev1,lev2,nchg )
    puffNos(pufno)%ipuf = pufno
    puffNos(pufno)%npuf = npuf
    puffNos(pufno)%frac = frac
  END IF

  ctot = ctot + frac
  zlow = ztop

END DO

!------ put remainder below current zi if necessary

IF( zbot < p%zi )THEN

  CALL check_newpuff()
  IF( nError /= NO_ERROR )GOTO 9999

  npuf = npuf + 1

  zlow = hp
  ztop = p%zi

  CALL chop_puff( p,sz,tfac,zlow,ztop,zrfl,hp,lrfl_top,lrfl_bot, &
                  frac,delz,szz )

  IF( frac == 0.0 )THEN
    npuf = npuf - 1
  ELSE
    CALL set_puff( p,szz,frac,delz,xmap,ymap,lev1,lev2,nchg )
    szz  = 0.65*(ztop-zlow)**2 ! well-mixed
    IF( puff(npuf)%szz > 0.8*szz )THEN
       puff(npuf)%szz  = szz
       puff(npuf)%idtl = p%idtl
    END IF
    puff(npuf)%zc = zlim
    puffNos(pufno)%ipuf = pufno
    puffNos(pufno)%npuf = npuf
    puffNos(pufno)%frac = frac
  END IF

  ctot = ctot + frac

END IF

!------ check mass

IF( ctot > 0.5 )THEN

  frac = 1.0/ctot
  DO jpuf = npufo+1,npuf
    CALL scale_psum( puff(jpuf),frac )
    CALL scale_paux( puff(jpuf),frac )
  END DO

ELSE

  nError   = UK_ERROR
  eMessage = 'Error partitioning puff mass about inversion'
  eInform  = 'Try reducing vertical resolution'
  GOTO 9999

END IF

RETURN

9999 CONTINUE

eRoutine ='splitz_nmw'

RETURN
END

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

SUBROUTINE splitz_wm( p,lev1,lev2,nchg )

USE scipuff_fi
USE step_p_fi
USE mpi_fi, ONLY: pufNo
USE StepMC_mpi_fi, ONLY: puffNos

!------ split part of puff above zi if it previously was below zi and zi has
!       decreased - assumes puff IS well-mixed

IMPLICIT NONE

TYPE( puff_str ), INTENT( INOUT ) :: p
INTEGER,          INTENT( IN    ) :: lev1
INTEGER,          INTENT( IN    ) :: lev2
INTEGER,          INTENT( INOUT ) :: nchg

INTEGER jpuf, npufo, i, nzp, ios
REAL    dzz, zbot, zlow, dzp, ctot, frac

INTEGER, EXTERNAL :: deallocatePuffAux

npufo = npuf

!------ determine vertical range for newly created puffs

ztop = p%zc
zbot = p%zi

dzz  = ztop - hp

!------ create puffs

IF( lsplitz )THEN
  nzp = MAX(INT((ztop-zbot)/dzg+0.5),1)
ELSE
  nzp = 1
END IF
dzp  = (ztop-zbot)/FLOAT(nzp)
zlow = zbot

ctot = 0.

DO i = 1,nzp

  ztop = zlow + dzp

  CALL check_newpuff()
  IF( nError /= NO_ERROR )GOTO 9999

  npuf = npuf + 1

  CALL chop_puffwm( p,dzz,zlow,ztop,xmap,ymap,lev1,lev2,nchg,frac )

  IF( frac == 0.0 )THEN
    ios = deallocatePuffAux( puff(npuf) )
    IF( ios /= 0 )THEN
      nError = SZ_ERROR
      eRoutine = 'splitz_wm'
      eMessage = 'Error deallocating puff aux array'
      WRITE(eInform,'(A,I0)')'IOS=',ios
      GOTO 9999
    END IF
    npuf  = npuf - 1
    puffNos(pufno)%ipuf = pufno
    puffNos(pufno)%npuf = npuf
    puffNos(pufno)%frac = frac
  END IF

  ctot = ctot + frac
  zlow = ztop

END DO

!------ put remainder below current zi

CALL check_newpuff()
IF( nError /= NO_ERROR )GOTO 9999

npuf = npuf + 1

zlow = hp
ztop = p%zi

CALL chop_puffwm( p,dzz,zlow,ztop,xmap,ymap,lev1,lev2,nchg,frac )

IF( frac == 0.0 )THEN
  ios = deallocatePuffAux( puff(npuf) )
  IF( ios /= 0 )THEN
    nError = SZ_ERROR
    eRoutine = 'splitz_wm'
    eMessage = 'Error deallcoating puff aux array'
    WRITE(eInform,'(A,I0)')'IOS=',ios
    GOTO 9999
  END IF
  npuf = npuf - 1

ELSE

  ctot = ctot + frac

  puff(npuf)%szz  = 0.65*(ztop-zlow)**2 ! set well-mixed
  puff(npuf)%idtl = p%idtl
  puff(npuf)%zc   = zlim

  puffNos(pufno)%ipuf = pufno
  puffNos(pufno)%npuf = npuf
  puffNos(pufno)%frac = frac
END IF

!------ check mass

IF( ctot > 0.5 )THEN

  frac = 1.0/ctot
  DO jpuf = npufo+1,npuf
    CALL scale_psum( puff(jpuf),frac )
    CALL scale_paux( puff(jpuf),frac )
  END DO

ELSE

  nError   = UK_ERROR
  eMessage = 'Error partitioning puff mass about inversion'
  eInform  = 'Try reducing vertical resolution'
  GOTO 9999

END IF

RETURN

9999 CONTINUE

eRoutine ='splitz_wm'

RETURN
END

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

SUBROUTINE chop_puff( p,sz,tfac,zlow,ztop,zrfl,hp, &
                      lrfl_top,lrfl_bot,frac,delz,szz )

USE scipuff_fi

IMPLICIT NONE

TYPE( puff_str ), INTENT( INOUT ) :: p
REAL,             INTENT( IN    ) :: sz
REAL,             INTENT( IN    ) :: tfac
REAL,             INTENT( IN    ) :: zlow
REAL,             INTENT( IN    ) :: ztop
REAL,             INTENT( IN    ) :: zrfl
REAL,             INTENT( IN    ) :: hp
LOGICAL,          INTENT( IN    ) :: lrfl_top, lrfl_bot
REAL,             INTENT( OUT   ) :: frac
REAL,             INTENT( OUT   ) :: delz
REAL,             INTENT( OUT   ) :: szz

REAL zbar_top, zbar_bot, zbar

!------ compute section mass fraction and vertical moments

CALL gauss_int( p%zbar,zlow,ztop,sz,p%szz,tfac,frac,delz,szz )

IF( frac == 0. )RETURN

!------ reflections if necessary

IF( lrfl_top )THEN

  zbar_top = 2.*zrfl - p%zbar
  CALL refl_int( zbar_top,p%zbar,zlow,ztop,sz,szz,delz,frac,tfac )

END IF

IF( lrfl_bot )THEN

  zbar_bot = 2.*hp - p%zbar
  CALL refl_int( zbar_bot,p%zbar,zlow,ztop,sz,szz,delz,frac,tfac )

  IF( lrfl_top )THEN

    zbar = 2.*zrfl - zbar_bot
    CALL refl_int( zbar,p%zbar,zlow,ztop,sz,szz,delz,frac,tfac )

    zbar = 2.*hp - zbar_top
    CALL refl_int( zbar,p%zbar,zlow,ztop,sz,szz,delz,frac,tfac )

  END IF

END IF

RETURN
END

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

SUBROUTINE chop_puffwm( p,dzz,zlow,ztop,xmap,ymap,lev1,lev2,nchg,frac )

USE scipuff_fi

IMPLICIT NONE

REAL, PARAMETER :: SZFAC = 0.16

TYPE( puff_str ), INTENT( INOUT ) :: p
REAL,             INTENT( IN    ) :: dzz
REAL,             INTENT( IN    ) :: zlow
REAL,             INTENT( IN    ) :: ztop
REAL,             INTENT( IN    ) :: xmap, ymap
INTEGER,          INTENT( IN    ) :: lev1, lev2
INTEGER,          INTENT( OUT   ) :: nchg
REAL,             INTENT( OUT   ) :: frac

REAL dzzo, szz, delz

dzzo = ztop - zlow
frac = dzzo/dzz
szz  = SZFAC*dzzo*dzzo
delz = zlow + 0.5*dzzo - p%zbar

CALL set_puff( p,szz,frac,delz,xmap,ymap,lev1,lev2,nchg )

RETURN
END

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

SUBROUTINE gauss_int( zbar,zlow,ztop,sz,szz,tfac,q0,delz,szzo )

USE constants_fd

!------ compute section mass fraction and vertical moments

IMPLICIT NONE

REAL, PARAMETER :: EPS  = 1.E-2

REAL, INTENT( IN  ) :: zbar
REAL, INTENT( IN  ) :: zlow
REAL, INTENT( IN  ) :: ztop
REAL, INTENT( IN  ) :: sz
REAL, INTENT( IN  ) :: szz
REAL, INTENT( IN  ) :: tfac
REAL, INTENT( OUT ) :: q0
REAL, INTENT( OUT ) :: delz
REAL, INTENT( OUT ) :: szzo

REAL t1, t2, e1, e2, q1, q2

REAL, EXTERNAL :: erfc

t1   = (zlow-zbar)*tfac
t2   = (ztop-zbar)*tfac
e1   = EXP(-t1**2)
e2   = EXP(-t2**2)
q0   = 0.5*(erfc( t1 ) - erfc( t2 ))

IF( q0 < EPS )THEN
  q0   = 0.
  delz = 0.
  szzo = 1.E+10
  GOTO 9999
END IF

q1   = (e1 - e2) / (SQRT2*SQRTPI)
q2   = (t1*e1 - t2*e2) / SQRTPI

delz = sz*q1/q0
delz = MAX(MIN(delz,ztop-zbar),zlow-zbar)
szzo = szz*(1.+q2/q0-(q1/q0)**2)
IF( szzo <= 0. )szzo = 0.0833333*(ztop-zlow)**2

9999 CONTINUE

RETURN
END

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

SUBROUTINE set_puff( p,szz,frac,delz,xmap,ymap,lev1,lev2,nchg )

USE scipuff_fi

!------ set new puff

IMPLICIT NONE

TYPE( puff_str ), INTENT( INOUT ) :: p
REAL,             INTENT( IN    ) :: szz
REAL,             INTENT( IN    ) :: frac
REAL,             INTENT( IN    ) :: delz
REAL,             INTENT( IN    ) :: xmap, ymap
INTEGER,          INTENT( IN    ) :: lev1, lev2
INTEGER,          INTENT( INOUT ) :: nchg

INTEGER il
REAL    delx, dely, temz, sxz, syz
REAL    dtk, del

puff(npuf)%sxx = p%sxx
puff(npuf)%sxy = p%sxy
puff(npuf)%syy = p%syy

CALL set_puff_xc( puff(npuf),p,frac )

delx = p%sxz/p%szz * delz
dely = p%syz/p%szz * delz

puff(npuf)%xbar = p%xbar + delx*xmap
puff(npuf)%ybar = p%ybar + dely*ymap
puff(npuf)%zbar = p%zbar + delz

IF( global_lon )CALL SetGlobalLon( puff(npuf)%xbar )

temz = SQRT(szz/p%szz)
sxz  = p%sxz*temz
syz  = p%syz*temz

puff(npuf)%sxz  = sxz
puff(npuf)%syz  = syz
puff(npuf)%szz  = szz

!------ check time level

dtk = 0.5*delz2/MAX(p%zwc/p%c,1.E-6)

del = delt
il  = 0

DO WHILE( dtk < del )
  il  = il + 1
  del = 0.5*del
END DO

il = MAX(il,lev1)
il = MAX(il,p%idtl)

!------ reduce diffusivity if time step is too big

IF( il > lev2 )THEN
  puff(npuf)%zwc = puff(npuf)%zwc*2.**(lev2-il)
  il = lev2
END IF

IF( il /= p%idtl )nchg = nchg + 1

puff(npuf)%ityp = p%ityp
puff(npuf)%inxt = 0
puff(npuf)%iprv = 0
puff(npuf)%idtl = il
puff(npuf)%idtn = 0
puff(npuf)%ipgd = p%ipgd

RETURN
END

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

SUBROUTINE set_puff_xc( p1,p2,frac )

USE scipuff_fi

IMPLICIT NONE

TYPE( puff_str), INTENT( INOUT ) :: p1
TYPE( puff_str), INTENT( IN    ) :: p2
REAL,            INTENT( IN    ) :: frac

TYPE( puff_str_xc) p1_xc, p2_xc

p1_xc = TRANSFER(p1,p1_xc)
p2_xc = TRANSFER(p2,p2_xc)

p1_xc%psum = frac*p2_xc%psum

p1_xc%prat = p2_xc%prat

p1 = TRANSFER(p1_xc,p1)

CALL copy_scale_paux( p1,p2,frac )

RETURN
END

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

SUBROUTINE refl_int( zbar,zbari,zlow,ztop,sz,szz,delz,frac,tfac )

!------ integrate reflected puff at z=zbar and combine with existing puff

IMPLICIT NONE

REAL, INTENT( IN    ) :: zbar
REAL, INTENT( IN    ) :: zbari
REAL, INTENT( IN    ) :: zlow, ztop
REAL, INTENT( IN    ) :: sz
REAL, INTENT( INOUT ) :: szz
REAL, INTENT( INOUT ) :: delz
REAL, INTENT( INOUT ) :: frac
REAL, INTENT( IN    ) :: tfac

REAL szzi, fracr, delr, szzr, ftot, temz

szzi = sz*sz

CALL gauss_int( zbar,zlow,ztop,sz,szzi,tfac,fracr,delr,szzr )

IF( fracr > 0. )THEN

  ftot = frac + fracr
  delr = delr + zbar - zbari
  temz = (frac*delz + fracr*delr)/ftot
  szz  = (frac*(delz*delz+szz) + fracr*(delr*delr+szzr))/ftot - temz**2
  delz = temz

  frac = ftot

END IF

RETURN
END

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

SUBROUTINE split_zi( p,zinv )

USE scipuff_fi
USE step_p_fi

IMPLICIT NONE

REAL, PARAMETER :: ZFAC = 3.0000000 ! set to match value in step_p
REAL, PARAMETER :: R2   = 0.7071068

TYPE( puff_str ), INTENT( INOUT ) :: p
REAL,             INTENT( IN    ) :: zinv

REAL    dzz, tfac, zrfl
REAL    temz, szz
LOGICAL lrfl_top, lrfl_bot

!------ check space for new puff

CALL check_newpuff()
IF( nError /= NO_ERROR )GOTO 9999

!------ determine vertical range for newly created puffs

zbar = p%zbar
dzz  = ZFAC*sz
tfac = R2/sz
zrfl = p%zc

!------ top-most extent

IF( zrfl >= zbar + dzz .OR. zrfl == 0. )THEN
  ztop = zbar + dzz
  lrfl_top = .FALSE.
ELSE
  ztop = zrfl
  lrfl_top = .TRUE.
END IF

!------ bottom reflection condition

lrfl_bot = (zbar - hp) < dzz

!------ create new puff above zi with remainder below (in original puff)

CALL chop_puff_zi( p,sz,tfac,zinv,ztop,zrfl,hp,lrfl_top,lrfl_bot )

!----- limit capped puff to well-mixed value

temz = WMX*(zinv-hp)
szz  = MIN(p%szz,temz*temz)

temz  = SQRT(szz/p%szz)
p%sxz = p%sxz*temz
p%syz = p%syz*temz
p%szz = szz
sz    = sz*temz

9999 CONTINUE

RETURN
END

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

SUBROUTINE chop_puff_zi( p,sz,tfac,zinv,ztop,zrfl,hp,lrfl_top,lrfl_bot )

USE scipuff_fi
USE mpi_fi, ONLY: pufNo, isSerial
USE StepMC_mpi_fi, ONLY: puffNos

USE files_fi

IMPLICIT NONE

TYPE( puff_str ), INTENT( INOUT ) :: p
REAL,             INTENT( IN    ) :: sz
REAL,             INTENT( IN    ) :: tfac
REAL,             INTENT( IN    ) :: zinv
REAL,             INTENT( IN    ) :: ztop
REAL,             INTENT( IN    ) :: zrfl
REAL,             INTENT( IN    ) :: hp
LOGICAL,          INTENT( IN    ) :: lrfl_top, lrfl_bot

REAL    zlow, frac
REAL    szz, delz, temz, sxz, syz
REAL    aspect, fracFrac, minFrac

REAL, PARAMETER :: MIN_FRAC = 0.05
REAL, PARAMETER :: MAX_FRAC = 0.5 - MIN_FRAC
REAL, PARAMETER :: MIN_ASPECT = 1000.0             !Aspect**2 => Aspect=31.6
REAL, PARAMETER :: MAX_ASPECT = 100.0*MIN_ASPECT   !Aspect**2 => Aspect=316.2

!----- set puff above zi

aspect = MAX(p%sxx,p%syy)/p%szz                    !Aspect**2
IF( aspect > 10.0*MAX_ASPECT )GOTO 9999            !Aspect**2 => Aspect = 1000.0

zlow = zinv
CALL chop_puff( p,sz,tfac,zlow,ztop,zrfl,hp,lrfl_top,lrfl_bot,frac,delz,szz )

fracFrac = 0.5*LOG10(MAX(MIN(aspect,MAX_ASPECT),MIN_ASPECT)/MIN_ASPECT)
minFrac = min_FRAC + fracFrac*MAX_FRAC

IF( frac <= minFrac )GOTO 9999

npuf = npuf + 1

CALL zero_puff( puff(npuf) )

IF( .NOT.isSerial )THEN
  WHERE( puffNos(:)%ipuf == pufNo)
    puffNos(:)%npuf = npuf
    puffNos(:)%frac = frac
  END WHERE
END IF
puff(npuf)%sxx = p%sxx
puff(npuf)%sxy = p%sxy
puff(npuf)%syy = p%syy

CALL set_puff_xc( puff(npuf),p,frac )

puff(npuf)%zc = p%zc
puff(npuf)%zi = p%zi

puff(npuf)%xbar = p%xbar
puff(npuf)%ybar = p%ybar
puff(npuf)%zbar = p%zbar + delz

temz = SQRT(szz/p%szz)
sxz  = p%sxz*temz
syz  = p%syz*temz

puff(npuf)%sxz  = sxz
puff(npuf)%syz  = syz
puff(npuf)%szz  = szz

puff(npuf)%ityp = p%ityp
puff(npuf)%inxt = 0
puff(npuf)%iprv = 0
puff(npuf)%idtl = p%idtl
puff(npuf)%idtn = p%idtn
puff(npuf)%ipgd = p%ipgd

!----- set puff below zi (don't move centroid or change second-moments)

frac = 1. - frac
CALL scale_puff( p,frac )

p%zc = zinv

9999 CONTINUE

RETURN
END
