Line data    Source code

       1             : module aer_drydep_mod
       2             : 
       3             :   use shr_kind_mod, only: r8 => shr_kind_r8
       4             :   use ppgrid
       5             : 
       6             :       ! Shared Data for dry deposition calculation.
       7             : 
       8             :       real(r8) rair                ! Gas constant for dry air (J/K/kg)
       9             :       real(r8) gravit              ! Gravitational acceleration
      10             : !      real(r8), allocatable :: phi(:)           ! grid latitudes (radians)11
      11             : 
      12             : contains
      13             : 
      14             : !##############################################################################
      15             : 
      16             : ! $Id$
      17             : 
      18           0 :       subroutine inidrydep( xrair, xgravit) !, xphi )
      19             : 
      20             : ! Initialize dry deposition parameterization.
      21             : 
      22             :       implicit none
      23             : 
      24             : ! Input arguments:
      25             :       real(r8), intent(in) :: xrair                ! Gas constant for dry air
      26             :       real(r8), intent(in) :: xgravit              ! Gravitational acceleration
      27             : !      real(r8), intent(in) :: xphi(:)           ! grid latitudes (radians)
      28             : 
      29             : ! Local variables:
      30             :       integer i, j, ncid, vid, ns
      31             : !-----------------------------------------------------------------------
      32             : !      ns = size(xphi)
      33             : !      allocate(phi(ns))
      34           0 :       rair = xrair
      35           0 :       gravit = xgravit
      36             : !      do j = 1, ns
      37             : !         phi(j) = xphi(j)
      38             : !      end do
      39             : 
      40           0 :       return
      41             :       end subroutine inidrydep
      42             : 
      43             : !##############################################################################
      44             : 
      45           0 :       subroutine setdvel( ncol, landfrac, icefrac, ocnfrac, vgl, vgo, vgsi, vg )
      46             : 
      47             : ! Set the deposition velocity depending on whether we are over
      48             : ! land, ocean, and snow/ice
      49             : 
      50             : 
      51             :       implicit none
      52             : 
      53             : ! Input arguments:
      54             : 
      55             :       integer, intent(in) :: ncol
      56             :       real (r8), intent(in) :: landfrac(pcols)       ! land fraction
      57             :       real (r8), intent(in) :: icefrac(pcols)       ! ice fraction
      58             :       real (r8), intent(in) :: ocnfrac(pcols)       ! ocean fraction
      59             : 
      60             :       real(r8), intent(in) :: vgl                  ! dry deposition velocity in m/s (land)
      61             :       real(r8), intent(in) :: vgo                  ! dry deposition velocity in m/s (ocean)
      62             :       real(r8), intent(in) :: vgsi                 ! dry deposition velocity in m/s (snow/ice)
      63             : 
      64             : ! Output arguments:
      65             :       real(r8), intent(out) ::  vg(pcols) ! dry deposition velocity in m/s
      66             : 
      67             : ! Local variables:
      68             : 
      69             :       integer i
      70             :       real(r8) a
      71             : 
      72             : 
      73           0 :       do i = 1, ncol
      74           0 :          vg(i) = landfrac(i)*vgl + ocnfrac(i)*vgo + icefrac(i)*vgsi
      75             : !         if (ioro(i).eq.0) then
      76             : !            vg(i) = vgo
      77             : !         else if (ioro(i).eq.1) then
      78             : !            vg(i) = vgl
      79             : !         else
      80             : !            vg(i) = vgsi
      81             : !         endif
      82             :       end do
      83             : 
      84           0 :       return
      85             :       end subroutine setdvel
      86             : 
      87             : !##############################################################################
      88             : 
      89           0 :       subroutine ddflux( ncol, vg, q, p, tv, flux )
      90             : 
      91             : ! Compute surface flux due to dry deposition processes.
      92             : 
      93             : 
      94             :       implicit none
      95             : 
      96             : ! Input arguments:
      97             :       integer , intent(in) :: ncol
      98             :       real(r8), intent(in) ::    vg(pcols)  ! dry deposition velocity in m/s
      99             :       real(r8), intent(in) ::    q(pcols)   ! tracer conc. in surface layer (kg tracer/kg moist air)
     100             :       real(r8), intent(in) ::    p(pcols)   ! midpoint pressure in surface layer (Pa)
     101             :       real(r8), intent(in) ::    tv(pcols)  ! midpoint virtual temperature in surface layer (K)
     102             : 
     103             : ! Output arguments:
     104             : 
     105             :       real(r8), intent(out) ::    flux(pcols) ! flux due to dry deposition in kg/m^s/sec
     106             : 
     107             : ! Local variables:
     108             : 
     109             :       integer i
     110             : 
     111           0 :       do i = 1, ncol
     112           0 :          flux(i) = -vg(i) * q(i) * p(i) /(tv(i) * rair)
     113             :       end do
     114             : 
     115           0 :       return
     116             :       end subroutine ddflux
     117             : 
     118             : !------------------------------------------------------------------------
     119             : !BOP
     120             : !
     121             : ! !IROUTINE: subroutine d3ddflux
     122             : !
     123             : ! !INTERFACE:
     124             : !
     125           0 :    subroutine  d3ddflux ( ncol, vlc_dry, q,pmid,pdel, tv, dep_dry,dep_dry_tend,dt)
     126             : ! Description:
     127             : !Do 3d- settling deposition calculations following Zender's dust codes, Dec 02.
     128             : !
     129             : ! Author: Natalie Mahowald
     130             : !
     131             :       implicit none
     132             : 
     133             : ! Input arguments:
     134             :       integer , intent(in) :: ncol
     135             :       real(r8), intent(in) ::    vlc_dry(pcols,pver)  ! dry deposition velocity in m/s
     136             :       real(r8), intent(in) ::    q(pcols,pver)   ! tracer conc. in surface layer (kg tracer/kg moist air)
     137             :       real(r8), intent(in) ::    pmid(pcols,pver)   ! midpoint pressure in surface layer (Pa)
     138             :       real(r8), intent(in) ::    pdel(pcols,pver)   ! delta pressure across level (Pa)
     139             :       real(r8), intent(in) ::    tv(pcols,pver)  ! midpoint virtual temperature in surface layer (K)
     140             :     real(r8),            intent(in)  :: dt             ! time step
     141             : 
     142             : ! Output arguments:
     143             : 
     144             :       real(r8), intent(out) ::    dep_dry(pcols) ! flux due to dry deposition in kg /m^s/sec
     145             :       real(r8), intent(out) ::    dep_dry_tend(pcols,pver) ! flux due to dry deposition in kg /m^s/sec
     146             : 
     147             : ! Local variables:
     148             : 
     149             :       real(r8) :: flux(pcols,0:pver)  ! downward flux at each level:  kg/m2/s 
     150             :       integer i,k
     151           0 :       do i=1,ncol
     152           0 :          flux(i,0)=0._r8
     153             :       enddo
     154           0 :       do k=1,pver
     155           0 :          do i = 1, ncol
     156           0 :             flux(i,k) = -min(vlc_dry(i,k) * q(i,k) * pmid(i,k) /(tv(i,k) * rair), &
     157           0 :                       q(i,k)*pdel(i,k)/gravit/dt)
     158           0 :             dep_dry_tend(i,k)=(flux(i,k)-flux(i,k-1))/pdel(i,k)*gravit  !kg/kg/s
     159             : 
     160             :          end do
     161             :       enddo
     162             : ! surface flux:
     163           0 :       do i=1,ncol
     164           0 :          dep_dry(i)=flux(i,pver)
     165             :       enddo
     166           0 :       return
     167             :       end subroutine d3ddflux
     168             : 
     169             : 
     170             : 
     171             : !------------------------------------------------------------------------
     172             : !BOP
     173             : !
     174             : ! !IROUTINE: subroutine Calcram
     175             : !
     176             : ! !INTERFACE:
     177             : !
     178             : 
     179           0 :       subroutine  calcram(ncol,landfrac,icefrac,ocnfrac,obklen,&
     180             :            ustar,ram1in,ram1,t,pmid,&
     181             :            pdel,fvin,fv)
     182             :         !
     183             :         ! !DESCRIPTION: 
     184             :         !  
     185             :         ! Calc aerodynamic resistance over oceans and sea ice (comes in from land model)
     186             :         ! from Seinfeld and Pandis, p.963.
     187             :         !  
     188             :         ! Author: Natalie Mahowald
     189             :         !
     190             :         implicit none
     191             :         integer, intent(in) :: ncol
     192             :         real(r8),intent(in) :: ram1in(pcols)         !aerodynamical resistance (s/m)
     193             :         real(r8),intent(in) :: fvin(pcols)                 ! sfc frc vel from land
     194             :         real(r8),intent(out) :: ram1(pcols)         !aerodynamical resistance (s/m)
     195             :         real(r8),intent(out) :: fv(pcols)                 ! sfc frc vel from land
     196             :         real(r8), intent(in) :: obklen(pcols)                 ! obklen
     197             :         real(r8), intent(in) :: ustar(pcols)                  ! sfc fric vel
     198             :         real(r8), intent(in) :: landfrac(pcols)               ! land fraction
     199             :         real(r8), intent(in) :: icefrac(pcols)                ! ice fraction
     200             :         real(r8), intent(in) :: ocnfrac(pcols)                ! ocean fraction
     201             :         real(r8), intent(in) :: t(pcols)       !atm temperature (K)
     202             :         real(r8), intent(in) :: pmid(pcols)    !atm pressure (Pa)
     203             :         real(r8), intent(in) :: pdel(pcols)    !atm pressure (Pa)
     204             :         real(r8), parameter :: zzocen = 0.0001_r8   ! Ocean aerodynamic roughness length
     205             :         real(r8), parameter :: zzsice = 0.0400_r8   ! Sea ice aerodynamic roughness length
     206             :         real(r8), parameter :: xkar   = 0.4_r8      ! Von Karman constant
     207             : 
     208             :         ! local variables
     209             :         real(r8) :: z,psi,psi0,nu,nu0,temp,ram
     210             :         integer :: i
     211             :         !    write(iulog,*) rair,zzsice,zzocen,gravit,xkar
     212             : 
     213             : 
     214           0 :         do i=1,ncol
     215           0 :            z=pdel(i)*rair*t(i)/pmid(i)/gravit/2.0_r8   !use half the layer height like Ganzefeld and Lelieveld, 1995
     216           0 :            if(obklen(i).eq.0) then
     217             :               psi=0._r8
     218             :               psi0=0._r8
     219             :            else
     220           0 :               psi=min(max(z/obklen(i),-1.0_r8),1.0_r8)
     221           0 :               psi0=min(max(zzocen/obklen(i),-1.0_r8),1.0_r8)
     222             :            endif
     223           0 :            temp=z/zzocen
     224           0 :            if(icefrac(i) > 0.5_r8) then 
     225           0 :               if(obklen(i).gt.0) then 
     226           0 :                  psi0=min(max(zzsice/obklen(i),-1.0_r8),1.0_r8)
     227             :               else
     228             :                  psi0=0.0_r8
     229             :               endif
     230           0 :               temp=z/zzsice
     231             :            endif
     232           0 :            if(psi> 0._r8) then
     233           0 :               ram=1/xkar/ustar(i)*(log(temp)+4.7_r8*(psi-psi0))
     234             :            else
     235           0 :               nu=(1.00_r8-15.000_r8*psi)**(.25_r8)
     236           0 :               nu0=(1.000_r8-15.000_r8*psi0)**(.25_r8)
     237           0 :               if(ustar(i).ne.0._r8) then
     238             :                  ram=1/xkar/ustar(i)*(log(temp) &
     239             :                       +log(((nu0**2+1.00_r8)*(nu0+1.0_r8)**2)/((nu**2+1.0_r8)*(nu+1.00_r8)**2)) &
     240           0 :                       +2.0_r8*(atan(nu)-atan(nu0)))
     241             :               else
     242             :                  ram=0._r8
     243             :               endif
     244             :            endif
     245           0 :            if(landfrac(i) < 0.000000001_r8) then
     246           0 :               fv(i)=ustar(i)
     247           0 :               ram1(i)=ram
     248             :            else
     249           0 :               fv(i)=fvin(i)
     250           0 :               ram1(i)=ram1in(i)
     251             :            endif
     252             :            !          write(iulog,*) i,pdel(i),t(i),pmid(i),gravit,obklen(i),psi,psi0,icefrac(i),nu,nu0,ram,ustar(i),&
     253             :            !             log(((nu0**2+1.00)*(nu0+1.0)**2)/((nu**2+1.0)*(nu+1.00)**2)),2.0*(atan(nu)-atan(nu0))
     254             : 
     255             :         enddo
     256             : 
     257             :         ! fvitt -- fv == 0 causes a floating point exception in 
     258             :         ! dry dep of sea salts and dust
     259           0 :         where ( fv(:ncol) == 0._r8 ) 
     260           0 :            fv(:ncol) = 1.e-12_r8
     261             :         endwhere
     262             : 
     263           0 :         return
     264             :       end subroutine calcram
     265             : 
     266             : 
     267             : !##############################################################################
     268             : end module aer_drydep_mod