LCOV - code coverage report
Current view: | top level - physics/cam - pkg_cld_sediment.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 12 | 174 | 6.9 %
Date: | 2024-12-17 22:39:59 |  | Functions: | 1 | 6 | 16.7 %

              Line data    Source code
           1             : #undef OLDLIQSED
           2             : module pkg_cld_sediment
           3             :
           4             : !---------------------------------------------------------------------------------
           5             : ! Purpose:
           6             : !
           7             : ! Contains routines to compute tendencies from sedimentation of cloud liquid and
           8             : ! ice particles
           9             : !
          10             : ! Author: Byron Boville  Sept 19, 2002 from code by P. J. Rasch
          11             : !
          12             : !---------------------------------------------------------------------------------
          13             :
          14             :   use shr_kind_mod,   only: r8=>shr_kind_r8
          15             :   use spmd_utils,     only: masterproc
          16             :   use ppgrid,         only: pcols, pver, pverp
          17             :   use physconst,      only: gravit, latvap, latice, rair, rhoh2o
          18             :   use cldwat,         only: icritc
          19             :   use pkg_cldoptics,  only: reitab, reltab
          20             :   use phys_control,   only: use_simple_phys
          21             :   use cam_abortutils, only: endrun
          22             :   use cam_logfile,    only: iulog
          23             :
          24             :   implicit none
          25             :   private
          26             :   save
          27             :
          28             :   public :: cld_sediment_readnl, cld_sediment_vel, cld_sediment_tend
          29             :
          30             :
          31             :   real (r8), parameter :: vland  = 1.5_r8            ! liquid fall velocity over land  (cm/s)
          32             :   real (r8), parameter :: vocean = 2.8_r8            ! liquid fall velocity over ocean (cm/s)
          33             :   real (r8), parameter :: mxsedfac   = 0.99_r8       ! maximum sedimentation flux factor
          34             :
          35             :   logical,   parameter :: stokes = .true.         ! use Stokes velocity instead of McFarquhar and Heymsfield
          36             :
          37             : ! parameter for modified McFarquhar and Heymsfield
          38             :   real (r8), parameter :: vice_small = 1._r8         ! ice fall velocity for small concentration (cm/s)
          39             :
          40             : ! parameters for Stokes velocity
          41             :   real (r8), parameter :: eta =  1.7e-5_r8           ! viscosity of air (kg m / s)
          42             :   real (r8), parameter :: r40 =  40._r8              !  40 micron radius
          43             :   real (r8), parameter :: r400= 400._r8              ! 400 micron radius
          44             :   real (r8), parameter :: v400= 1.00_r8              ! fall velocity of 400 micron sphere (m/s)
          45             :   real (r8)            :: v40 ! = (2._r8/9._r8) * rhoh2o * gravit/eta * r40**2 * 1.e-12_r8
          46             :                                                      ! Stokes fall velocity of 40 micron sphere (m/s)
          47             :   real (r8)            :: vslope !  = (v400 - v40)/(r400 -r40) ! linear slope for large particles m/s/micron
          48             :
          49             :   ! namelist variables
          50             :   real(r8) :: cldsed_ice_stokes_fac = huge(1._r8)    ! factor applied to the ice fall velocity computed from
          51             :                                                      ! stokes terminal velocity
          52             :
          53             : !===============================================================================
          54             : contains
          55             : !===============================================================================
          56             :
          57        1536 : subroutine cld_sediment_readnl(nlfile)
          58             :
          59             :    use namelist_utils,  only: find_group_name
          60             :    use units,           only: getunit, freeunit
          61             :    use mpishorthand
          62             :
          63             :    character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          64             :
          65             :    ! Local variables
          66             :    integer :: unitn, ierr
          67             :    character(len=*), parameter :: subname = 'cld_sediment_readnl'
          68             :
          69             :    namelist /cldsed_nl/ cldsed_ice_stokes_fac
          70             :    !-----------------------------------------------------------------------------
          71             :
          72        1536 :    if (masterproc) then
          73           2 :       unitn = getunit()
          74           2 :       open( unitn, file=trim(nlfile), status='old' )
          75           2 :       call find_group_name(unitn, 'cldsed_nl', status=ierr)
          76           2 :       if (ierr == 0) then
          77           0 :          read(unitn, cldsed_nl, iostat=ierr)
          78           0 :          if (ierr /= 0) then
          79           0 :             call endrun(subname // ':: ERROR reading namelist')
          80             :          end if
          81             :       end if
          82           2 :       close(unitn)
          83           2 :       call freeunit(unitn)
          84             :
          85             :    end if
          86             :
          87             : #ifdef SPMD
          88             :    ! Broadcast namelist variables
          89        1536 :    call mpibcast(cldsed_ice_stokes_fac, 1, mpir8, 0, mpicom)
          90             : #endif
          91             :
          92        1536 :    if (masterproc .and. .not. use_simple_phys) then
          93           2 :       write(iulog,*) subname//': cldsed_ice_stokes_fac = ', cldsed_ice_stokes_fac
          94             :    end if
          95             :
          96        1536 : end subroutine cld_sediment_readnl
          97             :
          98             : !===============================================================================
          99             :
         100           0 :   subroutine cld_sediment_vel (ncol,                               &
         101             :        icefrac , landfrac, ocnfrac , pmid    , pdel    , t       , &
         102             :        cloud   , cldliq  , cldice  , pvliq   , pvice   , landm, snowh)
         103             :
         104             : !----------------------------------------------------------------------
         105             :
         106             : ! Compute gravitational sedimentation velocities for cloud liquid water
         107             : ! and ice, based on Lawrence and Crutzen (1998).
         108             :
         109             : ! LIQUID
         110             :
         111             : ! The fall velocities assume that droplets have a gamma distribution
         112             : ! with effective radii for land and ocean as assessed by Han et al.;
         113             : ! see Lawrence and Crutzen (1998) for a derivation.
         114             :
         115             : ! ICE
         116             :
         117             : ! The fall velocities are based on data from McFarquhar and Heymsfield
         118             : ! or on Stokes terminal velocity for spheres and the effective radius.
         119             :
         120             : ! NEED TO BE CAREFUL - VELOCITIES SHOULD BE AT THE *LOWER* INTERFACE
         121             : ! (THAT IS, FOR K+1), FLUXES ARE ALSO AT THE LOWER INTERFACE (K+1),
         122             : ! BUT MIXING RATIOS ARE AT THE MIDPOINTS (K)...
         123             :
         124             : ! NOTE THAT PVEL IS ON PVERP (INTERFACES), WHEREAS VFALL IS FOR THE CELL
         125             : ! AVERAGES (I.E., MIDPOINTS); ASSUME THE FALL VELOCITY APPLICABLE TO THE
         126             : ! LOWER INTERFACE (K+1) IS THE SAME AS THAT APPLICABLE FOR THE CELL (V(K))
         127             :
         128             : !-----------------------------------------------------------------------
         129             : !     MATCH-MPIC version 2.0, Author: mgl, March 1998
         130             : ! adapted by P. J. Rasch
         131             : !            B. A. Boville, September 19, 2002
         132             : !            P. J. Rasch    May 22, 2003 (added stokes flow calc for liquid
         133             : !                                         drops based on effect radii)
         134             : !-----------------------------------------------------------------------
         135             :
         136             :
         137             : ! Arguments
         138             :     integer, intent(in) :: ncol                     ! number of colums to process
         139             :
         140             :     real(r8), intent(in)  :: icefrac (pcols)        ! sea ice fraction (fraction)
         141             :     real(r8), intent(in)  :: landfrac(pcols)        ! land fraction (fraction)
         142             :     real(r8), intent(in)  :: ocnfrac (pcols)        ! ocean fraction (fraction)
         143             :     real(r8), intent(in)  :: pmid  (pcols,pver)     ! pressure of midpoint levels (Pa)
         144             :     real(r8), intent(in)  :: pdel  (pcols,pver)     ! pressure diff across layer (Pa)
         145             :     real(r8), intent(in)  :: cloud (pcols,pver)     ! cloud fraction (fraction)
         146             :     real(r8), intent(in)  :: t     (pcols,pver)     ! temperature (K)
         147             :     real(r8), intent(in)  :: cldliq(pcols,pver)     ! cloud water, liquid (kg/kg)
         148             :     real(r8), intent(in)  :: cldice(pcols,pver)     ! cloud water, ice    (kg/kg)
         149             :     real(r8), intent(in) :: snowh(pcols)         ! Snow depth over land, water equivalent (m)
         150             :
         151             :     real(r8), intent(out) :: pvliq (pcols,pverp)    ! vertical velocity of cloud liquid drops (Pa/s)
         152             :     real(r8), intent(out) :: pvice (pcols,pverp)    ! vertical velocity of cloud ice particles (Pa/s)
         153             :     real(r8), intent(in) :: landm(pcols)            ! land fraction ramped over water
         154             : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
         155             :
         156             : ! Local variables
         157             :     real (r8) :: rho(pcols,pver)                    ! air density in kg/m3
         158             :     real (r8) :: vfall                              ! settling velocity of cloud particles (m/s)
         159             :     real (r8) :: icice                              ! in cloud ice water content (kg/kg)
         160             :     real (r8) :: iciwc                              ! in cloud ice water content in g/m3
         161             :     real (r8) :: icefac
         162             :     real (r8) :: logiwc
         163             :
         164             :     real (r8) :: rei(pcols,pver)                    ! effective radius of ice particles (microns)
         165             :     real (r8) :: rel(pcols,pver)                    ! effective radius of liq particles (microns)
         166             :     real(r8)  pvliq2 (pcols,pverp)    ! vertical velocity of cloud liquid drops (Pa/s)
         167             :
         168             :     integer i,k
         169             :
         170             :     real (r8) :: lbound, ac, bc, cc
         171             :
         172             : !-----------------------------------------------------------------------
         173             : !------- initialize linear ramp variables for fall velocity ------------
         174             : !-----------------------------------------------------------------------
         175             :
         176           0 :   v40 = (2._r8/9._r8) * rhoh2o * gravit/eta * r40**2 * 1.e-12_r8
         177           0 :   vslope = (v400 - v40)/(r400 -r40)
         178             :
         179             : !-----------------------------------------------------------------------
         180             : !--------------------- liquid fall velocity ----------------------------
         181             : !-----------------------------------------------------------------------
         182             :
         183             :     ! compute air density
         184           0 :     rho(:ncol,:) = pmid(:ncol,:) / (rair * t(:ncol,:))
         185             :
         186           0 :     pvliq(:ncol,:) = 0._r8
         187             :
         188             :     ! get effective radius of liquid drop
         189           0 :     call reltab(ncol, t, landfrac, landm, icefrac, rel, snowh)
         190             :
         191           0 :     do k = 1,pver
         192           0 :        do i = 1,ncol
         193           0 :           if (cloud(i,k) > 0._r8 .and. cldliq(i,k) > 0._r8) then
         194             :
         195             : #ifdef OLDLIQSED
         196             : ! oldway
         197             :              ! merge the liquid fall velocities for land and ocean (cm/s)
         198             :              ! SHOULD ALSO ACCOUNT FOR ICEFRAC
         199             :              vfall = vland*landfrac(i) + vocean*(1._r8-landfrac(i))
         200             : !!$          vfall = vland*landfrac(i) + vocean*ocnfrac(i) + vseaice*icefrac(i)
         201             :
         202             :              ! convert the fall speed to pressure units, but do not apply the traditional
         203             :              ! negative convention for pvel.
         204             :              pvliq(i,k+1) = vfall     &
         205             :                   * 0.01_r8                 & ! cm to meters
         206             :                   * rho(i,k)*gravit        ! meters/sec to pascals/sec
         207             : #else
         208             :
         209             : ! newway
         210           0 :              if (rel(i,k) < 40._r8 ) then
         211           0 :                 vfall = 2._r8/9._r8 * rhoh2o * gravit * rel(i,k)**2 / eta  * 1.e-12_r8  ! micons^2 -> m^2
         212             :              else
         213           0 :                 vfall = v40 + vslope * (rel(i,k)-r40)      ! linear above 40 microns
         214             :              end if
         215             :              ! convert the fall speed to pressure units
         216             :              ! but do not apply the traditional
         217             :              ! negative convention for pvel.
         218             : !             pvliq2(i,k+1) = vfall * rho(i,k)*gravit        ! meters/sec to pascals/sec
         219           0 :              pvliq(i,k+1) = vfall * rho(i,k)*gravit        ! meters/sec to pascals/sec
         220             : #endif
         221             :           end if
         222             :        end do
         223             :     end do
         224             :
         225             :     !-----------------------------------------------------------------------
         226             :     !--------------------- ice fall velocity -------------------------------
         227             :     !-----------------------------------------------------------------------
         228             :
         229           0 :     pvice(:ncol,:) = 0._r8
         230             :
         231             :     if (stokes) then
         232             :
         233             :        !-----------------------------------------------------------------------
         234             :        !--------------------- stokes terminal velocity < 40 microns -----------
         235             :        !-----------------------------------------------------------------------
         236             :
         237             :        ! get effective radius
         238           0 :        call reitab(ncol, t, rei)
         239             :
         240           0 :        do k = 1,pver
         241           0 :           do i = 1,ncol
         242           0 :              if (cloud(i,k) > 0._r8 .and. cldice(i,k) > 0._r8) then
         243           0 :                 if (rei(i,k) < 40._r8 ) then
         244           0 :                    vfall = 2._r8/9._r8 * rhoh2o * gravit * rei(i,k)**2 / eta  * 1.e-12_r8  ! micons^2 -> m^2
         245           0 :                    vfall = vfall * cldsed_ice_stokes_fac
         246             :                 else
         247           0 :                    vfall = v40 + vslope * (rei(i,k)-r40)      ! linear above 40 microns
         248             :                 end if
         249             :
         250             :                 ! convert the fall speed to pressure units, but do not apply the traditional
         251             :                 ! negative convention for pvel.
         252           0 :                 pvice(i,k+1) = vfall * rho(i,k)*gravit        ! meters/sec to pascals/sec
         253             :              end if
         254             :           end do
         255             :        end do
         256             :
         257             :     else
         258             :
         259             :        !-----------------------------------------------------------------------
         260             :        !--------------------- McFarquhar and Heymsfield > icritc --------------
         261             :        !-----------------------------------------------------------------------
         262             :
         263             :        ! lower bound for iciwc
         264             :
         265             :        cc = 128.64_r8
         266             :        bc = 53.242_r8
         267             :        ac = 5.4795_r8
         268             :        lbound = (-bc + sqrt(bc*bc-4*ac*cc))/(2*ac)
         269             :        lbound = 10._r8**lbound
         270             :
         271             :        do k = 1,pver
         272             :           do i = 1,ncol
         273             :              if (cloud(i,k) > 0._r8 .and. cldice(i,k) > 0._r8) then
         274             :
         275             :                 ! compute the in-cloud ice concentration (kg/kg)
         276             :                 icice = cldice(i,k) / cloud(i,k)
         277             :
         278             :                 ! compute the ice water content in g/m3
         279             :                 iciwc = icice * rho(i,k) * 1.e3_r8
         280             :
         281             :                 ! set the fall velocity (cm/s) to depend on the ice water content in g/m3,
         282             :                 if (iciwc > lbound) then ! need this because of log10
         283             :                    logiwc = log10(iciwc)
         284             :                    !          Median -
         285             :                    vfall = 128.64_r8 + 53.242_r8*logiwc + 5.4795_r8*logiwc**2
         286             :                    !          Average -
         287             :                 !!$             vfall = 122.63 + 44.111*logiwc + 4.2144*logiwc**2
         288             :                 else
         289             :                    vfall = 0._r8
         290             :                 end if
         291             :
         292             :                 ! set ice velocity to 1 cm/s if ice mixing ratio < icritc, ramp to value
         293             :                 ! calculated above at 2*icritc
         294             :                 if (icice <= icritc) then
         295             :                    vfall = vice_small
         296             :                 else if(icice < 2*icritc) then
         297             :                    icefac = (icice-icritc)/icritc
         298             :                    vfall = vice_small * (1._r8-icefac) + vfall * icefac
         299             :                 end if
         300             :
         301             :                 ! bound the terminal velocity of ice particles at high concentration
         302             :                 vfall = min(100.0_r8, vfall)
         303             :
         304             :                 ! convert the fall speed to pressure units, but do not apply the traditional
         305             :                 ! negative convention for pvel.
         306             :                 pvice(i,k+1) = vfall     &
         307             :                    * 0.01_r8                 & ! cm to meters
         308             :                    * rho(i,k)*gravit        ! meters/sec to pascals/sec
         309             :              end if
         310             :           end do
         311             :        end do
         312             :
         313             :     end if
         314             :
         315           0 :  end subroutine cld_sediment_vel
         316             :
         317             :
         318             : !===============================================================================
         319           0 :   subroutine cld_sediment_tend (ncol, dtime  ,               &
         320             :        pint   , pmid   , pdel   , t      ,                   &
         321             :        cloud  , cldliq , cldice , pvliq  , pvice  ,          &
         322             :        liqtend, icetend, wvtend , htend  , sfliq  , sfice   )
         323             :
         324             : !----------------------------------------------------------------------
         325             : !     Apply Cloud Particle Gravitational Sedimentation to Condensate
         326             : !----------------------------------------------------------------------
         327             :
         328             :
         329             : ! Arguments
         330             :     integer,  intent(in)  :: ncol                      ! number of colums to process
         331             :
         332             :     real(r8), intent(in)  :: dtime                     ! time step
         333             :     real(r8), intent(in)  :: pint  (pcols,pverp)       ! interfaces pressure (Pa)
         334             :     real(r8), intent(in)  :: pmid  (pcols,pver)        ! midpoint pressures (Pa)
         335             :     real(r8), intent(in)  :: pdel  (pcols,pver)        ! pressure diff across layer (Pa)
         336             :     real(r8), intent(in)  :: cloud (pcols,pver)        ! cloud fraction (fraction)
         337             :     real(r8), intent(in)  :: t     (pcols,pver)        ! temperature (K)
         338             :     real(r8), intent(in)  :: cldliq(pcols,pver)        ! cloud liquid water (kg/kg)
         339             :     real(r8), intent(in)  :: cldice(pcols,pver)        ! cloud ice water    (kg/kg)
         340             :     real(r8), intent(in)  :: pvliq (pcols,pverp)       ! vertical velocity of liquid drops  (Pa/s)
         341             :     real(r8), intent(in)  :: pvice (pcols,pverp)       ! vertical velocity of ice particles (Pa/s)
         342             : ! -> note that pvel is at the interfaces (loss from cell is based on pvel(k+1))
         343             :
         344             :     real(r8), intent(out) :: liqtend(pcols,pver)       ! liquid condensate tend
         345             :     real(r8), intent(out) :: icetend(pcols,pver)       ! ice condensate tend
         346             :     real(r8), intent(out) :: wvtend (pcols,pver)       ! water vapor tend
         347             :     real(r8), intent(out) :: htend  (pcols,pver)       ! heating rate
         348             :     real(r8), intent(out) :: sfliq  (pcols)            ! surface flux of liquid (rain, kg/m/s)
         349             :     real(r8), intent(out) :: sfice  (pcols)            ! surface flux of ice    (snow, kg/m/s)
         350             :
         351             : ! Local variables
         352             :     real(r8) :: fxliq(pcols,pverp)                     ! fluxes at the interfaces, liquid (positive = down)
         353             :     real(r8) :: fxice(pcols,pverp)                     ! fluxes at the interfaces, ice    (positive = down)
         354             :     real(r8) :: cldab(pcols)                           ! cloud in layer above
         355             :     real(r8) :: evapliq                                ! evaporation of cloud liquid into environment
         356             :     real(r8) :: evapice                                ! evaporation of cloud ice into environment
         357             :     real(r8) :: cldovrl                                ! cloud overlap factor
         358             :
         359             :     integer :: i,k
         360             : !----------------------------------------------------------------------
         361             :
         362             : ! initialize variables
         363           0 :     fxliq  (:ncol,:) = 0._r8 ! flux at interfaces (liquid)
         364           0 :     fxice  (:ncol,:) = 0._r8 ! flux at interfaces (ice)
         365           0 :     liqtend(:ncol,:) = 0._r8 ! condensate tend (liquid)
         366           0 :     icetend(:ncol,:) = 0._r8 ! condensate tend (ice)
         367           0 :     wvtend(:ncol,:)  = 0._r8 ! environmental moistening
         368           0 :     htend(:ncol,:)   = 0._r8 ! evaporative cooling
         369           0 :     sfliq(:ncol)     = 0._r8 ! condensate sedimentation flux out bot of column (liquid)
         370           0 :     sfice(:ncol)     = 0._r8 ! condensate sedimentation flux out bot of column (ice)
         371             :
         372             : ! fluxes at interior points
         373           0 :     call getflx(ncol, pint, cldliq, pvliq, dtime, fxliq)
         374           0 :     call getflx(ncol, pint, cldice, pvice, dtime, fxice)
         375             :
         376             : ! calculate fluxes at boundaries
         377           0 :     do i = 1,ncol
         378           0 :        fxliq(i,1) = 0._r8
         379           0 :        fxice(i,1) = 0._r8
         380             : ! surface flux by upstream scheme
         381           0 :        fxliq(i,pverp) = cldliq(i,pver) * pvliq(i,pverp) * dtime
         382           0 :        fxice(i,pverp) = cldice(i,pver) * pvice(i,pverp) * dtime
         383             :     end do
         384             :
         385             : ! filter out any negative fluxes from the getflx routine
         386             : ! (typical fluxes are of order > 1e-3 when clouds are present)
         387           0 :     do k = 2,pver
         388           0 :        fxliq(:ncol,k) = max(0._r8, fxliq(:ncol,k))
         389           0 :        fxice(:ncol,k) = max(0._r8, fxice(:ncol,k))
         390             :     end do
         391             :
         392             : ! Limit the flux out of the bottom of each cell to the water content in each phase.
         393             : ! Apply mxsedfac to prevent generating very small negative cloud water/ice
         394             : ! NOTE, REMOVED CLOUD FACTOR FROM AVAILABLE WATER. ALL CLOUD WATER IS IN CLOUDS.
         395             : ! ***Should we include the flux in the top, to allow for thin surface layers?
         396             : ! ***Requires simple treatment of cloud overlap, already included below.
         397           0 :     do k = 1,pver
         398           0 :        do i = 1,ncol
         399           0 :           fxliq(i,k+1) = min( fxliq(i,k+1), mxsedfac * cldliq(i,k) * pdel(i,k) )
         400           0 :           fxice(i,k+1) = min( fxice(i,k+1), mxsedfac * cldice(i,k) * pdel(i,k) )
         401             : !!$        fxliq(i,k+1) = min( fxliq(i,k+1), cldliq(i,k) * pdel(i,k) + fxliq(i,k))
         402             : !!$        fxice(i,k+1) = min( fxice(i,k+1), cldice(i,k) * pdel(i,k) + fxice(i,k))
         403             : !!$        fxliq(i,k+1) = min( fxliq(i,k+1), cloud(i,k) * cldliq(i,k) * pdel(i,k) )
         404             : !!$        fxice(i,k+1) = min( fxice(i,k+1), cloud(i,k) * cldice(i,k) * pdel(i,k) )
         405             :        end do
         406             :     end do
         407             :
         408             : ! Now calculate the tendencies assuming that condensate evaporates when
         409             : ! it falls into environment, and does not when it falls into cloud.
         410             : ! All flux out of the layer comes from the cloudy part.
         411             : ! Assume maximum overlap for stratiform clouds
         412             : !  if cloud above < cloud,  all water falls into cloud below
         413             : !  if cloud above > cloud,  water split between cloud  and environment
         414           0 :     do k = 1,pver
         415           0 :        cldab(:ncol) = 0._r8
         416           0 :        do i = 1,ncol
         417             : ! cloud overlap cloud factor
         418           0 :           cldovrl  = min( cloud(i,k) / (cldab(i)+.0001_r8), 1._r8 )
         419           0 :           cldab(i) = cloud(i,k)
         420             : ! evaporation into environment cause moistening and cooling
         421           0 :           evapliq = fxliq(i,k) * (1._r8-cldovrl) / (dtime * pdel(i,k))  ! into env (kg/kg/s)
         422           0 :           evapice = fxice(i,k) * (1._r8-cldovrl) / (dtime * pdel(i,k))  ! into env (kg/kg/s)
         423           0 :           wvtend(i,k) = evapliq + evapice                          ! evaporation into environment (kg/kg/s)
         424           0 :           htend (i,k) = -latvap*evapliq -(latvap+latice)*evapice   ! evaporation (W/kg)
         425             : ! net flux into cloud changes cloud liquid/ice (all flux is out of cloud)
         426           0 :           liqtend(i,k)  = (fxliq(i,k)*cldovrl - fxliq(i,k+1)) / (dtime * pdel(i,k))
         427           0 :           icetend(i,k)  = (fxice(i,k)*cldovrl - fxice(i,k+1)) / (dtime * pdel(i,k))
         428             :        end do
         429             :     end do
         430             :
         431             : ! convert flux out the bottom to mass units Pa -> kg/m2/s
         432           0 :     sfliq(:ncol) = fxliq(:ncol,pverp) / (dtime*gravit)
         433           0 :     sfice(:ncol) = fxice(:ncol,pverp) / (dtime*gravit)
         434             :
         435           0 :     return
         436             :   end subroutine cld_sediment_tend
         437             :
         438             : !===============================================================================
         439           0 :   subroutine getflx(ncol, xw, phi, vel, deltat, flux)
         440             :
         441             : !.....xw1.......xw2.......xw3.......xw4.......xw5.......xw6
         442             : !....psiw1.....psiw2.....psiw3.....psiw4.....psiw5.....psiw6
         443             : !....velw1.....velw2.....velw3.....velw4.....velw5.....velw6
         444             : !.........phi1......phi2.......phi3.....phi4.......phi5.......
         445             :
         446             :
         447             :
         448             :     integer, intent(in) :: ncol                      ! number of colums to process
         449             :
         450             :     integer i
         451             :     integer k
         452             :
         453             :     real (r8), intent(in) :: vel(pcols,pverp)
         454             :     real (r8) flux(pcols,pverp)
         455             :     real (r8) xw(pcols,pverp)
         456             :     real (r8) psi(pcols,pverp)
         457             :     real (r8), intent(in) :: phi(pcols,pverp-1)
         458             :     real (r8) fdot(pcols,pverp)
         459             :     real (r8) xx(pcols)
         460             :     real (r8) fxdot(pcols)
         461             :     real (r8) fxdd(pcols)
         462             :
         463             :     real (r8) psistar(pcols)
         464             :     real (r8) deltat
         465             :
         466             :     real (r8) xxk(pcols,pver)
         467             :
         468           0 :     do i = 1,ncol
         469             : !        integral of phi
         470           0 :        psi(i,1) = 0._r8
         471             : !        fluxes at boundaries
         472           0 :        flux(i,1) = 0._r8
         473           0 :        flux(i,pverp) = 0._r8
         474             :     end do
         475             :
         476             : !     integral function
         477           0 :     do k = 2,pverp
         478           0 :        do i = 1,ncol
         479           0 :           psi(i,k) = phi(i,k-1)*(xw(i,k)-xw(i,k-1)) + psi(i,k-1)
         480             :        end do
         481             :     end do
         482             :
         483             :
         484             : !     calculate the derivatives for the interpolating polynomial
         485           0 :     call cfdotmc_pro (ncol, xw, psi, fdot)
         486             :
         487             : !  NEW WAY
         488             : !     calculate fluxes at interior pts
         489           0 :     do k = 2,pver
         490           0 :        do i = 1,ncol
         491           0 :           xxk(i,k) = xw(i,k)-vel(i,k)*deltat
         492             :        end do
         493             :     end do
         494           0 :     do k = 2,pver
         495           0 :        call cfint2(ncol, xw, psi, fdot, xxk(1,k), fxdot, fxdd, psistar)
         496           0 :        do i = 1,ncol
         497           0 :           flux(i,k) = (psi(i,k)-psistar(i))
         498             :        end do
         499             :     end do
         500             :
         501             :
         502           0 :     return
         503             :   end subroutine getflx
         504             :
         505             :
         506             :
         507             : !##############################################################################
         508             :
         509           0 :   subroutine cfint2 (ncol, x, f, fdot, xin, fxdot, fxdd, psistar)
         510             :
         511             :
         512             :
         513             : ! input
         514             :     integer ncol                      ! number of colums to process
         515             :
         516             :     real (r8) x(pcols, pverp)
         517             :     real (r8) f(pcols, pverp)
         518             :     real (r8) fdot(pcols, pverp)
         519             :     real (r8) xin(pcols)
         520             :
         521             : ! output
         522             :     real (r8) fxdot(pcols)
         523             :     real (r8) fxdd(pcols)
         524             :     real (r8) psistar(pcols)
         525             :
         526             :     integer i
         527             :     integer k
         528             :     integer intz(pcols)
         529             :     real (r8) dx
         530             :     real (r8) s
         531             :     real (r8) c2
         532             :     real (r8) c3
         533             :     real (r8) xx
         534             :     real (r8) xinf
         535             :     real (r8) psi1, psi2, psi3, psim
         536             :     real (r8) cfint
         537             :     real (r8) cfnew
         538             :     real (r8) xins(pcols)
         539             :
         540             : !     the minmod function
         541             :     real (r8) a, b, c
         542             :     real (r8) minmod
         543             :     real (r8) medan
         544             :     logical found_error
         545             :
         546             :     minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
         547             :     medan(a,b,c) = a + minmod(b-a,c-a)
         548             :
         549           0 :     do i = 1,ncol
         550           0 :        xins(i) = medan(x(i,1), xin(i), x(i,pverp))
         551           0 :        intz(i) = 0
         552             :     end do
         553             :
         554             : ! first find the interval
         555           0 :     do k =  1,pverp-1
         556           0 :        do i = 1,ncol
         557           0 :           if ((xins(i)-x(i,k))*(x(i,k+1)-xins(i)).ge.0) then
         558           0 :              intz(i) = k
         559             :           endif
         560             :        end do
         561             :     end do
         562             :
         563           0 :     found_error=.false.
         564           0 :     do i = 1,ncol
         565           0 :        if (intz(i).eq.0._r8) found_error=.true.
         566             :     end do
         567           0 :     if(found_error) then
         568           0 :        do i = 1,ncol
         569           0 :           if (intz(i).eq.0._r8) then
         570           0 :              write(iulog,*) ' interval was not found for col i ', i
         571           0 :              call endrun('CFINT2')
         572             :           endif
         573             :        end do
         574             :     endif
         575             :
         576             : ! now interpolate
         577           0 :     do i = 1,ncol
         578           0 :        k = intz(i)
         579           0 :        dx = (x(i,k+1)-x(i,k))
         580           0 :        s = (f(i,k+1)-f(i,k))/dx
         581           0 :        c2 = (3*s-2*fdot(i,k)-fdot(i,k+1))/dx
         582           0 :        c3 = (fdot(i,k)+fdot(i,k+1)-2*s)/dx**2
         583           0 :        xx = (xins(i)-x(i,k))
         584           0 :        fxdot(i) =  (3*c3*xx + 2*c2)*xx + fdot(i,k)
         585           0 :        fxdd(i) = 6*c3*xx + 2*c2
         586           0 :        cfint = ((c3*xx + c2)*xx + fdot(i,k))*xx + f(i,k)
         587             :
         588             : !        limit the interpolant
         589           0 :        psi1 = f(i,k)+(f(i,k+1)-f(i,k))*xx/dx
         590           0 :        if (k.eq.1) then
         591           0 :           psi2 = f(i,1)
         592             :        else
         593           0 :           psi2 = f(i,k) + (f(i,k)-f(i,k-1))*xx/(x(i,k)-x(i,k-1))
         594             :        endif
         595           0 :        if (k+1.eq.pverp) then
         596           0 :           psi3 = f(i,pverp)
         597             :        else
         598           0 :           psi3 = f(i,k+1) - (f(i,k+2)-f(i,k+1))*(dx-xx)/(x(i,k+2)-x(i,k+1))
         599             :        endif
         600           0 :        psim = medan(psi1, psi2, psi3)
         601           0 :        cfnew = medan(cfint, psi1, psim)
         602             :        if (abs(cfnew-cfint)/(abs(cfnew)+abs(cfint)+1.e-36_r8)  .gt..03_r8) then
         603             : !     CHANGE THIS BACK LATER!!!
         604             : !     $        .gt..1) then
         605             :
         606             :
         607             : !     UNCOMMENT THIS LATER!!!
         608             : !            write(iulog,*) ' cfint2 limiting important ', cfint, cfnew
         609             :
         610             :
         611             :        endif
         612           0 :        psistar(i) = cfnew
         613             :     end do
         614             :
         615           0 :     return
         616             :   end subroutine cfint2
         617             :
         618             :
         619             :
         620             : !##############################################################################
         621             :
         622           0 :   subroutine cfdotmc_pro (ncol, x, f, fdot)
         623             :
         624             : !     prototype version; eventually replace with final SPITFIRE scheme
         625             :
         626             : !     calculate the derivative for the interpolating polynomial
         627             : !     multi column version
         628             :
         629             :
         630             :
         631             : ! input
         632             :     integer ncol                      ! number of colums to process
         633             :
         634             :     real (r8) x(pcols, pverp)
         635             :     real (r8) f(pcols, pverp)
         636             : ! output
         637             :     real (r8) fdot(pcols, pverp)          ! derivative at nodes
         638             :
         639             : ! assumed variable distribution
         640             : !     x1.......x2.......x3.......x4.......x5.......x6     1,pverp points
         641             : !     f1.......f2.......f3.......f4.......f5.......f6     1,pverp points
         642             : !     ...sh1.......sh2......sh3......sh4......sh5....     1,pver points
         643             : !     .........d2.......d3.......d4.......d5.........     2,pver points
         644             : !     .........s2.......s3.......s4.......s5.........     2,pver points
         645             : !     .............dh2......dh3......dh4.............     2,pver-1 points
         646             : !     .............eh2......eh3......eh4.............     2,pver-1 points
         647             : !     ..................e3.......e4..................     3,pver-1 points
         648             : !     .................ppl3......ppl4................     3,pver-1 points
         649             : !     .................ppr3......ppr4................     3,pver-1 points
         650             : !     .................t3........t4..................     3,pver-1 points
         651             : !     ................fdot3.....fdot4................     3,pver-1 points
         652             :
         653             :
         654             : ! work variables
         655             :
         656             :
         657             :     integer i
         658             :     integer k
         659             :
         660             :     real (r8) a                    ! work var
         661             :     real (r8) b                    ! work var
         662             :     real (r8) c                    ! work var
         663             :     real (r8) s(pcols,pverp)             ! first divided differences at nodes
         664             :     real (r8) sh(pcols,pverp)            ! first divided differences between nodes
         665             :     real (r8) d(pcols,pverp)             ! second divided differences at nodes
         666             :     real (r8) dh(pcols,pverp)            ! second divided differences between nodes
         667             :     real (r8) e(pcols,pverp)             ! third divided differences at nodes
         668             :     real (r8) eh(pcols,pverp)            ! third divided differences between nodes
         669             :     real (r8) pp                   ! p prime
         670             :     real (r8) ppl(pcols,pverp)           ! p prime on left
         671             :     real (r8) ppr(pcols,pverp)           ! p prime on right
         672             :     real (r8) qpl
         673             :     real (r8) qpr
         674             :     real (r8) ttt
         675             :     real (r8) t
         676             :     real (r8) tmin
         677             :     real (r8) tmax
         678             :     real (r8) delxh(pcols,pverp)
         679             :
         680             :
         681             : !     the minmod function
         682             :     real (r8) minmod
         683             :     real (r8) medan
         684             :     minmod(a,b) = 0.5_r8*(sign(1._r8,a) + sign(1._r8,b))*min(abs(a),abs(b))
         685             :     medan(a,b,c) = a + minmod(b-a,c-a)
         686             :
         687           0 :     do k = 1,pver
         688             :
         689             :
         690             : !        first divided differences between nodes
         691           0 :        do i = 1, ncol
         692           0 :           delxh(i,k) = (x(i,k+1)-x(i,k))
         693           0 :           sh(i,k) = (f(i,k+1)-f(i,k))/delxh(i,k)
         694             :        end do
         695             :
         696             : !        first and second divided differences at nodes
         697           0 :        if (k.ge.2) then
         698           0 :           do i = 1,ncol
         699           0 :              d(i,k) = (sh(i,k)-sh(i,k-1))/(x(i,k+1)-x(i,k-1))
         700           0 :              s(i,k) = minmod(sh(i,k),sh(i,k-1))
         701             :           end do
         702             :        endif
         703             :     end do
         704             :
         705             : !     second and third divided diffs between nodes
         706           0 :     do k = 2,pver-1
         707           0 :        do i = 1, ncol
         708           0 :           eh(i,k) = (d(i,k+1)-d(i,k))/(x(i,k+2)-x(i,k-1))
         709           0 :           dh(i,k) = minmod(d(i,k),d(i,k+1))
         710             :        end do
         711             :     end do
         712             :
         713             : !     treat the boundaries
         714           0 :     do i = 1,ncol
         715           0 :        e(i,2) = eh(i,2)
         716           0 :        e(i,pver) = eh(i,pver-1)
         717             : !        outside level
         718             :        fdot(i,1) = sh(i,1) - d(i,2)*delxh(i,1)  &
         719           0 :             - eh(i,2)*delxh(i,1)*(x(i,1)-x(i,3))
         720           0 :        fdot(i,1) = minmod(fdot(i,1),3*sh(i,1))
         721             :        fdot(i,pverp) = sh(i,pver) + d(i,pver)*delxh(i,pver)  &
         722           0 :             + eh(i,pver-1)*delxh(i,pver)*(x(i,pverp)-x(i,pver-1))
         723           0 :        fdot(i,pverp) = minmod(fdot(i,pverp),3*sh(i,pver))
         724             : !        one in from boundary
         725           0 :        fdot(i,2) = sh(i,1) + d(i,2)*delxh(i,1) - eh(i,2)*delxh(i,1)*delxh(i,2)
         726           0 :        fdot(i,2) = minmod(fdot(i,2),3*s(i,2))
         727             :        fdot(i,pver) = sh(i,pver) - d(i,pver)*delxh(i,pver)   &
         728           0 :             - eh(i,pver-1)*delxh(i,pver)*delxh(i,pver-1)
         729           0 :        fdot(i,pver) = minmod(fdot(i,pver),3*s(i,pver))
         730             :     end do
         731             :
         732             :
         733           0 :     do k = 3,pver-1
         734           0 :        do i = 1,ncol
         735           0 :           e(i,k) = minmod(eh(i,k),eh(i,k-1))
         736             :        end do
         737             :     end do
         738             :
         739             :
         740             :
         741           0 :     do k = 3,pver-1
         742             :
         743           0 :        do i = 1,ncol
         744             :
         745             : !           p prime at k-0.5
         746           0 :           ppl(i,k)=sh(i,k-1) + dh(i,k-1)*delxh(i,k-1)
         747             : !           p prime at k+0.5
         748           0 :           ppr(i,k)=sh(i,k)   - dh(i,k)  *delxh(i,k)
         749             :
         750           0 :           t = minmod(ppl(i,k),ppr(i,k))
         751             :
         752             : !           derivate from parabola thru f(i,k-1), f(i,k), and f(i,k+1)
         753           0 :           pp = sh(i,k-1) + d(i,k)*delxh(i,k-1)
         754             :
         755             : !           quartic estimate of fdot
         756             :           fdot(i,k) = pp                            &
         757             :                - delxh(i,k-1)*delxh(i,k)            &
         758           0 :                *(  eh(i,k-1)*(x(i,k+2)-x(i,k  ))    &
         759           0 :                + eh(i,k  )*(x(i,k  )-x(i,k-2))      &
         760           0 :                )/(x(i,k+2)-x(i,k-2))
         761             :
         762             : !           now limit it
         763             :           qpl = sh(i,k-1)       &
         764             :                + delxh(i,k-1)*minmod(d(i,k-1)+e(i,k-1)*(x(i,k)-x(i,k-2)), &
         765           0 :                d(i,k)  -e(i,k)*delxh(i,k))
         766             :           qpr = sh(i,k)         &
         767             :                + delxh(i,k  )*minmod(d(i,k)  +e(i,k)*delxh(i,k-1),        &
         768           0 :                d(i,k+1)+e(i,k+1)*(x(i,k)-x(i,k+2)))
         769             :
         770           0 :           fdot(i,k) = medan(fdot(i,k), qpl, qpr)
         771             :
         772           0 :           ttt = minmod(qpl, qpr)
         773           0 :           tmin = min(0._r8,3*s(i,k),1.5_r8*t,ttt)
         774           0 :           tmax = max(0._r8,3*s(i,k),1.5_r8*t,ttt)
         775             :
         776           0 :           fdot(i,k) = fdot(i,k) + minmod(tmin-fdot(i,k), tmax-fdot(i,k))
         777             :
         778             :        end do
         779             :
         780             :     end do
         781             :
         782           0 :     return
         783             :   end subroutine cfdotmc_pro
         784             : end module pkg_cld_sediment

Generated by: LCOV version 1.14