MODULE MODULE_GOCART_SETTLING
!
! Original GOCART
! KYKang, MKlose, CLWu various changes 
! YShao [2011/09/30] complete update
!
! 11/29/18 - A. Ukhov, bug fix: dust (and sea salt) mass balance was
! violated in the "settling" routine 

CONTAINS

SUBROUTINE gocart_settling_driver(dt,config_flags,t_phy,moist,      &
           chem,rho_phy,dz8w,p8w,p_phy,                             &
           dustin,seasin,dx,g,                                      &
           dustgraset_1,dustgraset_2,dustgraset_3,                  & 
           dustgraset_4,dustgraset_5,                               & 
           setvel_1,setvel_2,setvel_3,setvel_4,setvel_5,imod,       & 
           ids,ide,jds,jde,kds,kde,                                 &
           ims,ime,jms,jme,kms,kme,                                 &
           its,ite,jts,jte,kts,kte)

  USE module_configure
  USE module_state_description
  USE module_data_gocart_dust
  USE module_data_gocart_seas
  USE module_model_constants, ONLY: mwdry
  IMPLICIT NONE

  TYPE(grid_config_rec_type), INTENT(IN) :: config_flags

  INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde,                  &
                         ims,ime, jms,jme, kms,kme,                  &
                         its,ite, jts,jte, kts,kte

  REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_moist), INTENT(IN) :: moist
  REAL, DIMENSION(ims:ime,kms:kme,jms:jme,num_chem), INTENT(INOUT) :: chem
  REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN) :: t_phy,p_phy,dz8w,p8w,rho_phy
  REAL,  DIMENSION( ims:ime , jms:jme, 5 ),                          &
          INTENT(IN   ) ::  dustin,seasin  

  INTEGER, INTENT(IN) :: imod          
  REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT)::                  &
        dustgraset_1,dustgraset_2,dustgraset_3,                      &
        dustgraset_4,dustgraset_5,                                   &
        setvel_1,setvel_2,setvel_3,setvel_4,setvel_5 
                       
  REAL*8,  DIMENSION (1,1,5)          :: graset_dust,grasetvel_dust 
  REAL*8,  DIMENSION (1,1,4)          :: graset_ss,grasetvel_ss

  REAL, INTENT(IN) :: dt,dx,g
  INTEGER          :: kkk,nmx,i,j,k,kk,lmx,iseas,idust
  REAL*8, DIMENSION (1,1,kte-kts+1)   :: tmp,airden,p_mid,delz,rh
  REAL*8, DIMENSION (1,1,kte-kts+1,5) :: ddust
  REAL*8, DIMENSION (1,1,kte-kts+1,4) :: sea_salt
  
  INTEGER :: uoc_flag  ! flag for UoC dust schemes 
!
!   REAL*8, DIMENSION (5), PARAMETER :: den_dust(5)=(/2500.,2650.,2650.,2650.,2650./)
!   REAL*8, DIMENSION (5), PARAMETER :: reff_dust(5)=(/0.73D-6,1.4D-6,2.4D-6,4.5D-6,8.0D-6/)
!   REAL*8, DIMENSION (4), PARAMETER :: den_seas(4)=(/2200.,2200.,2200.,2290./)
!   REAL*8, DIMENSION (4), PARAMETER :: reff_seas(4)=(/0.30D-6,1.00D-6,3.25D-6,7.50D-6/)
!  
  REAL*8 conver, converi
  conver=1.e-9
  converi=1.e9
    
  uoc_flag = 0
  if (config_flags%dust_opt .eq. 4) uoc_flag = 1

  lmx=kte-kts+1

  do j=jts,jte
    do i=its,ite

       IF(config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11) THEN
         do kkk=1,5
            ddust(1,1,kts,kkk)=dustin(i,j,kkk)*conver
         enddo
         kk=0
         do k=kts,kte
            kk=kk+1
            
            p_mid(1,1,kk) =.01*p_phy(i,kte-k+kts,j)
            delz(1,1,kk)  =dz8w(i,kte-k+kts,j)
            airden(1,1,kk)=rho_phy(i,k,j)
            tmp(1,1,kk)= t_phy(i,k,j)
            rh(1,1,kk) = .95
            rh(1,1,kk) = MIN( .95, moist(i,k,j,p_qv) / &
                 (3.80*exp(17.27*(t_phy(i,k,j)-273.)/ &
                 (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))))
            rh(1,1,kk) = max(1.0D-1,rh(1,1,kk))
         enddo
       ELSE  
         kk=0       
         DO k=kts,kte
            kk=kk+1
