MODULE MODULE_VASH_SETTLING

CONTAINS

SUBROUTINE vash_settling_driver(dt,config_flags,t_phy,moist,               &
         chem,rho_phy,dz8w,p8w,p_phy,                                      &
         ash_fall,dx,g,                                                    &
         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 ),                        &
          INTENT(INOUT   ) ::  ash_fall

  REAL, INTENT(IN   ) :: dt,dx,g
  integer :: nmx,i,j,k,kk,lmx,iseas,idust
  real*8, DIMENSION (1,1,kte-kts+1) :: tmp,airden,airmas,p_mid,delz,rh
  real*8, DIMENSION (1,1,kte-kts+1,4) :: sea_salt
!srf
  real*8, DIMENSION (1,1,kte-kts+1,10) :: ash
  real*8, DIMENSION (10), PARAMETER :: den_ash(10)=(/2500.,2500.,2500.,2500.,2500., &
                                                     2500.,2500.,2500.,2500.,2500. /)
  real*8, DIMENSION (10), PARAMETER :: reff_ash(10)=(/0.5000D-3,&! 1.00 mm diameter 
                                                      0.3750D-3,&! 0.75 mm
						      0.1875D-3,&!
						      93.750D-6,&!
						      46.875D-6,&!
						      23.437D-6,&!
						      11.719D-6,&!
						      05.859D-6,&!
						      02.930D-6,&!
						      00.975D-6 /)! 3.9 um
  real*8, DIMENSION (10) :: bstl_ash
  integer iash
!srf

!
! bstl is for budgets
!

  real*8 conver,converi
       conver=1.e-9
       converi=1.e9
       lmx=kte-kts+1
       do j=jts,jte
       do i=its,ite
          kk=0
	  bstl_ash(:)=0.
	  ash(:,:,:,:)=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)
          airmas(1,1,kk)=-(p8w(i,k+1,j)-p8w(i,k,j))/g
          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

!ash settling

          iseas=0
          idust=0
	  iash =1
	  
          kk=0
!         write(0,*)'1',chem(i,1,j,p_dust_4)
          do k=kts,kte
          kk=kk+1
          ash(1,1,kk,7)=chem(i,k,j,p_vash_7)*conver
          ash(1,1,kk,8)=chem(i,k,j,p_vash_8)*conver
          ash(1,1,kk,9)=chem(i,k,j,p_vash_9)*conver
          ash(1,1,kk,10)=chem(i,k,j,p_vash_10)*conver
          enddo
          if(config_flags%chem_opt == 400  .or. config_flags%chem_opt == 402 ) then
          do k=kts,kte
          kk=kk+1
          ash(1,1,kk,1)=chem(i,k,j,p_vash_1)*conver
          ash(1,1,kk,2)=chem(i,k,j,p_vash_2)*conver
          ash(1,1,kk,3)=chem(i,k,j,p_vash_3)*conver
          ash(1,1,kk,4)=chem(i,k,j,p_vash_4)*conver
          ash(1,1,kk,5)=chem(i,k,j,p_vash_5)*conver
          ash(1,1,kk,6)=chem(i,k,j,p_vash_6)*conver
          enddo
          endif

          call vsettling(1, 1, lmx, 10, g,&
                    ash, tmp, p_mid, delz, airmas, &
                    den_ash, reff_ash, dt, bstl_ash, rh, idust, iseas,iash)
          kk=0
          ash_fall(i,j)=ash_fall(i,j)+sum(bstl_ash(1:10))
          do k=kts,kte
          kk=kk+1
            chem(i,k,j,p_vash_7)=ash(1,1,kk,7)*converi
            chem(i,k,j,p_vash_8)=ash(1,1,kk,8)*converi
            chem(i,k,j,p_vash_9)=ash(1,1,kk,9)*converi
            chem(i,k,j,p_vash_10)=ash(1,1,kk,10)*converi
          enddo
          if(config_flags%chem_opt == 400  .or. config_flags%chem_opt == 402 ) then
          do k=kts,kte
          kk=kk+1
            chem(i,k,j,p_vash_1)=ash(1,1,kk,1)*converi
            chem(i,k,j,p_vash_2)=ash(1,1,kk,2)*converi
            chem(i,k,j,p_vash_3)=ash(1,1,kk,3)*converi
            chem(i,k,j,p_vash_4)=ash(1,1,kk,4)*converi
            chem(i,k,j,p_vash_5)=ash(1,1,kk,5)*converi
            chem(i,k,j,p_vash_6)=ash(1,1,kk,6)*converi
          enddo
          endif

!ash settling end

       enddo
       enddo
END SUBROUTINE vash_settling_driver


          subroutine vsettling(imx,jmx, lmx, nmx,g0, &
                    tc, tmp, p_mid, delz, airmas, &
                    den, reff, dt, bstl, rh, idust, iseas,iash)
