LCOV - code coverage report
Current view: | top level - chemistry/aerosol - dust_sediment_mod.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 118 | 0.0 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 0 | 5 | 0.0 %

              Line data    Source code
           1             : module dust_sediment_mod
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : ! Purpose:
           5             : !
           6             : ! Contains routines to compute tendencies from sedimentation of dust
           7             : !
           8             : ! Author: Phil Rasch
           9             : !
          10             : !---------------------------------------------------------------------------------
          11             :
          12             :   use shr_kind_mod,      only: r8=>shr_kind_r8
          13             :   use ppgrid,            only: pcols, pver, pverp
          14             :   use physconst,         only: gravit, rair
          15             :   use cam_logfile,       only: iulog
          16             :   use cam_abortutils,    only: endrun
          17             :
          18             :   private
          19             :   public :: dust_sediment_vel, dust_sediment_tend
          20             :
          21             :
          22             :   real (r8), parameter :: vland  = 2.8_r8            ! dust fall velocity over land  (cm/s)
          23             :   real (r8), parameter :: vocean = 1.5_r8            ! dust fall velocity over ocean (cm/s)
          24             :   real (r8), parameter :: mxsedfac   = 0.99_r8       ! maximum sedimentation flux factor
          25             :
          26             : contains
          27             :
          28             : !===============================================================================
          29           0 :   subroutine dust_sediment_vel (ncol,                               &
          30             :        icefrac , landfrac, ocnfrac , pmid    , pdel    , t       , &
          31             :        dustmr  , pvdust   )
          32             :
          33             : !----------------------------------------------------------------------
          34             :
          35             : ! Compute gravitational sedimentation velocities for dust
          36             :
          37             :     implicit none
          38             :
          39             : ! Arguments
          40             :     integer, intent(in) :: ncol                     ! number of colums to process
          41             :
          42             :     real(r8), intent(in)  :: icefrac (pcols)        ! sea ice fraction (fraction)
          43             :     real(r8), intent(in)  :: landfrac(pcols)        ! land fraction (fraction)
          44             :     real(r8), intent(in)  :: ocnfrac (pcols)        ! ocean fraction (fraction)
          45             :     real(r8), intent(in)  :: pmid  (pcols,pver)     ! pressure of midpoint levels (Pa)
          46             :     real(r8), intent(in)  :: pdel  (pcols,pver)     ! pressure diff across layer (Pa)
          47             :     real(r8), intent(in)  :: t     (pcols,pver)     ! temperature (K)
          48             :     real(r8), intent(in)  :: dustmr(pcols,pver)     ! dust (kg/kg)
          49             :
          50             :     real(r8), intent(out) :: pvdust (pcols,pverp)    ! vertical velocity of dust (Pa/s)
          51             : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
          52             :
          53             : ! Local variables
          54             :     real (r8) :: rho(pcols,pver)                    ! air density in kg/m3
          55             :     real (r8) :: vfall(pcols)                       ! settling velocity of dust particles (m/s)
          56             :
          57             :     integer i,k
          58             :
          59             :     real (r8) :: lbound, ac, bc, cc
          60             :
          61             : !-----------------------------------------------------------------------
          62             : !--------------------- dust fall velocity ----------------------------
          63             : !-----------------------------------------------------------------------
          64             :
          65           0 :     do k = 1,pver
          66           0 :        do i = 1,ncol
          67             :
          68             :           ! merge the dust fall velocities for land and ocean (cm/s)
          69             :           ! SHOULD ALSO ACCOUNT FOR ICEFRAC
          70           0 :           vfall(i) = vland*landfrac(i) + vocean*(1._r8-landfrac(i))
          71             :           !!         vfall(i) = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)
          72             :
          73             :           ! fall velocity (assume positive downward)
          74           0 :           pvdust(i,k+1) = vfall(i)
          75             :        end do
          76             :     end do
          77             :
          78           0 :     return
          79             :   end subroutine dust_sediment_vel
          80             :
          81             :
          82             : !===============================================================================
          83           0 :   subroutine dust_sediment_tend ( &
          84             :        ncol,   dtime,  pint,     pmid,    pdel,  t,   &
          85             :        dustmr ,pvdust, dusttend, sfdust )
          86             :
          87             : !----------------------------------------------------------------------
          88             : !     Apply Particle Gravitational Sedimentation
          89             : !----------------------------------------------------------------------
          90             :
          91             :     implicit none
          92             :
          93             : ! Arguments
          94             :     integer,  intent(in)  :: ncol                      ! number of colums to process
          95             :
          96             :     real(r8), intent(in)  :: dtime                     ! time step
          97             :     real(r8), intent(in)  :: pint  (pcols,pverp)       ! interfaces pressure (Pa)
          98             :     real(r8), intent(in)  :: pmid  (pcols,pver)        ! midpoint pressures (Pa)
          99             :     real(r8), intent(in)  :: pdel  (pcols,pver)        ! pressure diff across layer (Pa)
         100             :     real(r8), intent(in)  :: t     (pcols,pver)        ! temperature (K)
         101             :     real(r8), intent(in)  :: dustmr(pcols,pver)        ! dust (kg/kg)
         102             :     real(r8), intent(in)  :: pvdust (pcols,pverp)      ! vertical velocity of dust drops  (Pa/s)
         103             : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
         104             :
         105             :     real(r8), intent(out) :: dusttend(pcols,pver)      ! dust tend
         106             :     real(r8), intent(out) :: sfdust  (pcols)           ! surface flux of dust (rain, kg/m/s)
         107             :
         108             : ! Local variables
         109             :     real(r8) :: fxdust(pcols,pverp)                     ! fluxes at the interfaces, dust (positive = down)
         110             :
         111             :     integer :: i,k
         112             : !----------------------------------------------------------------------
         113             :
         114             : ! initialize variables
         115           0 :     fxdust  (:ncol,:) = 0._r8 ! flux at interfaces (dust)
         116           0 :     dusttend(:ncol,:) = 0._r8 ! tend (dust)
         117           0 :     sfdust(:ncol)     = 0._r8 ! sedimentation flux out bot of column (dust)
         118             :
         119             : ! fluxes at interior points
         120           0 :     call getflx(ncol, pint, dustmr, pvdust, dtime, fxdust)
         121             :
         122             : ! calculate fluxes at boundaries
         123           0 :     do i = 1,ncol
         124           0 :        fxdust(i,1) = 0
         125             : ! surface flux by upstream scheme
         126           0 :        fxdust(i,pverp) = dustmr(i,pver) * pvdust(i,pverp) * dtime
         127             :     end do
         128             :
         129             : ! filter out any negative fluxes from the getflx routine
         130           0 :     do k = 2,pver
         131           0 :        fxdust(:ncol,k) = max(0._r8, fxdust(:ncol,k))
         132             :     end do
         133             :
         134             : ! Limit the flux out of the bottom of each cell to the water content in each phase.
         135             : ! Apply mxsedfac to prevent generating very small negative cloud water/ice
         136             : ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
         137             : ! ***Should we include the flux in the top, to allow for thin surface layers?
         138             : ! ***Requires simple treatment of cloud overlap, already included below.
         139           0 :     do k = 1,pver
         140           0 :        do i = 1,ncol
         141           0 :           fxdust(i,k+1) = min( fxdust(i,k+1), mxsedfac * dustmr(i,k) * pdel(i,k) )
         142             : !!$        fxdust(i,k+1) = min( fxdust(i,k+1), dustmr(i,k) * pdel(i,k) + fxdust(i,k))
         143             :        end do
         144             :     end do
         145             :
         146             : ! Now calculate the tendencies
         147           0 :     do k = 1,pver
         148           0 :        do i = 1,ncol
         149             : ! net flux into cloud changes cloud dust/ice (all flux is out of cloud)
         150           0 :           dusttend(i,k)  = (fxdust(i,k) - fxdust(i,k+1)) / (dtime * pdel(i,k))
         151             :        end do
         152             :     end do
         153             :
         154             : ! convert flux out the bottom to mass units Pa -> kg/m2/s
         155           0 :     sfdust(:ncol) = fxdust(:ncol,pverp) / (dtime*gravit)
         156             :
         157           0 :     return
         158             :   end subroutine dust_sediment_tend
         159             :
         160             : !===============================================================================
         161           0 :   subroutine getflx(ncol, xw, phi, vel, deltat, flux)
         162             :
         163             : !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
         164             : !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
         165             : !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
         166             : !.........phi1......phi2.......phi3.....phi4.......phi5.......
         167             :
         168             :
         169             :     implicit none
         170             :
         171             :     integer ncol                      ! number of colums to process
         172             :
         173             :     integer i
         174             :     integer k
         175             :
         176             :     real (r8) vel(pcols,pverp)
         177             :     real (r8) flux(pcols,pverp)
         178             :     real (r8) xw(pcols,pverp)
         179             :     real (r8) psi(pcols,pverp)
         180             :     real (r8) phi(pcols,pverp-1)
         181             :     real (r8) fdot(pcols,pverp)
         182             :     real (r8) xx(pcols)
         183             :     real (r8) fxdot(pcols)
         184             :     real (r8) fxdd(pcols)
         185             :
         186             :     real (r8) psistar(pcols)
         187             :     real (r8) deltat
         188             :
         189             :     real (r8) xxk(pcols,pver)
         190             :
         191           0 :     do i = 1,ncol
         192             : !        integral of phi
         193           0 :        psi(i,1) = 0._r8
         194             : !        fluxes at boundaries
         195           0 :        flux(i,1) = 0
         196           0 :        flux(i,pverp) = 0._r8
         197             :     end do
         198             :
         199             : !     integral function
         200           0 :     do k = 2,pverp
         201           0 :        do i = 1,ncol
         202           0 :           psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
         203             :        end do
         204             :     end do
         205             :
         206             :
         207             : !     calculate the derivatives for the interpolating polynomial
         208           0 :     call cfdotmc_pro (ncol, xw, psi, fdot)
         209             :
         210             : !  NEW WAY
         211             : !     calculate fluxes at interior pts
         212           0 :     do k = 2,pver
         213           0 :        do i = 1,ncol
         214           0 :           xxk(i,k) = xw(i,k)-vel(i,k)*deltat
         215             :        end do
         216             :     end do
         217           0 :     do k = 2,pver
         218           0 :        call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
         219           0 :        do i = 1,ncol
         220           0 :           flux(i,k) = (psi(i,k)-psistar(i))
         221             :        end do
         222             :     end do
         223             :
         224             :
         225           0 :     return
         226             :   end subroutine getflx
         227             :
         228             :
         229             :
         230             : !##############################################################################
         231             :
         232           0 :   subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar)
         233             :
         234             :
         235             :     implicit none
         236             :
         237             : ! input
         238             :     integer ncol                      ! number of colums to process
         239             :
         240             :     real (r8) x(pcols, pverp)
         241             :     real (r8) f(pcols, pverp)
         242             :     real (r8) fdot(pcols, pverp)
         243             :     real (r8) xin(pcols)
         244             :
         245             : ! output
         246             :     real (r8) fxdot(pcols)
         247             :     real (r8) fxdd(pcols)
         248             :     real (r8) psistar(pcols)
         249             :
         250             :     integer i
         251             :     integer k
         252             :     integer intz(pcols)
         253             :     real (r8) dx
         254             :     real (r8) s
         255             :     real (r8) c2
         256             :     real (r8) c3
         257             :     real (r8) xx
         258             :     real (r8) xinf
         259             :     real (r8) psi1, psi2, psi3, psim
         260             :     real (r8) cfint
         261             :     real (r8) cfnew
         262             :     real (r8) xins(pcols)
         263             :
         264             : !     the minmod function
         265             :     real (r8) a, b, c
         266             :     real (r8) minmod
         267             :     real (r8) medan
         268             :     minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
         269             :     medan(a,b,c) = a + minmod(b-a,c-a)
         270             :
         271           0 :     do i = 1,ncol
         272           0 :        xins(i) = medan(x(i,1), xin(i), x(i,pverp))
         273           0 :        intz(i) = 0
         274             :     end do
         275             :
         276             : ! first find the interval
         277           0 :     do k =  1,pverp-1
         278           0 :        do i = 1,ncol
         279           0 :           if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0._r8) then
         280           0 :              intz(i) = k
         281             :           endif
         282             :        end do
         283             :     end do
         284             :
         285           0 :     do i = 1,ncol
         286           0 :        if (intz(i).eq.0) then
         287           0 :           write(iulog,*) ' interval was not found for col i ', i
         288           0 :           call endrun('DUST_SEDIMENT_MOD:cfint2 -- interval was not found ')
         289             :        endif
         290             :     end do
         291             :
         292             : ! now interpolate
         293           0 :     do i = 1,ncol
         294           0 :        k = intz(i)
         295           0 :        dx = (x(i,k+1)-x(i,k))
         296           0 :        s = (f(i,k+1)-f(i,k))/dx
         297           0 :        c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
         298           0 :        c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
         299           0 :        xx = (xins(i)-x(i,k))
         300           0 :        fxdot(i) =  (3*c3*xx + 2*c2)*xx + fdot(i,k)
         301           0 :        fxdd(i) = 6*c3*xx + 2*c2
         302           0 :        cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)
         303             :
         304             : !        limit the interpolant
         305           0 :        psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
         306           0 :        if (k.eq.1) then
         307           0 :           psi2 = f(i,1)
         308             :        else
         309           0 :           psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
         310             :        endif
         311           0 :        if (k+1.eq.pverp) then
         312           0 :           psi3 = f(i,pverp)
         313             :        else
         314           0 :           psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
         315             :        endif
         316           0 :        psim = medan(psi1, psi2, psi3)
         317           0 :        cfnew = medan(cfint, psi1, psim)
         318             :        if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36_r8)  .gt..03_r8) then
         319             : !     CHANGE THIS BACK LATER!!!
         320             : !     $        .gt..1) then
         321             :
         322             :
         323             : !     UNCOMMENT THIS LATER!!!
         324             : !            write(iulog,*) ' cfint2 limiting important ', cfint, cfnew
         325             :
         326             :
         327             :        endif
         328           0 :        psistar(i) = cfnew
         329             :     end do
         330             :
         331           0 :     return
         332             :   end subroutine cfint2
         333             :
         334             :
         335             :
         336             : !##############################################################################
         337             :
         338           0 :   subroutine cfdotmc_pro (ncol, x, f, fdot)
         339             :
         340             : !     prototype version; eventually replace with final SPITFIRE scheme
         341             :
         342             : !     calculate the derivative for the interpolating polynomial
         343             : !     multi column version
         344             :
         345             :
         346             :     implicit none
         347             :
         348             : ! input
         349             :     integer ncol                      ! number of colums to process
         350             :
         351             :     real (r8) x(pcols, pverp)
         352             :     real (r8) f(pcols, pverp)
         353             : ! output
         354             :     real (r8) fdot(pcols, pverp)          ! derivative at nodes
         355             :
         356             : ! assumed variable distribution
         357             : !     x1.......x2.......x3.......x4.......x5.......x6     1,pverp points
         358             : !     f1.......f2.......f3.......f4.......f5.......f6     1,pverp points
         359             : !     ...sh1.......sh2......sh3......sh4......sh5....     1,pver points
         360             : !     .........d2.......d3.......d4.......d5.........     2,pver points
         361             : !     .........s2.......s3.......s4.......s5.........     2,pver points
         362             : !     .............dh2......dh3......dh4.............     2,pver-1 points
         363             : !     .............eh2......eh3......eh4.............     2,pver-1 points
         364             : !     ..................e3.......e4..................     3,pver-1 points
         365             : !     .................ppl3......ppl4................     3,pver-1 points
         366             : !     .................ppr3......ppr4................     3,pver-1 points
         367             : !     .................t3........t4..................     3,pver-1 points
         368             : !     ................fdot3.....fdot4................     3,pver-1 points
         369             :
         370             :
         371             : ! work variables
         372             :
         373             :
         374             :     integer i
         375             :     integer k
         376             :
         377             :     real (r8) a                    ! work var
         378             :     real (r8) b                    ! work var
         379             :     real (r8) c                    ! work var
         380             :     real (r8) s(pcols,pverp)             ! first divided differences at nodes
         381             :     real (r8) sh(pcols,pverp)            ! first divided differences between nodes
         382             :     real (r8) d(pcols,pverp)             ! second divided differences at nodes
         383             :     real (r8) dh(pcols,pverp)            ! second divided differences between nodes
         384             :     real (r8) e(pcols,pverp)             ! third divided differences at nodes
         385             :     real (r8) eh(pcols,pverp)            ! third divided differences between nodes
         386             :     real (r8) pp                   ! p prime
         387             :     real (r8) ppl(pcols,pverp)           ! p prime on left
         388             :     real (r8) ppr(pcols,pverp)           ! p prime on right
         389             :     real (r8) qpl
         390             :     real (r8) qpr
         391             :     real (r8) ttt
         392             :     real (r8) t
         393             :     real (r8) tmin
         394             :     real (r8) tmax
         395             :     real (r8) delxh(pcols,pverp)
         396             :
         397             :
         398             : !     the minmod function
         399             :     real (r8) minmod
         400             :     real (r8) medan
         401             :     minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
         402             :     medan(a,b,c) = a + minmod(b-a,c-a)
         403             :
         404           0 :     do k = 1,pver
         405             :
         406             :
         407             : !        first divided differences between nodes
         408           0 :        do i = 1, ncol
         409           0 :           delxh(i,k) = (x(i,k+1)-x(i,k))
         410           0 :           sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
         411             :        end do
         412             :
         413             : !        first and second divided differences at nodes
         414           0 :        if (k.ge.2) then
         415           0 :           do i = 1,ncol
         416           0 :              d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
         417           0 :              s(i,k) = minmod(sh(i,k),sh(i,k-1))
         418             :           end do
         419             :        endif
         420             :     end do
         421             :
         422             : !     second and third divided diffs between nodes
         423           0 :     do k = 2,pver-1
         424           0 :        do i = 1, ncol
         425           0 :           eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
         426           0 :           dh(i,k) = minmod(d(i,k),d(i,k+1))
         427             :        end do
         428             :     end do
         429             :
         430             : !     treat the boundaries
         431           0 :     do i = 1,ncol
         432           0 :        e(i,2) = eh(i,2)
         433           0 :        e(i,pver) = eh(i,pver-1)
         434             : !        outside level
         435             :        fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1)  &
         436           0 :             - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
         437           0 :        fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
         438             :        fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver)  &
         439           0 :             + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
         440           0 :        fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
         441             : !        one in from boundary
         442           0 :        fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
         443           0 :        fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
         444             :        fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver)   &
         445           0 :             - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
         446           0 :        fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
         447             :     end do
         448             :
         449             :
         450           0 :     do k = 3,pver-1
         451           0 :        do i = 1,ncol
         452           0 :           e(i,k) = minmod(eh(i,k),eh(i,k-1))
         453             :        end do
         454             :     end do
         455             :
         456             :
         457             :
         458           0 :     do k = 3,pver-1
         459             :
         460           0 :        do i = 1,ncol
         461             :
         462             : !           p prime at k-0.5
         463           0 :           ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)
         464             : !           p prime at k+0.5
         465           0 :           ppr(i,k)=sh(i,k)   - dh(i,k)  *delxh(i,k)
         466             :
         467           0 :           t = minmod(ppl(i,k),ppr(i,k))
         468             :
         469             : !           derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
         470           0 :           pp = sh(i,k-1) + d(i,k)*delxh(i,k-1)
         471             :
         472             : !           quartic estimate of fdot
         473             :           fdot(i,k) = pp                            &
         474             :                - delxh(i,k-1)*delxh(i,k)            &
         475           0 :                *(  eh(i,k-1)*(x(i,k+2)-x(i,k  ))    &
         476           0 :                + eh(i,k  )*(x(i,k  )-x(i,k-2))      &
         477           0 :                )/(x(i,k+2)-x(i,k-2))
         478             :
         479             : !           now limit it
         480             :           qpl = sh(i,k-1)       &
         481             :                + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
         482           0 :                d(i,k)  -e(i,k)*delxh(i,k))
         483             :           qpr = sh(i,k)         &
         484             :                + delxh(i,k  )*minmod(d(i,k)  +e(i,k)*delxh(i,k-1),        &
         485           0 :                d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))
         486             :
         487           0 :           fdot(i,k) = medan(fdot(i,k), qpl, qpr)
         488             :
         489           0 :           ttt = minmod(qpl, qpr)
         490           0 :           tmin = min(0._r8,3*s(i,k),1.5_r8*t,ttt)
         491           0 :           tmax = max(0._r8,3*s(i,k),1.5_r8*t,ttt)
         492             :
         493           0 :           fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))
         494             :
         495             :        end do
         496             :
         497             :     end do
         498             :
         499           0 :     return
         500             :   end subroutine cfdotmc_pro
         501             : end module dust_sediment_mod

Generated by: LCOV version 1.14