!
            ddust(1,1,kk,1)=chem(i,k,j,p_dust_1)                         ! chem and dust in [ug/kg]
            ddust(1,1,kk,2)=chem(i,k,j,p_dust_2)
            ddust(1,1,kk,3)=chem(i,k,j,p_dust_3)
            ddust(1,1,kk,4)=chem(i,k,j,p_dust_4)
            ddust(1,1,kk,5)=chem(i,k,j,p_dust_5)
            
            p_mid(1,1,kk)=.01*p_phy(i,kte-k+kts,j)
            delz(1,1,kk)=dz8w(i,kte-k+kts,j)                            ! delz(1) = dz8w(kte), delz(lmx)=dz8w(kts)
            airden(1,1,kk)=rho_phy(i,k,j)
            tmp(1,1,kk)   =t_phy(i,k,j)
            rh(1,1,kk)    = .95
            rh(1,1,kk)    = MIN( .95, moist(i,k,j,p_qv) /               &
                 (3.80*exp(17.27*(t_phy(i,k,j)-273.)/                   &
                 (t_phy(i,k,j)-36.))/(.01*p_phy(i,k,j))))
            rh(1,1,kk)=max(1.0D-1,rh(1,1,kk))
         ENDDO
       ENDIF
       graset_dust(1,1,:)=0.
       graset_ss(1,1,:)=0.                                           

       iseas=0
       idust=1
       CALL settling(1,1,lmx,5,g,dyn_visc,ddust,tmp,p_mid,delz,       &
                     imod,graset_dust,grasetvel_dust, uoc_flag,       & 
                     den_dust,reff_dust,dt,rh,idust,iseas,airden)
          
       IF (config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) THEN
          kk=1
          do kkk=1,5
             if (kkk .le. 4) sea_salt(1,1,kts,kkk)=seasin(i,j,kkk)*conver
             if(ddust(1,1,kk,kkk) .ge. dustin(i,j,kkk)) ddust(1,1,kk,kkk)=dustin(i,j,kkk)
          enddo
          chem(i,kts,j,p_p25i)=chem(i,kts,j,p_p25i) &
                       +.25*(ddust(1,1,kk,1)+.286*ddust(1,1,kk,2))*converi
          chem(i,kts,j,p_p25i)=max(chem(i,kts,j,p_p25i),1.e-16)
          chem(i,kts,j,p_p25j)=chem(i,kts,j,p_p25j) &
                       +.75*(ddust(1,1,kk,1)+.286*ddust(1,1,kk,2))*converi
          chem(i,kts,j,p_p25j)=max(chem(i,kts,j,p_p25j),1.e-16)
          chem(i,kts,j,p_soila)=chem(i,kts,j,p_soila) &
                       +(.714*ddust(1,1,kk,2)+ddust(1,1,kk,3))*converi
          chem(i,kts,j,p_soila)=max(chem(i,kts,j,p_soila),1.e-16)
       ELSE                     
          kk = 0
          DO k = kts,kte
             kk = kk+1
             chem(i,k,j,p_dust_1)=ddust(1,1,kk,1)                         ! dust for size bin 1 [ug/kg]
             chem(i,k,j,p_dust_2)=ddust(1,1,kk,2)                         ! ...
             chem(i,k,j,p_dust_3)=ddust(1,1,kk,3)                         ! ...
             chem(i,k,j,p_dust_4)=ddust(1,1,kk,4)                         ! ...
             chem(i,k,j,p_dust_5)=ddust(1,1,kk,5)                         ! dust for size bin 5 (dust_opt 3: for all size bins) [ug/kg]
             
             sea_salt(1,1,kk,1)=chem(i,k,j,p_seas_1)                    ! salt [ug/kg]
             sea_salt(1,1,kk,2)=chem(i,k,j,p_seas_2)
             sea_salt(1,1,kk,3)=chem(i,k,j,p_seas_3)
             sea_salt(1,1,kk,4)=chem(i,k,j,p_seas_4)              
          ENDDO
       ENDIF
