LCOV - code coverage report
Current view: | top level - physics/cam - ndrop_bam.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 181 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 4 | 0.0 %

              Line data    Source code
           1             : module ndrop_bam
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : !
           5             : !   CAM Interface for droplet activation by bulk aerosols.
           6             : !
           7             : !---------------------------------------------------------------------------------
           8             :
           9             : use shr_kind_mod,     only: r8=>shr_kind_r8
          10             : use spmd_utils,       only: masterproc
          11             : use ppgrid,           only: pcols, pver, pverp
          12             : use physconst,        only: gravit, rair, tmelt, cpair, rh2o, &
          13             :      r_universal, mwh2o, rhoh2o, latvap
          14             :
          15             : use rad_constituents, only: rad_cnst_get_info, rad_cnst_get_aer_props
          16             :
          17             : use shr_spfn_mod,     only: erf => shr_spfn_erf, &
          18             :                             erfc => shr_spfn_erfc
          19             : use wv_saturation,    only: qsat
          20             : use cam_history,      only: addfld, add_default, outfld
          21             : use cam_logfile,      only: iulog
          22             : use cam_abortutils,   only: endrun
          23             : use ref_pres,         only: top_lev=>trop_cloud_top_lev
          24             :
          25             : implicit none
          26             : private
          27             : save
          28             :
          29             : public :: ndrop_bam_init, ndrop_bam_run, ndrop_bam_ccn
          30             :
          31             : ! activate parameters
          32             :
          33             : real(r8) :: pi     ! pi
          34             : integer,  parameter :: psat = 6           ! number of supersaturations to calc ccn concentration
          35             : real(r8), parameter :: supersat(psat) = & ! supersaturation (%) to determine ccn concentration
          36             :                        (/0.02_r8,0.05_r8,0.1_r8,0.2_r8,0.5_r8,1.0_r8/)
          37             : real(r8)              :: super(psat)
          38             : real(r8), allocatable :: ccnfact(:,:)
          39             :
          40             : real(r8), allocatable :: alogsig(:)       ! natl log of geometric standard dev of aerosol
          41             : real(r8), allocatable :: exp45logsig(:)
          42             : real(r8), allocatable :: argfactor(:)
          43             : real(r8), allocatable :: amcube(:)        ! cube of dry mode radius (m)
          44             : real(r8), allocatable :: smcrit(:)        ! critical supersatuation for activation
          45             : real(r8), allocatable :: lnsm(:)          ! ln(smcrit)
          46             : real(r8), allocatable :: amcubefactor(:)  ! factors for calculating mode radius
          47             : real(r8), allocatable :: smcritfactor(:)  ! factors for calculating critical supersaturation
          48             : real(r8), allocatable :: f1(:), f2(:)     ! abdul-razzak functions of width
          49             :
          50             : real(r8) :: aten
          51             : real(r8) :: third, sixth
          52             : real(r8) :: sq2, sqpi
          53             : real(r8) :: alogten, alog2, alog3, alogaten
          54             :
          55             : real(r8) :: pref = 1013.25e2_r8 ! reference pressure (Pa)
          56             :
          57             : ! aerosol properties
          58             : character(len=20), allocatable :: aername(:)
          59             : real(r8), allocatable :: dryrad_aer(:)
          60             : real(r8), allocatable :: density_aer(:)
          61             : real(r8), allocatable :: hygro_aer(:)
          62             : real(r8), allocatable :: dispersion_aer(:)
          63             : real(r8), allocatable :: num_to_mass_aer(:)
          64             :
          65             : integer :: naer_all      ! number of aerosols affecting climate
          66             : integer :: idxsul   = -1 ! index in aerosol list for sulfate
          67             :
          68             : !===============================================================================
          69             : contains
          70             : !===============================================================================
          71             :
          72           0 : subroutine ndrop_bam_init
          73             :
          74             :   use phys_control, only: phys_getopts
          75             :
          76             :    !-----------------------------------------------------------------------
          77             :    !
          78             :    ! Initialize constants for droplet activation by bulk aerosols
          79             :    !
          80             :    !-----------------------------------------------------------------------
          81             :
          82             :    integer  :: l, m, iaer
          83             :    real(r8) :: surften       ! surface tension of water w/respect to air (N/m)
          84             :    real(r8) :: arg
          85             :    logical  :: history_amwg
          86             :    !-------------------------------------------------------------------------------
          87             :
          88           0 :    call phys_getopts(history_amwg_out=history_amwg)
          89             :
          90             :    ! Access the physical properties of the bulk aerosols that are affecting the climate
          91             :    ! by using routines from the rad_constituents module.
          92             :
          93           0 :    call rad_cnst_get_info(0, naero=naer_all)
          94             :    allocate( &
          95           0 :       aername(naer_all),        &
          96           0 :       dryrad_aer(naer_all),     &
          97           0 :       density_aer(naer_all),    &
          98           0 :       hygro_aer(naer_all),      &
          99           0 :       dispersion_aer(naer_all), &
         100           0 :       num_to_mass_aer(naer_all) )
         101             :
         102           0 :    do iaer = 1, naer_all
         103             :       call rad_cnst_get_aer_props(0, iaer, &
         104           0 :          aername         = aername(iaer), &
         105           0 :          dryrad_aer      = dryrad_aer(iaer), &
         106           0 :          density_aer     = density_aer(iaer), &
         107           0 :          hygro_aer       = hygro_aer(iaer), &
         108           0 :          dispersion_aer  = dispersion_aer(iaer), &
         109           0 :          num_to_mass_aer = num_to_mass_aer(iaer) )
         110             :
         111             :       ! Look for sulfate aerosol in this list (Bulk aerosol only)
         112           0 :       if (trim(aername(iaer)) == 'SULFATE') idxsul   = iaer
         113             :
         114             :       ! aerosol number concentration
         115           0 :       call addfld(trim(aername(iaer))//'_m3', (/ 'lev' /), 'A', 'm-3', 'aerosol number concentration')
         116             :
         117             :    end do
         118             :
         119           0 :    if (masterproc) then
         120           0 :       write(iulog,*) 'ndrop_bam_init: iaer, name, dryrad, density, hygro, dispersion, num_to_mass'
         121           0 :       do iaer = 1, naer_all
         122           0 :          write(iulog,*) iaer, aername(iaer), dryrad_aer(iaer), density_aer(iaer), hygro_aer(iaer), &
         123           0 :             dispersion_aer(iaer), num_to_mass_aer(iaer)
         124             :       end do
         125           0 :       if (idxsul < 1) then
         126           0 :          write(iulog,*) 'ndrop_bam_init: SULFATE aerosol properties NOT FOUND'
         127             :       else
         128           0 :          write(iulog,*) 'ndrop_bam_init: SULFATE aerosol properties FOUND at index ',idxsul
         129             :       end if
         130             :    end if
         131             :
         132           0 :    call addfld ('CCN1',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.02%')
         133           0 :    call addfld ('CCN2',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.05%')
         134           0 :    call addfld ('CCN3',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.1%')
         135           0 :    call addfld ('CCN4',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.2%')
         136           0 :    call addfld ('CCN5',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=0.5%')
         137           0 :    call addfld ('CCN6',(/ 'lev' /), 'A','#/cm3','CCN concentration at S=1.0%')
         138             :
         139           0 :    if (history_amwg) then
         140           0 :       call add_default('CCN3', 1, ' ')
         141             :    end if
         142             :
         143             :    ! set parameters for droplet activation, following abdul-razzak and ghan 2000, JGR
         144             :
         145           0 :    third   = 1._r8/3._r8
         146           0 :    sixth   = 1._r8/6._r8
         147           0 :    sq2     = sqrt(2._r8)
         148           0 :    pi      = 4._r8*atan(1.0_r8)
         149           0 :    sqpi    = sqrt(pi)
         150           0 :    surften = 0.076_r8
         151           0 :    aten    = 2._r8*mwh2o*surften/(r_universal*tmelt*rhoh2o)
         152           0 :    alogaten= log(aten)
         153           0 :    alog2   = log(2._r8)
         154           0 :    alog3   = log(3._r8)
         155           0 :    super(:)= 0.01_r8*supersat(:)
         156             :
         157             :    allocate(                  &
         158           0 :       alogsig(naer_all),      &
         159           0 :       exp45logsig(naer_all),  &
         160           0 :       argfactor(naer_all),    &
         161           0 :       f1(naer_all),           &
         162           0 :       f2(naer_all),           &
         163           0 :       amcubefactor(naer_all), &
         164           0 :       smcritfactor(naer_all), &
         165           0 :       amcube(naer_all),       &
         166           0 :       smcrit(naer_all),       &
         167           0 :       lnsm(naer_all),         &
         168           0 :       ccnfact(psat,naer_all)  )
         169             :
         170             :
         171           0 :    do m = 1, naer_all
         172             :
         173             :       ! Skip aerosols that don't have a dispersion defined.
         174           0 :       if (dispersion_aer(m) == 0._r8) cycle
         175             :
         176           0 :       alogsig(m)     = log(dispersion_aer(m))
         177           0 :       exp45logsig(m) = exp(4.5_r8*alogsig(m)*alogsig(m))
         178           0 :       argfactor(m)   = 2._r8/(3._r8*sqrt(2._r8)*alogsig(m))
         179           0 :       f1(m)          = 0.5_r8*exp(2.5_r8*alogsig(m)*alogsig(m))
         180           0 :       f2(m)          = 1._r8 + 0.25_r8*alogsig(m)
         181           0 :       amcubefactor(m)= 3._r8/(4._r8*pi*exp45logsig(m)*density_aer(m))
         182           0 :       smcritfactor(m)= 2._r8*aten*sqrt(aten/(27._r8*max(1.e-10_r8,hygro_aer(m))))
         183           0 :       amcube(m)      = amcubefactor(m)/num_to_mass_aer(m)
         184             :
         185           0 :       if (hygro_aer(m) .gt. 1.e-10_r8) then
         186           0 :          smcrit(m) = smcritfactor(m)/sqrt(amcube(m))
         187             :       else
         188           0 :          smcrit(m) = 100._r8
         189             :       endif
         190           0 :       lnsm(m) = log(smcrit(m))
         191             :
         192           0 :       do l = 1, psat
         193           0 :          arg = argfactor(m)*log(smcrit(m)/super(l))
         194           0 :          if (arg < 2) then
         195           0 :             if (arg < -2) then
         196           0 :                ccnfact(l,m) = 1.e-6_r8
         197             :             else
         198           0 :                ccnfact(l,m) = 1.e-6_r8*0.5_r8*erfc(arg)
         199             :             endif
         200             :          else
         201           0 :             ccnfact(l,m) = 0._r8
         202             :          endif
         203             :       enddo
         204             :
         205             :    end do
         206             :
         207           0 : end subroutine ndrop_bam_init
         208             :
         209             : !===============================================================================
         210             :
         211           0 : subroutine ndrop_bam_run( &
         212           0 :    wbar, tair, rhoair, na, pmode, &
         213           0 :    nmode, ma, nact)
         214             :
         215             :    ! calculates number fraction of aerosols activated as CCN
         216             :    ! assumes an internal mixture within each of up to pmode multiple aerosol modes
         217             :    ! a gaussian spectrum of updrafts can be treated.
         218             :
         219             :    !      mks units
         220             :
         221             :    !      Abdul-Razzak and Ghan, A parameterization of aerosol activation.
         222             :    !      2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
         223             :
         224             :    ! input
         225             :    integer,  intent(in) :: pmode         ! dimension of modes
         226             :    integer,  intent(in) :: nmode         ! number of aerosol modes
         227             :    real(r8), intent(in) :: wbar          ! grid cell mean vertical velocity (m/s)
         228             :    real(r8), intent(in) :: tair          ! air temperature (K)
         229             :    real(r8), intent(in) :: rhoair        ! air density (kg/m3)
         230             :    real(r8), intent(in) :: na(pmode)     ! aerosol number concentration (1/m3)
         231             :    real(r8), intent(in) :: ma(pmode)     ! aerosol mass concentration (kg/m3)
         232             :
         233             :    ! output
         234             :    real(r8), intent(out) :: nact         ! number fraction of aerosols activated
         235             :
         236             :    ! local variables
         237             :    integer :: maxmodes
         238             :
         239           0 :    real(r8), allocatable :: volc(:) ! total aerosol volume  concentration (m3/m3)
         240           0 :    real(r8), allocatable :: eta(:)
         241           0 :    real(r8), allocatable :: smc(:)
         242           0 :    real(r8), allocatable :: etafactor2(:)
         243           0 :    real(r8), allocatable :: amcubeloc(:)
         244           0 :    real(r8), allocatable :: lnsmloc(:)
         245             :
         246             :    real(r8) :: pres     ! pressure (Pa)
         247             :    real(r8) :: diff0
         248             :    real(r8) :: conduct0 ! thermal conductivity (Joule/m/sec/deg)
         249             :    real(r8) :: qs       ! water vapor saturation mixing ratio
         250             :    real(r8) :: dqsdt    ! change in qs with temperature
         251             :    real(r8) :: gloc     ! thermodynamic function (m2/s)
         252             :    real(r8) :: zeta
         253             :    real(r8) :: lnsmax   ! ln(smax)
         254             :    real(r8) :: alpha
         255             :    real(r8) :: gammaloc
         256             :    real(r8) :: beta
         257             :    real(r8) :: sqrtg
         258             :    real(r8) :: wnuc
         259             :    real(r8) :: alw
         260             :    real(r8) :: sqrtalw
         261             :    real(r8) :: smax
         262             :    real(r8) :: x
         263             :    real(r8) :: etafactor1
         264             :    real(r8) :: etafactor2max
         265             :    real(r8) :: es
         266             :    integer  :: m
         267             :    !-------------------------------------------------------------------------------
         268             :
         269           0 :    maxmodes = naer_all
         270             :    allocate( &
         271             :       volc(maxmodes),       &
         272             :       eta(maxmodes),        &
         273             :       smc(maxmodes),        &
         274             :       etafactor2(maxmodes), &
         275             :       amcubeloc(maxmodes),  &
         276           0 :       lnsmloc(maxmodes)     )
         277             :
         278           0 :    if (maxmodes < pmode) then
         279           0 :       write(iulog,*)'ndrop_bam_run: maxmodes,pmode=',maxmodes,pmode
         280           0 :       call endrun('ndrop_bam_run')
         281             :    endif
         282             :
         283           0 :    nact = 0._r8
         284             :
         285           0 :    if (nmode .eq. 1  .and.  na(1) .lt. 1.e-20_r8) return
         286             :
         287           0 :    if (wbar .le. 0._r8) return
         288             :
         289           0 :    pres     = rair*rhoair*tair
         290           0 :    diff0    = 0.211e-4_r8*(pref/pres)*(tair/tmelt)**1.94_r8
         291           0 :    conduct0 = (5.69_r8 + 0.017_r8*(tair-tmelt))*4.186e2_r8*1.e-5_r8 ! convert to J/m/s/deg
         292           0 :    call qsat(tair, pres, es, qs)
         293           0 :    dqsdt    = latvap/(rh2o*tair*tair)*qs
         294           0 :    alpha    = gravit*(latvap/(cpair*rh2o*tair*tair) - 1._r8/(rair*tair))
         295           0 :    gammaloc = (1 + latvap/cpair*dqsdt)/(rhoair*qs)
         296             :    ! growth coefficent Abdul-Razzak & Ghan 1998 eqn 16
         297             :    ! should depend on mean radius of mode to account for gas kinetic effects
         298             :    gloc     = 1._r8/(rhoh2o/(diff0*rhoair*qs)                               &
         299           0 :               + latvap*rhoh2o/(conduct0*tair)*(latvap/(rh2o*tair) - 1._r8))
         300           0 :    sqrtg    = sqrt(gloc)
         301           0 :    beta     = 4._r8*pi*rhoh2o*gloc*gammaloc
         302           0 :    etafactor2max = 1.e10_r8/(alpha*wbar)**1.5_r8 ! this should make eta big if na is very small.
         303             :
         304           0 :    do m = 1, nmode
         305             :       ! skip aerosols with no dispersion, since they aren't meant to be CCN
         306           0 :       if (dispersion_aer(m) == 0._r8) then
         307           0 :          smc(m)=100._r8
         308           0 :          cycle
         309             :       end if
         310             :       ! internal mixture of aerosols
         311           0 :       volc(m) = ma(m)/(density_aer(m)) ! only if variable size dist
         312           0 :       if (volc(m) > 1.e-39_r8 .and. na(m) > 1.e-39_r8) then
         313           0 :          etafactor2(m) = 1._r8/(na(m)*beta*sqrtg)  !fixed or variable size dist
         314             :          ! number mode radius (m)
         315           0 :          amcubeloc(m) = (3._r8*volc(m)/(4._r8*pi*exp45logsig(m)*na(m)))  ! only if variable size dist
         316           0 :          smc(m) = smcrit(m) ! only for prescribed size dist
         317             :
         318           0 :          if (hygro_aer(m) > 1.e-10_r8) then   ! loop only if variable size dist
         319           0 :             smc(m) = 2._r8*aten*sqrt(aten/(27._r8*hygro_aer(m)*amcubeloc(m)))
         320             :          else
         321           0 :             smc(m) = 100._r8
         322             :          endif
         323             :       else
         324           0 :          smc(m) = 1._r8
         325           0 :          etafactor2(m) = etafactor2max ! this should make eta big if na is very small.
         326             :       end if
         327           0 :       lnsmloc(m) = log(smc(m)) ! only if variable size dist
         328             :    end do
         329             :
         330             :    ! single  updraft
         331           0 :    wnuc    = wbar
         332           0 :    alw     = alpha*wnuc
         333           0 :    sqrtalw = sqrt(alw)
         334           0 :    zeta    = 2._r8*sqrtalw*aten/(3._r8*sqrtg)
         335           0 :    etafactor1 = 2._r8*alw*sqrtalw
         336             :
         337           0 :    do m = 1, nmode
         338             :       ! skip aerosols with no dispersion, since they aren't meant to be CCN
         339           0 :       if (dispersion_aer(m) /= 0._r8) eta(m) = etafactor1*etafactor2(m)
         340             :    end do
         341             :
         342           0 :    call maxsat(zeta, eta, nmode, smc, smax)
         343           0 :    lnsmax = log(smax)
         344             :
         345           0 :    nact = 0._r8
         346           0 :    do m = 1, nmode
         347             :       ! skip aerosols with no dispersion, since they aren't meant to be CCN
         348           0 :       if (dispersion_aer(m) == 0._r8) cycle
         349           0 :       x = 2*(lnsmloc(m) - lnsmax)/(3*sq2*alogsig(m))
         350           0 :       nact = nact + 0.5_r8*(1._r8 - erf(x))*na(m)
         351             :    end do
         352           0 :    nact = nact/rhoair ! convert from #/m3 to #/kg
         353             :
         354           0 :    deallocate( &
         355             :       volc,       &
         356           0 :       eta,        &
         357           0 :       smc,        &
         358           0 :       etafactor2, &
         359           0 :       amcubeloc,  &
         360           0 :       lnsmloc     )
         361             :
         362           0 : end subroutine ndrop_bam_run
         363             :
         364             : !===============================================================================
         365             :
         366           0 : subroutine ndrop_bam_ccn(lchnk, ncol, maerosol, naer2)
         367             :
         368             :    !-------------------------------------------------------------------------------
         369             :    !
         370             :    ! Write diagnostic bulk aerosol ccn concentration
         371             :    !
         372             :    !-------------------------------------------------------------------------------
         373             :
         374             :    ! Input arguments
         375             :    integer,  intent(in) :: lchnk           ! chunk identifier
         376             :    integer,  intent(in) :: ncol            ! number of columns
         377             :    real(r8), intent(in) :: naer2(:,:,:)    ! aerosol number concentration (1/m3)
         378             :    real(r8), intent(in) :: maerosol(:,:,:) ! aerosol mass conc (kg/m3)
         379             :
         380             :    ! Local variables
         381             :    integer :: i, k, l, m
         382             :    real(r8) :: arg
         383             :
         384             :    character*8, parameter :: ccn_name(psat)=(/'CCN1','CCN2','CCN3','CCN4','CCN5','CCN6'/)
         385             :    real(r8) :: amcubesulfate(pcols)  ! cube of dry mode radius (m) of sulfate
         386             :    real(r8) :: smcritsulfate(pcols)  ! critical supersatuation for activation of sulfate
         387             :    real(r8) :: ccnfactsulfate
         388             :    real(r8) :: ccn(pcols,pver,psat)  ! number conc of aerosols activated at supersat
         389             :    !-------------------------------------------------------------------------------
         390             :
         391           0 :    ccn(:ncol,:,:) = 0._r8
         392             :
         393           0 :    do k = top_lev, pver
         394             :
         395           0 :       do m = 1, naer_all
         396             :
         397           0 :          if (m == idxsul) then
         398             :             ! Lohmann treatment for sulfate has variable size distribution
         399           0 :             do i = 1, ncol
         400           0 :                if (naer2(i,k,m) > 0._r8) then
         401           0 :                   amcubesulfate(i) = amcubefactor(m)*maerosol(i,k,m)/(naer2(i,k,m))
         402           0 :                   smcritsulfate(i) = smcritfactor(m)/sqrt(amcubesulfate(i))
         403             :                else
         404           0 :                   smcritsulfate(i) = 1._r8
         405             :                end if
         406             :             end do
         407             :          end if
         408             :
         409           0 :          do l = 1, psat
         410             :
         411           0 :             if (m == idxsul) then
         412             :                ! This code is modifying ccnfact for sulfate only.
         413           0 :                do i = 1, ncol
         414           0 :                   arg = argfactor(m)*log(smcritsulfate(i)/super(l))
         415           0 :                   if (arg < 2) then
         416           0 :                      if (arg < -2) then
         417             :                         ccnfactsulfate = 1.0e-6_r8
         418             :                      else
         419           0 :                         ccnfactsulfate = 0.5e-6_r8*erfc(arg)
         420             :                      end if
         421             :                   else
         422             :                      ccnfactsulfate = 0.0_r8
         423             :                   end if
         424           0 :                   ccn(i,k,l) = ccn(i,k,l) + naer2(i,k,m)*ccnfactsulfate
         425             :                end do
         426             :             else
         427             :                ! Non-sulfate species use ccnfact computed by the init routine
         428           0 :                ccn(:ncol,k,l) = ccn(:ncol,k,l) + naer2(:ncol,k,m)*ccnfact(l,m)
         429             :             end if
         430             :
         431             :          end do   ! supersaturation
         432             :       end do      ! bulk aerosol
         433             :    end do         ! level
         434             :
         435           0 :    do l = 1, psat
         436           0 :       call outfld(ccn_name(l), ccn(1,1,l), pcols, lchnk)
         437             :    end do
         438             :
         439           0 :    do l = 1, naer_all
         440           0 :       call outfld(trim(aername(l))//'_m3', naer2(:,:,l), pcols, lchnk)
         441             :    end do
         442             :
         443           0 : end subroutine ndrop_bam_ccn
         444             :
         445             : !===============================================================================
         446             :
         447           0 : subroutine maxsat(zeta, eta, nmode, smc, smax)
         448             :
         449             :    ! calculates maximum supersaturation for multiple
         450             :    ! competing aerosol modes.
         451             :
         452             :    ! Abdul-Razzak and Ghan, A parameterization of aerosol activation.
         453             :    ! 2. Multiple aerosol types. J. Geophys. Res., 105, 6837-6844.
         454             :
         455             :    real(r8), intent(in) :: zeta
         456             :    integer,  intent(in) :: nmode ! number of modes
         457             :    real(r8), intent(in) :: smc(:) ! critical supersaturation for number mode radius
         458             :    real(r8), intent(in) :: eta(:)
         459             :
         460             :    real(r8), intent(out) :: smax ! maximum supersaturation
         461             :
         462             :    integer :: m  ! mode index
         463             :    real(r8) :: sum, g1, g2
         464             :
         465           0 :    do m=1,nmode
         466           0 :       if(zeta.gt.1.e5_r8*eta(m).or.smc(m)*smc(m).gt.1.e5_r8*eta(m))then
         467             :          ! weak forcing. essentially none activated
         468           0 :          smax=1.e-20_r8
         469             :       else
         470             :          ! significant activation of this mode. calc activation all modes.
         471             :          go to 1
         472             :       endif
         473             :    enddo
         474             :
         475           0 :    return
         476             :
         477             : 1  continue
         478             :
         479             :    sum=0
         480           0 :    do m=1,nmode
         481           0 :       if(eta(m).gt.1.e-20_r8)then
         482           0 :          g1=sqrt(zeta/eta(m))
         483           0 :          g1=g1*g1*g1
         484           0 :          g2=smc(m)/sqrt(eta(m)+3*zeta)
         485           0 :          g2=sqrt(g2)
         486           0 :          g2=g2*g2*g2
         487           0 :          sum=sum+(f1(m)*g1+f2(m)*g2)/(smc(m)*smc(m))
         488             :       else
         489             :          sum=1.e20_r8
         490             :       endif
         491             :    enddo
         492             :
         493           0 :    smax=1._r8/sqrt(sum)
         494             :
         495             : end subroutine maxsat
         496             :
         497             : !===============================================================================
         498             :
         499             : end module ndrop_bam

Generated by: LCOV version 1.14