! ****************************************************************************
! *                                                                          *
! *  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,iash
  INTEGER :: ntdt
  REAL, INTENT(IN) :: dt,g0 ! ,dyn_visc
  REAL*8,    INTENT(IN) :: tmp(imx,jmx,lmx), delz(imx,jmx,lmx),  &
                         airmas(imx,jmx,lmx), rh(imx,jmx,lmx), &
                         den(nmx), reff(nmx), p_mid(imx,jmx,lmx)
  REAL*8, INTENT(INOUT) :: tc(imx,jmx,lmx,nmx)
  REAL*8, INTENT(OUT)   :: bstl(imx,jmx,nmx)

  REAL*8    :: tc1(imx,jmx,lmx,nmx), 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    :: addmass,c_stokes, free_path, c_cun, viscosity, vd_cor, growth_fac
  REAL,    PARAMETER :: dyn_visc = 1.5E-5
  INTEGER :: k, n, i, j, l, l2
  ! 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)

  ! executable statements

! IF (type) /= 'dust' .AND. TRIM(aero_type) /= 'sea_salt') RETURN
  if(idust.ne.1.and.iseas.ne.1.and.iash.ne.1)return

  WHERE (tc(:,:,:,:) < 0.0) tc(:,:,:,:) = 1.0d-32

  dzmin = MINVAL(delz(:,:,:))
  IF (idust == 1)     growth_fac = 1.0
  IF (iseas == 1)     growth_fac = 3.0
  IF (iash  == 1)     growth_fac = 1.0

  DO k = 1,nmx

     ! 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

     tc1(:,:,:,k) = tc(:,:,:,k)
     vsettl = 2.0/9.0 * g0 * den(k) * (growth_fac*reff(k))**2 / &
              (0.5*dyn_visc)

     ! Determine the maximum time-step satisying the CFL condition:
     ! 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
!srf     IF (idust.eq.1)then
     IF (idust.eq.1   .or. iash==1)then
          rwet(k) = reff(k)
          ratio_r(k) = 1.0
          rho(k) = den(k)
      endif
  END DO

  ! 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
 
     DO k = 1,nmx
        IF (idust.eq.1 .or. iash==1) THEN
           rwet_priv(k) = rwet(k)
           rho_priv(k)  = rho(k)
        END IF

        DO n = 1,ndt_settl(k)

           ! Solve each vertical layer successively (layer l)
      
           DO l = lmx,1,-1
              l2 = lmx - l + 1

!           DO j = 1,jmx
              DO i = 1,imx

                 ! Dynamic viscosity
                 c_stokes = 1.458E-6 * tmp(i,j,l)**1.5/(tmp(i,j,l) + 110.4) 

                 ! Mean free path as a function of pressure (mb) and 
                 ! temperature (K)
                 ! order of p_mid is top->sfc
                 free_path = 1.1E-3/p_mid(i,j,l2)/SQRT(tmp(i,j,l))
!!!                 free_path = 1.1E-3/p_edge(i,j,l2)/SQRT(tmp(i,j,l))

                 ! Slip Correction Factor
                 c_cun = 1.0+ free_path/rwet_priv(k)* &
                      (1.257 + 0.4*EXP(-1.1*rwet_priv(k)/free_path))

                 ! Corrected dynamic viscosity (kg/m/s)
                 viscosity = c_stokes / c_cun

                 ! Settling velocity
!                IF (iseas.eq.1) THEN
!                   rho_priv(k) = ratio_r(k)*den(k) + (1.0 - ratio_r(k))*1000.0
!                END IF

                 vd_cor = 2.0/9.0*g0*rho_priv(k)*rwet_priv(k)**2/viscosity

                 ! Update mixing ratio
                 ! Order of delz is top->sfc
                 IF (l == lmx) THEN
                    tc(i,j,l,k) = tc(i,j,l,k) / &
                         (1.0 + dt_settl(k)*vd_cor/delz(i,j,l2))
                 ELSE
                    tc(i,j,l,k) = 1.0/(1.0+dt_settl(k)*vd_cor/delz(i,j,l2))&
                         *(tc(i,j,l,k) + dt_settl(k)*vd_cor /delz(i,j,l2-1) &
                         * tc(i,j,l+1,k))
                 END IF
              END DO   !i
!           END DO   !j
        END DO  !l

     END DO  !n
  END DO  !k

  END DO   !j
!$OMP END PARALLEL DO

  DO n = 1,nmx
     DO i = 1,imx
        DO j = 1,jmx
           bstl(i,j,n) = 0.0
           addmass=0.
           DO l = 1,lmx
              addmass=addmass+(tc(i,j,l,n) - tc1(i,j,l,n)) * airmas(i,j,l)
              IF (tc(i,j,l,n) < 0.0) tc(i,j,l,n) = 1.0D-32
           END DO
           if(addmass.gt.0.)addmass=0
           bstl(i,j,n) = bstl(i,j,n) - addmass
        END DO
     END DO
  END DO
  
END SUBROUTINE vsettling

END MODULE MODULE_VASH_SETTLING