!
! gravitional settling in [ug/m2/s]; from settling, graset_dust in [ug/kg][m/s]
!
       dustgraset_1(i,j)=graset_dust(1,1,1)*airden(1,1,1)*(-1.d0)
       dustgraset_2(i,j)=graset_dust(1,1,2)*airden(1,1,1)*(-1.d0)
       dustgraset_3(i,j)=graset_dust(1,1,3)*airden(1,1,1)*(-1.d0)
       dustgraset_4(i,j)=graset_dust(1,1,4)*airden(1,1,1)*(-1.d0)
       dustgraset_5(i,j)=graset_dust(1,1,5)*airden(1,1,1)*(-1.d0)

       setvel_1(i,j)=grasetvel_dust(1,1,1)                            ! settling velocity [m/s]
       setvel_2(i,j)=grasetvel_dust(1,1,2)
       setvel_3(i,j)=grasetvel_dust(1,1,3)
       setvel_4(i,j)=grasetvel_dust(1,1,4)
       setvel_5(i,j)=grasetvel_dust(1,1,5)

       iseas=1
       idust=0

       CALL settling(1,1,lmx,4,g,dyn_visc,sea_salt,tmp,p_mid,delz,        &
                     imod,graset_ss,grasetvel_ss, uoc_flag,               &
                     den_seas,reff_seas,dt,rh,idust,iseas,airden)
       IF (config_flags%chem_opt == 2 .or. config_flags%chem_opt == 11 ) THEN
          kk=1
          do kkk=1,4
             if(sea_salt(1,1,kk,kkk) .ge. seasin(i,j,kkk))sea_salt(1,1,kk,kkk)=seasin(i,j,kkk)
          enddo
          chem(i,kts,j,p_naai)=chem(i,kts,j,p_naai) &
                       +.25*(sea_salt(1,1,kk,1)+.942*sea_salt(1,1,kk,2))*converi
          chem(i,kts,j,p_naai)=max(1.e-16,chem(i,kts,j,p_naai))
          chem(i,kts,j,p_naaj)=chem(i,kts,j,p_naaj) &
                       +.75*(sea_salt(1,1,kk,1)+.942*sea_salt(1,1,kk,2))*converi
          chem(i,kts,j,p_naaj)=max(1.e-16,chem(i,kts,j,p_naaj))
          chem(i,kts,j,p_seas)=chem(i,kts,j,p_seas) &
                       +(.058*sea_salt(1,1,kk,2)+sea_salt(1,1,kk,3))*converi
          chem(i,kts,j,p_seas)=max(1.e-16,chem(i,kts,j,p_seas))
       ELSE
          kk=0
          DO k=kts,kte
             kk=kk+1
             chem(i,k,j,p_seas_1)=sea_salt(1,1,kk,1)
             chem(i,k,j,p_seas_2)=sea_salt(1,1,kk,2)
             chem(i,k,j,p_seas_3)=sea_salt(1,1,kk,3)
             chem(i,k,j,p_seas_4)=sea_salt(1,1,kk,4)
          ENDDO
       ENDIF 
       
    enddo  ! i
  enddo  ! j 

 END SUBROUTINE gocart_settling_driver


 subroutine settling(imx,jmx,lmx,nmx,g0,dyn_visc,tc,tmp,p_mid,delz, &
                    imod,graset, grasetvel, uoc,                    & 
                    den_in,reff_in,dt,rh,idust,iseas,airden)

