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 kk=0 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 kk=0 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