! ****************************************************************************
! *                                                                          *
! *  Calculate the loss by settling, using an implicit method                *
! *                                                                          *
! *  Input variables:                                                        *
! *    SIGE(k)         - sigma coordinate of the vertical edges              *
! *    PS(i,j)         - Surface pressure (mb)                               *
! *    TMP(i,j,k)      - Air temperature  (K)                                *
! *    CT(i,j)         - Surface exchange coeff for moisture
! *                                                                          *
! **************************************************************************** 

  IMPLICIT  NONE

  INTEGER, INTENT(IN) :: imx, jmx, lmx, nmx,iseas,idust
  INTEGER             :: ntdt
  REAL,    INTENT(IN) :: dt,g0,dyn_visc
  REAL*8,  INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx),  &
                         rh(imx,jmx,lmx), p_mid(imx,jmx,lmx),airden(imx,jmx,lmx)
! 
  REAL*8, INTENT(IN)  :: den_in(nmx), reff_in(nmx)

  REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)

  INTEGER, INTENT(IN)   :: imod, uoc
  REAL*8, INTENT(INOUT) :: graset(imx,jmx,nmx)
  REAL*8, INTENT(OUT)   :: grasetvel(imx,jmx,nmx)
  
  REAL*8    :: den(nmx), reff(nmx)                    ! local variables here
  REAL*8    :: dt_settl(nmx), rcm(nmx), rho(nmx)
  INTEGER   :: ndt_settl(nmx)
  REAL*8    :: dzmin, vsettl, dtmax, pres, rhb, rwet(nmx), ratio_r(nmx)
  REAL*8    :: c_stokes, free_path, c_cun, viscosity, growth_fac
  REAL*8    :: vd_cor(lmx), vd_wk1
  INTEGER   :: k, n, i, j, l, l2

  REAL*8    :: transfer_to_below_level,temp_tc

! for sea-salt:
  REAL*8, PARAMETER :: c1=0.7674, c2=3.079, c3=2.573E-11, c4=-1.424 

! for OMP:
  REAL*8 :: rwet_priv(nmx), rho_priv(nmx)

! IF (type) /= 'dust' .AND. TRIM(aero_type) /= 'sea_salt') RETURN
  IF ( idust.ne.1 .and. iseas.ne.1 ) RETURN
  WHERE ( tc(:,:,:,:) < 0.0 ) tc(:,:,:,:) = 1.0D-32

  den = den_in  
  reff = reff_in  
  dzmin = MINVAL(delz(:,:,:))
  IF (idust == 1) growth_fac = 1.0
  IF (iseas == 1) growth_fac = 3.0
!
  IF (idust == 1 .and. uoc == 1) then
      den(1) = 2650.            ! constant density is use in UoC dust emission schemes; 
!                               dust radii are now consistent with GOCART dust schemes.  [mklose, 03082015]
  ENDIF

  DO k = 1, nmx                                                   ! k for different size bins
!
! Settling velocity (m/s) for each tracer (Stokes Law)
! DEN         density                        (kg/m3)
! REFF        effective radius               (m)
! dyn_visc    dynamic viscosity              (kg/m/s)
! g0          gravity                        (m/s2)
! 3.0         corresponds to a growth of a factor 3 of radius with 100% RH
! 0.5         upper limit with temp correction
!
     vsettl = 4.0/9.0 * g0 * den(k) * (growth_fac*reff(k))**2 / dyn_visc
!
! Determine the maximum time-step satisying CFL: dt <= (dz)_min / v_settl
!
     ntdt = INT(dt)
     dtmax = dzmin / vsettl
     ndt_settl(k) = MAX( 1,INT(ntdt/dtmax) )

! Limit maximum number of iterations
     IF (ndt_settl(k) > 12) ndt_settl(k) = 12
     dt_settl(k) = REAL(ntdt) / REAL(ndt_settl(k))

! Particles radius in centimeters
     IF (iseas.eq.1)rcm(k) = reff(k)*100.0
     IF (idust.eq.1) THEN
      rwet(k) = reff(k)
      ratio_r(k) = 1.0
      rho(k) = den(k)
     ENDIF
  ENDDO

! Solve the bidiagonal matrix (l,l)

!$OMP PARALLEL DO &
!$OMP DEFAULT( SHARED ) &
!$OMP PRIVATE( i,   j,   l,   l2, n,   k,   rhb, rwet_priv, ratio_r, c_stokes)&
!$OMP PRIVATE( free_path, c_cun, viscosity, rho_priv, vd_cor )

! Loop over latitudes

  DO j = 1,jmx                      ! lat loop
    DO i = 1,imx                    ! lon loop
      DO k = 1,nmx                  ! bin loop
        graset(i,j,k)=0.
        grasetvel(i,j,k)=0.
    
        IF (idust.eq.1) THEN
           rwet_priv(k) = rwet(k)
           rho_priv(k)  = rho(k)
        END IF

        DO n = 1,ndt_settl(k)        ! time loop 

          transfer_to_below_level=0

          DO l = lmx,1,-1            ! height loop, from top
             l2 = lmx - l + 1

             IF (iseas.eq.1) THEN
                rhb = MIN(9.9D-1, rh(i,j,l))                             ! Aerosol growth with relative humidity (Gerber, 1985)
                rwet_priv(k) = 0.01*(c1*rcm(k)**c2/(c3*rcm(k)**c4 -  &
                LOG10(rhb)) + rcm(k)**3)**0.33                           ! td changed to LOG10
                ratio_r(k) = (reff(k)/rwet_priv(k))**3.0
             END IF

             c_stokes = 1.458E-6*tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4)    ! Dynamic viscosity
             free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l))           ! Free path as func of pres(mb) and temp(K); order of p_mid: top->sfc

             c_cun = 1.0+free_path/rwet_priv(k)*                     &
             (1.257 + 0.4*EXP(-1.1*rwet_priv(k)/free_path))              ! Slip correction
             viscosity = c_stokes / c_cun                                ! Corrected dynamic viscosity (kg/m/s)

             IF (iseas.eq.1) THEN
                rho_priv(k) = ratio_r(k)*den(k) + (1.-ratio_r(k))*1000.
             END IF

             vd_cor(l) = 2./9.*g0*rho_priv(k)*rwet_priv(k)**2/viscosity  ! Settling velocity, depends on temp

             ! Update mixing ratio; order of delz: top->sfc
             temp_tc=tc(i,j,l,k)                                          ! temp_tc - for temporal storage [ug/kg]            
             vd_wk1 = dt_settl(k)*vd_cor(l)/delz(i,j,l2)                  ! fraction to leave level

             tc(i,j,l,k)   =  tc(i,j,l,k)*(1.- vd_wk1)+transfer_to_below_level          ! [ug/kg]
            
             IF (l==1) THEN
                graset(i,j,k)=graset(i,j,k)+vd_cor(l)*temp_tc/ndt_settl(k) ! [ug/kg][m/s]
                grasetvel(i,j,k)=vd_cor(l)                                 ! [m/s]
             ELSE
                transfer_to_below_level = (temp_tc*vd_wk1)*((delz(i,j,l2)*airden(i,j,l))/(delz(i,j,l2+1)*airden(i,j,l-1))) ! [ug/kg]
             ENDIF
          ENDDO                 !l, height
        ENDDO                   !n, time
      ENDDO                     !k, bin
!
    ENDDO                     !i
  END DO                      !j
!$OMP END PARALLEL DO

END SUBROUTINE settling

END MODULE MODULE_GOCART_SETTLING
