LCOV - code coverage report
Current view: top level - icepack/columnphysics - icepack_atmo.F90 (source / functions) Hit Total Coverage
Test: 210526-003153:a09e51fb16:8:first,base,travis,decomp,reprosum,io,quick,unittest Lines: 240 250 96.00 %
Date: 2021-05-26 02:35:17 Functions: 6 6 100.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : 
       3             : ! Atmospheric boundary interface (stability based flux calculations)
       4             : 
       5             : ! author: Elizabeth C. Hunke, LANL
       6             : !
       7             : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
       8             : ! 2004: Block structure added by William Lipscomb
       9             : ! 2006: Converted to free source form (F90) by Elizabeth Hunke
      10             : ! 2013: Form drag routine added (neutral_drag_coeffs) by David Schroeder 
      11             : ! 2014: Adjusted form drag and added high frequency coupling by Andrew Roberts
      12             : 
      13             :       module icepack_atmo
      14             : 
      15             :       use icepack_kinds
      16             :       use icepack_parameters,  only: c0, c1, c2, c4, c5, c8, c10
      17             :       use icepack_parameters,  only: c16, c20, p001, p01, p2, p4, p5, p75, puny
      18             :       use icepack_parameters,  only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
      19             :       use icepack_parameters,  only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
      20             :       use icepack_parameters,  only: pih, dragio, rhoi, rhos, rhow
      21             :       use icepack_parameters, only: atmbndy, calc_strair, formdrag
      22             :       use icepack_tracers, only: n_iso
      23             :       use icepack_tracers, only: tr_iso
      24             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      25             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      26             : 
      27             :       implicit none
      28             : 
      29             :       private
      30             :       public :: neutral_drag_coeffs, &
      31             :                 icepack_atm_boundary
      32             : 
      33             : !=======================================================================
      34             : 
      35             :       contains
      36             : 
      37             : !=======================================================================
      38             : 
      39             : ! Compute coefficients for atm/ice fluxes, stress, and reference
      40             : ! temperature and humidity. NOTE:
      41             : ! (1) All fluxes are positive downward, 
      42             : ! (2) Here, tstar = (WT)/U*, and qstar = (WQ)/U*,
      43             : ! (3a) wind speeds should all be above a minimum speed (eg. 1.0 m/s).
      44             : !
      45             : ! ASSUME:
      46             : !  The saturation humidity of air at T(K): qsat(T)  (kg/m**3)
      47             : !
      48             : ! Code originally based on CSM1
      49             : 
      50  2137677905 :       subroutine atmo_boundary_layer (sfctype,            &
      51             :                                       calc_strair, formdrag, &   ! LCOV_EXCL_LINE
      52             :                                       Tsf,      potT,     &   ! LCOV_EXCL_LINE
      53             :                                       uatm,     vatm,     &     ! LCOV_EXCL_LINE
      54             :                                       wind,     zlvl,     &     ! LCOV_EXCL_LINE
      55             :                                       Qa,       rhoa,     &   ! LCOV_EXCL_LINE
      56             :                                       strx,     stry,     &      ! LCOV_EXCL_LINE
      57             :                                       Tref,     Qref,     &   ! LCOV_EXCL_LINE
      58             :                                       delt,     delq,     &   ! LCOV_EXCL_LINE
      59             :                                       lhcoef,   shcoef,   &   ! LCOV_EXCL_LINE
      60             :                                       Cdn_atm,            &   ! LCOV_EXCL_LINE
      61             :                                       Cdn_atm_ratio_n,    &   ! LCOV_EXCL_LINE
      62             :                                       Qa_iso,   Qref_iso, &   ! LCOV_EXCL_LINE
      63             :                                       iso_flag,           &   ! LCOV_EXCL_LINE
      64             :                                       uvel,     vvel,     &   ! LCOV_EXCL_LINE
      65             :                                       Uref                )     
      66             : 
      67             :       use icepack_parameters, only: highfreq, natmiter, atmiter_conv
      68             : 
      69             :       character (len=3), intent(in) :: &
      70             :          sfctype      ! ice or ocean
      71             : 
      72             :       logical (kind=log_kind), intent(in) :: &
      73             :          calc_strair, &  ! if true, calculate wind stress components   ! LCOV_EXCL_LINE
      74             :          formdrag        ! if true, calculate form drag
      75             : 
      76             :       real (kind=dbl_kind), intent(in) :: &
      77             :          Tsf      , & ! surface temperature of ice or ocean   ! LCOV_EXCL_LINE
      78             :          potT     , & ! air potential temperature  (K)   ! LCOV_EXCL_LINE
      79             :          uatm     , & ! x-direction wind speed (m/s)   ! LCOV_EXCL_LINE
      80             :          vatm     , & ! y-direction wind speed (m/s)   ! LCOV_EXCL_LINE
      81             :          wind     , & ! wind speed (m/s)   ! LCOV_EXCL_LINE
      82             :          zlvl     , & ! atm level height (m)   ! LCOV_EXCL_LINE
      83             :          Qa       , & ! specific humidity (kg/kg)   ! LCOV_EXCL_LINE
      84             :          rhoa         ! air density (kg/m^3)
      85             : 
      86             :       real (kind=dbl_kind), intent(inout) :: &
      87             :          Cdn_atm      ! neutral drag coefficient
      88             :  
      89             :       real (kind=dbl_kind), intent(inout) :: &
      90             :          Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff
      91             : 
      92             :       real (kind=dbl_kind), &
      93             :          intent(inout) :: &   ! LCOV_EXCL_LINE
      94             :          strx     , & ! x surface stress (N)   ! LCOV_EXCL_LINE
      95             :          stry         ! y surface stress (N)
      96             : 
      97             :       real (kind=dbl_kind), intent(inout) :: &
      98             :          Tref     , & ! reference height temperature  (K)   ! LCOV_EXCL_LINE
      99             :          Qref     , & ! reference height specific humidity (kg/kg)   ! LCOV_EXCL_LINE
     100             :          delt     , & ! potential T difference   (K)   ! LCOV_EXCL_LINE
     101             :          delq     , & ! humidity difference      (kg/kg)   ! LCOV_EXCL_LINE
     102             :          shcoef   , & ! transfer coefficient for sensible heat   ! LCOV_EXCL_LINE
     103             :          lhcoef       ! transfer coefficient for latent heat
     104             : 
     105             :       logical (kind=log_kind), intent(in), optional :: &
     106             :          iso_flag     ! flag to trigger iso calculations
     107             : 
     108             :       real (kind=dbl_kind), intent(in), optional, dimension(:) :: &
     109             :          Qa_iso       ! specific isotopic humidity (kg/kg)
     110             : 
     111             :       real (kind=dbl_kind), intent(inout), optional, dimension(:) :: &
     112             :          Qref_iso     ! reference specific isotopic humidity (kg/kg)
     113             : 
     114             :       real (kind=dbl_kind), intent(in) :: &
     115             :          uvel     , & ! x-direction ice speed (m/s)   ! LCOV_EXCL_LINE
     116             :          vvel         ! y-direction ice speed (m/s)
     117             : 
     118             :       real (kind=dbl_kind), intent(out) :: &
     119             :          Uref         ! reference height wind speed (m/s)
     120             : 
     121             :       ! local variables
     122             : 
     123             :       integer (kind=int_kind) :: &
     124             :          k,n         ! iteration index
     125             : 
     126             :       real (kind=dbl_kind) :: &
     127             :          TsfK  , & ! surface temperature in Kelvin (K)   ! LCOV_EXCL_LINE
     128             :          xqq   , & ! temporary variable   ! LCOV_EXCL_LINE
     129             :          psimh , & ! stability function at zlvl   (momentum)   ! LCOV_EXCL_LINE
     130             :          tau   , & ! stress at zlvl   ! LCOV_EXCL_LINE
     131             :          fac   , & ! interpolation factor   ! LCOV_EXCL_LINE
     132             :          al2   , & ! ln(z10   /zTrf)   ! LCOV_EXCL_LINE
     133             :          psix2 , & ! stability function at zTrf   (heat and water)   ! LCOV_EXCL_LINE
     134             :          psimhs, & ! stable profile   ! LCOV_EXCL_LINE
     135             :          ssq   , & ! sat surface humidity     (kg/kg)   ! LCOV_EXCL_LINE
     136             :          qqq   , & ! for qsat, dqsfcdt   ! LCOV_EXCL_LINE
     137             :          TTT   , & ! for qsat, dqsfcdt   ! LCOV_EXCL_LINE
     138             :          qsat  , & ! the saturation humidity of air (kg/m^3)   ! LCOV_EXCL_LINE
     139             :          Lheat , & ! Lvap or Lsub, depending on surface type   ! LCOV_EXCL_LINE
     140   216108905 :          umin      ! minimum wind speed (m/s)
     141             : 
     142             :       real (kind=dbl_kind) :: &
     143             :          ustar , & ! ustar (m/s)   ! LCOV_EXCL_LINE
     144             :          ustar_prev , & ! ustar_prev (m/s)   ! LCOV_EXCL_LINE
     145             :          vscl  , & ! vscl   ! LCOV_EXCL_LINE
     146             :          tstar , & ! tstar   ! LCOV_EXCL_LINE
     147             :          qstar , & ! qstar   ! LCOV_EXCL_LINE
     148             :          ratio , & ! ratio   ! LCOV_EXCL_LINE
     149             :          rdn   , & ! sqrt of neutral exchange coefficient (momentum)   ! LCOV_EXCL_LINE
     150             :          rhn   , & ! sqrt of neutral exchange coefficient (heat)   ! LCOV_EXCL_LINE
     151             :          ren   , & ! sqrt of neutral exchange coefficient (water)   ! LCOV_EXCL_LINE
     152             :          rd    , & ! sqrt of exchange coefficient (momentum)   ! LCOV_EXCL_LINE
     153             :          re    , & ! sqrt of exchange coefficient (water)   ! LCOV_EXCL_LINE
     154             :          rh    , & ! sqrt of exchange coefficient (heat)   ! LCOV_EXCL_LINE
     155             :          vmag  , & ! surface wind magnitude   (m/s)   ! LCOV_EXCL_LINE
     156             :          alz   , & ! ln(zlvl  /z10)   ! LCOV_EXCL_LINE
     157             :          thva  , & ! virtual temperature      (K)   ! LCOV_EXCL_LINE
     158             :          cp    , & ! specific heat of moist air   ! LCOV_EXCL_LINE
     159             :          hol   , & ! H (at zlvl  ) over L   ! LCOV_EXCL_LINE
     160             :          stable, & ! stability factor   ! LCOV_EXCL_LINE
     161             :          cpvir , & ! defined as cp_wv/cp_air - 1.   ! LCOV_EXCL_LINE
     162   216108905 :          psixh     ! stability function at zlvl   (heat and water)
     163             : 
     164             :       real (kind=dbl_kind), parameter :: &
     165             :          zTrf  = c2     ! reference height for air temp (m)
     166             : 
     167             :       logical (kind=log_kind) :: &
     168             :          l_iso_flag     ! local flag to trigger iso calculations
     169             : 
     170             :       character(len=*),parameter :: subname='(atmo_boundary_layer)'
     171             : 
     172  2137677905 :       l_iso_flag = .false.
     173  2137677905 :       if (present(iso_flag)) then
     174  2137677905 :         l_iso_flag = iso_flag
     175             :       endif
     176             : 
     177  2137677905 :       al2 = log(zref/zTrf)
     178             : 
     179             :       !------------------------------------------------------------
     180             :       ! Initialize
     181             :       !------------------------------------------------------------
     182             : 
     183  2137677905 :       cpvir = cp_wv/cp_air-c1   ! defined as cp_wv/cp_air - 1.
     184             : 
     185  2137677905 :       if (highfreq) then       
     186   128464665 :        umin  = p5 ! minumum allowable wind-ice speed difference of 0.5 m/s
     187             :       else
     188  2009213240 :        umin  = c1 ! minumum allowable wind speed of 1m/s
     189             :       endif
     190             : 
     191  2137677905 :       Tref = c0
     192  2137677905 :       Qref = c0
     193  2137677905 :       Uref = c0
     194  2137677905 :       delt = c0
     195  2137677905 :       delq = c0
     196  2137677905 :       shcoef = c0
     197  2137677905 :       lhcoef = c0
     198             : 
     199             :       !------------------------------------------------------------
     200             :       ! Compute turbulent flux coefficients, wind stress, and
     201             :       ! reference temperature and humidity.
     202             :       !------------------------------------------------------------
     203             : 
     204             :       !------------------------------------------------------------
     205             :       ! define variables that depend on surface type
     206             :       !------------------------------------------------------------
     207             : 
     208  2137677905 :       if (sfctype(1:3)=='ice') then
     209             : 
     210   975413849 :          qqq   = qqqice          ! for qsat
     211   975413849 :          TTT   = TTTice          ! for qsat
     212   975413849 :          Lheat = Lsub            ! ice to vapor
     213             : 
     214   975413849 :             if (highfreq) then
     215             :                vmag = max(umin, sqrt( (uatm-uvel)**2 + &
     216    98806785 :                                       (vatm-vvel)**2) )
     217             :             else
     218   876607064 :                vmag = max(umin, wind)
     219             :             endif
     220             : 
     221   975413849 :             if (formdrag .and. Cdn_atm > puny) then 
     222    23322166 :                rdn = sqrt(Cdn_atm)               
     223             :             else
     224   952091683 :                rdn  = vonkar/log(zref/iceruf) ! neutral coefficient
     225   952091683 :                Cdn_atm = rdn * rdn
     226             :             endif
     227             : 
     228  1162264056 :       elseif (sfctype(1:3)=='ocn') then
     229             : 
     230  1162264056 :          qqq   = qqqocn
     231  1162264056 :          TTT   = TTTocn
     232  1162264056 :          Lheat = Lvap           ! liquid to vapor
     233  1162264056 :          vmag = max(umin, wind)
     234             :          rdn  = sqrt(0.0027_dbl_kind/vmag &
     235  1162264056 :               + .000142_dbl_kind + .0000764_dbl_kind*vmag)
     236             : 
     237             :       endif   ! sfctype
     238             : 
     239             :       !------------------------------------------------------------
     240             :       ! define some more needed variables
     241             :       !------------------------------------------------------------
     242             : 
     243  2137677905 :       TsfK = Tsf + Tffresh     ! surface temp (K)
     244  2137677905 :       delt = potT - TsfK       ! pot temp diff (K)
     245  2137677905 :       qsat = qqq * exp(-TTT/TsfK)   ! saturation humidity (kg/m^3)
     246  2137677905 :       ssq  = qsat / rhoa       ! sat surf hum (kg/kg)
     247             :       
     248  2137677905 :       thva = potT * (c1 + zvir * Qa) ! virtual pot temp (K)
     249  2137677905 :       delq = Qa - ssq          ! spec hum dif (kg/kg)
     250  2137677905 :       alz  = log(zlvl/zref)
     251  2137677905 :       cp   = cp_air*(c1 + cpvir*ssq)
     252             :       
     253             :       !------------------------------------------------------------
     254             :       ! first estimate of Z/L and ustar, tstar and qstar
     255             :       !------------------------------------------------------------
     256             : 
     257             :       ! neutral coefficients, z/L = 0.0
     258  2137677905 :       rhn = rdn
     259  2137677905 :       ren = rdn
     260             :       
     261             :       ! ustar,tstar,qstar
     262  2137677905 :       ustar = rdn * vmag
     263  2137677905 :       tstar = rhn * delt
     264  2137677905 :       qstar = ren * delq
     265             : 
     266             :       !------------------------------------------------------------
     267             :       ! iterate to converge on Z/L, ustar, tstar and qstar
     268             :       !------------------------------------------------------------
     269             : 
     270  2137677905 :       ustar_prev = c2 * ustar
     271             : 
     272  2137677905 :       k = 1
     273 13050468559 :       do while (abs(ustar - ustar_prev)/ustar > atmiter_conv .and. k <= natmiter)
     274 10912790654 :          k = k + 1
     275 10912790654 :          ustar_prev = ustar
     276             : 
     277             :          ! compute stability & evaluate all stability functions
     278             :          hol = vonkar * gravit * zlvl &   ! LCOV_EXCL_LINE
     279             :                  * (tstar/thva &   ! LCOV_EXCL_LINE
     280             :                  + qstar/(c1/zvir+Qa)) &   ! LCOV_EXCL_LINE
     281 10912790654 :                  / ustar**2
     282 10912790654 :          hol    = sign( min(abs(hol),c10), hol )
     283 10912790654 :          stable = p5 + sign(p5 , hol)
     284 10912790654 :          xqq    = max(sqrt(abs(c1 - c16*hol)) , c1)
     285 10912790654 :          xqq    = sqrt(xqq)
     286             : 
     287             :          ! Jordan et al 1999
     288             :          psimhs = -(0.7_dbl_kind*hol &
     289             :                 + 0.75_dbl_kind*(hol-14.3_dbl_kind) &   ! LCOV_EXCL_LINE
     290 10912790654 :                 * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
     291             :          psimh  = psimhs*stable &
     292 10912790654 :                 + (c1 - stable)*psimhu(xqq)
     293             :          psixh  = psimhs*stable &
     294 10912790654 :                 + (c1 - stable)*psixhu(xqq)
     295             : 
     296             :          ! shift all coeffs to measurement height and stability
     297 10912790654 :          rd = rdn / (c1+rdn/vonkar*(alz-psimh))
     298 10912790654 :          rh = rhn / (c1+rhn/vonkar*(alz-psixh))
     299 10912790654 :          re = ren / (c1+ren/vonkar*(alz-psixh))
     300             :       
     301             :          ! update ustar, tstar, qstar using updated, shifted coeffs
     302 10912790654 :          ustar = rd * vmag
     303 10912790654 :          tstar = rh * delt
     304 13050468559 :          qstar = re * delq
     305             : 
     306             :       enddo                     ! end iteration
     307             : 
     308  2137677905 :       if (calc_strair) then
     309             : 
     310             :          ! initialize
     311  2137677905 :          strx = c0
     312  2137677905 :          stry = c0
     313             : 
     314  2137677905 :          if (highfreq .and. sfctype(1:3)=='ice') then
     315             : 
     316             :             !------------------------------------------------------------
     317             :             ! momentum flux for high frequency coupling (RASM/CESM) 
     318             :             !------------------------------------------------------------
     319             :             ! tau = rhoa * rd * rd
     320             :             ! strx = tau * |Uatm-U| * (uatm-u)
     321             :             ! stry = tau * |Uatm-U| * (vatm-v)
     322             :             !------------------------------------------------------------
     323             : 
     324    98806785 :             tau = rhoa * rd * rd ! not the stress at zlvl
     325             : 
     326             :             ! high frequency momentum coupling following Roberts et al. (2014)
     327    98806785 :             strx = tau * sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * (uatm-uvel)
     328    98806785 :             stry = tau * sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * (vatm-vvel)
     329             : 
     330             :          else
     331             : 
     332             :             !------------------------------------------------------------
     333             :             ! momentum flux
     334             :             !------------------------------------------------------------
     335             :             ! tau = rhoa * ustar * ustar
     336             :             ! strx = tau * uatm / vmag
     337             :             ! stry = tau * vatm / vmag
     338             :             !------------------------------------------------------------
     339             : 
     340  2038871120 :             tau = rhoa * ustar * rd ! not the stress at zlvl
     341  2038871120 :             strx = tau * uatm
     342  2038871120 :             stry = tau * vatm
     343             : 
     344             :          endif
     345             : 
     346  2137677905 :          Cdn_atm_ratio_n = rd * rd / rdn / rdn
     347             : 
     348             :       endif                     ! calc_strair
     349             : 
     350             :       !------------------------------------------------------------
     351             :       ! coefficients for turbulent flux calculation
     352             :       !------------------------------------------------------------
     353             :       ! add windless coefficient for sensible heat flux
     354             :       ! as in Jordan et al (JGR, 1999)
     355             :       !------------------------------------------------------------
     356             : 
     357  2137677905 :       shcoef = rhoa * ustar * cp * rh + c1
     358  2137677905 :       lhcoef = rhoa * ustar * Lheat  * re
     359             : 
     360             :       !------------------------------------------------------------
     361             :       ! Compute diagnostics: 2m ref T, Q, U
     362             :       !------------------------------------------------------------
     363             : 
     364  2137677905 :       hol   = hol*zTrf/zlvl
     365  2137677905 :       xqq   = max( c1, sqrt(abs(c1-c16*hol)) )
     366  2137677905 :       xqq   = sqrt(xqq)
     367  2137677905 :       psix2 = -c5*hol*stable + (c1-stable)*psixhu(xqq)
     368             :       fac   = (rh/vonkar) &
     369  2137677905 :             * (alz + al2 - psixh + psix2)
     370  2137677905 :       Tref  = potT - delt*fac
     371  2137677905 :       Tref  = Tref - p01*zTrf ! pot temp to temp correction
     372             :       fac   = (re/vonkar) &
     373  2137677905 :             * (alz + al2 - psixh + psix2)
     374  2137677905 :       Qref  = Qa - delq*fac
     375             : 
     376  2137677905 :       if (highfreq .and. sfctype(1:3)=='ice') then
     377    98806785 :          Uref = sqrt((uatm-uvel)**2 + (vatm-vvel)**2) * rd / rdn
     378             :       else
     379  2038871120 :          Uref = vmag * rd / rdn
     380             :       endif
     381             : 
     382  2137677905 :       if (l_iso_flag) then
     383   975413849 :        if (present(Qref_iso) .and. present(Qa_iso)) then
     384  1041260231 :          Qref_iso(:) = c0 
     385   975413849 :          if (tr_iso) then
     386  1063209025 :             do n = 1, n_iso
     387    65846382 :                ratio = c0
     388    65846382 :                if (Qa_iso(2) > puny) ratio = Qa_iso(n)/Qa_iso(2)
     389    87795176 :                Qref_iso(n) = Qa_iso(n) - ratio*delq*fac
     390             :             enddo
     391             :          endif
     392             :        else
     393           0 :          call icepack_warnings_add(subname//' l_iso_flag true but optional arrays missing')
     394           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     395           0 :          return
     396             :        endif
     397             :       endif
     398             : 
     399             :       end subroutine atmo_boundary_layer
     400             : 
     401             : !=======================================================================
     402             : 
     403             : ! Compute coefficients for atm/ice fluxes, stress
     404             : ! NOTE: \\
     405             : ! (1) all fluxes are positive downward,  \\
     406             : ! (2) reference temperature and humidity are NOT computed
     407             : 
     408     2079360 :       subroutine atmo_boundary_const (sfctype,  calc_strair, &
     409             :                                       uatm,     vatm,     &     ! LCOV_EXCL_LINE
     410             :                                       wind,     rhoa,     &   ! LCOV_EXCL_LINE
     411             :                                       strx,     stry,     &      ! LCOV_EXCL_LINE
     412             :                                       Tsf,      potT,     &   ! LCOV_EXCL_LINE
     413             :                                       Qa,                 &   ! LCOV_EXCL_LINE
     414             :                                       delt,     delq,     &   ! LCOV_EXCL_LINE
     415             :                                       lhcoef,   shcoef    )
     416             : 
     417             :       character (len=3), intent(in) :: &
     418             :          sfctype      ! ice or ocean
     419             : 
     420             :       logical (kind=log_kind), intent(in) :: &
     421             :          calc_strair  ! if true, calculate wind stress components
     422             : 
     423             :       real (kind=dbl_kind), intent(in) :: &
     424             :          Tsf      , & ! surface temperature of ice or ocean   ! LCOV_EXCL_LINE
     425             :          potT     , & ! air potential temperature  (K)   ! LCOV_EXCL_LINE
     426             :          Qa       , & ! specific humidity (kg/kg)   ! LCOV_EXCL_LINE
     427             :          uatm     , & ! x-direction wind speed (m/s)   ! LCOV_EXCL_LINE
     428             :          vatm     , & ! y-direction wind speed (m/s)   ! LCOV_EXCL_LINE
     429             :          wind     , & ! wind speed (m/s)   ! LCOV_EXCL_LINE
     430             :          rhoa         ! air density (kg/m^3)
     431             : 
     432             :       real (kind=dbl_kind), intent(inout):: &
     433             :          strx     , & ! x surface stress (N)   ! LCOV_EXCL_LINE
     434             :          stry         ! y surface stress (N)
     435             : 
     436             :       real (kind=dbl_kind), intent(out):: &
     437             :          delt     , & ! potential T difference   (K)   ! LCOV_EXCL_LINE
     438             :          delq     , & ! humidity difference      (kg/kg)   ! LCOV_EXCL_LINE
     439             :          shcoef   , & ! transfer coefficient for sensible heat   ! LCOV_EXCL_LINE
     440             :          lhcoef       ! transfer coefficient for latent heat
     441             : 
     442             :        ! local variables
     443             : 
     444             :       real (kind=dbl_kind) :: &
     445             :          TsfK, & ! surface temperature in Kelvin (K)   ! LCOV_EXCL_LINE
     446             :          qsat, & ! the saturation humidity of air (kg/m^3)   ! LCOV_EXCL_LINE
     447             :          ssq , & ! sat surface humidity     (kg/kg)   ! LCOV_EXCL_LINE
     448             :          tau, &  ! stress at zlvl   ! LCOV_EXCL_LINE
     449           0 :          Lheat   ! Lvap or Lsub, depending on surface type
     450             : 
     451             :       character(len=*),parameter :: subname='(atmo_boundary_const)'
     452             : 
     453             :       !------------------------------------------------------------
     454             :       ! Initialize
     455             :       !------------------------------------------------------------
     456             : 
     457     2079360 :       delt = c0
     458     2079360 :       delq = c0
     459     2079360 :       shcoef = c0
     460     2079360 :       lhcoef = c0
     461             :       
     462     2079360 :       if (calc_strair) then
     463             : 
     464           0 :          strx = c0
     465           0 :          stry = c0
     466             : 
     467             :       !------------------------------------------------------------
     468             :       ! momentum flux
     469             :       !------------------------------------------------------------
     470           0 :          tau = rhoa * 0.0012_dbl_kind * wind
     471             : !AOMIP         tau = rhoa * (1.10_dbl_kind + c4*p01*wind) &
     472             : !AOMIP                         * wind * p001
     473           0 :          strx = tau * uatm
     474           0 :          stry = tau * vatm
     475             : 
     476             :       endif                     ! calc_strair
     477             : 
     478             :       !------------------------------------------------------------
     479             :       ! define variables that depend on surface type
     480             :       !------------------------------------------------------------
     481             : 
     482     2079360 :       if (sfctype(1:3)=='ice') then
     483           0 :          Lheat = Lsub           ! ice to vapor
     484     2079360 :       elseif (sfctype(1:3)=='ocn') then
     485     2079360 :          Lheat = Lvap           ! liquid to vapor
     486             :       endif   ! sfctype
     487             : 
     488             :       !------------------------------------------------------------
     489             :       ! potential temperature and specific humidity differences
     490             :       !------------------------------------------------------------
     491             : 
     492     2079360 :       TsfK     = Tsf + Tffresh    ! surface temp (K)
     493     2079360 :       qsat     = qqqocn * exp(-TTTocn/TsfK) ! sat humidity (kg/m^3)
     494     2079360 :       ssq      = qsat / rhoa      ! sat surf hum (kg/kg)
     495             :       
     496     2079360 :       delt= potT - TsfK      ! pot temp diff (K)
     497     2079360 :       delq= Qa - ssq         ! spec hum dif (kg/kg)
     498             :       
     499             :       !------------------------------------------------------------
     500             :       ! coefficients for turbulent flux calculation
     501             :       !------------------------------------------------------------
     502             : 
     503     2079360 :       shcoef = (1.20e-3_dbl_kind)*cp_air*rhoa*wind
     504     2079360 :       lhcoef = (1.50e-3_dbl_kind)*Lheat *rhoa*wind
     505             : 
     506     2079360 :       end subroutine atmo_boundary_const
     507             : 
     508             : !=======================================================================
     509             : 
     510             : ! Neutral drag coefficients for ocean and atmosphere also compute the 
     511             : ! intermediate necessary variables ridge height, distance, floe size
     512             : ! based upon Tsamados et al. (2014), JPO, DOI: 10.1175/JPO-D-13-0215.1.
     513             : ! Places where the code varies from the paper are commented.
     514             : !
     515             : ! authors: Michel Tsamados, CPOM
     516             : !          David Schroeder, CPOM
     517             : !
     518             : ! changes: Andrew Roberts, NPS (RASM/CESM coupling and documentation)
     519             : 
     520    20367264 :       subroutine neutral_drag_coeffs (apnd,     hpnd,     &
     521             :                                       ipnd,               &   ! LCOV_EXCL_LINE
     522             :                                       alvl,     vlvl,     &   ! LCOV_EXCL_LINE
     523             :                                       aice,     vice,     &   ! LCOV_EXCL_LINE
     524             :                                       vsno,     aicen,    &   ! LCOV_EXCL_LINE
     525             :                                       vicen, &   ! LCOV_EXCL_LINE
     526             :                                       Cdn_ocn,  Cdn_ocn_skin,    &   ! LCOV_EXCL_LINE
     527             :                                       Cdn_ocn_floe, Cdn_ocn_keel,&   ! LCOV_EXCL_LINE
     528             :                                       Cdn_atm,  Cdn_atm_skin,    &   ! LCOV_EXCL_LINE
     529             :                                       Cdn_atm_floe, Cdn_atm_pond,&   ! LCOV_EXCL_LINE
     530             :                                       Cdn_atm_rdg, hfreebd,      &   ! LCOV_EXCL_LINE
     531             :                                       hdraft,   hridge,          &   ! LCOV_EXCL_LINE
     532             :                                       distrdg,  hkeel,           &   ! LCOV_EXCL_LINE
     533             :                                       dkeel,    lfloe,           &   ! LCOV_EXCL_LINE
     534             :                                       dfloe,    ncat)
     535             : 
     536             :       use icepack_tracers, only: tr_pond, tr_pond_lvl, tr_pond_topo
     537             : 
     538             :       integer (kind=int_kind), intent(in) :: &
     539             :          ncat
     540             : 
     541             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     542             :          apnd     ,& ! melt pond fraction of sea ice    ! LCOV_EXCL_LINE
     543             :          hpnd     ,& ! mean melt pond depth over sea ice    ! LCOV_EXCL_LINE
     544             :          ipnd     ,& ! mean ice pond depth over sea ice in cat n   ! LCOV_EXCL_LINE
     545             :          alvl     ,& ! level ice area fraction (of grid cell ?)   ! LCOV_EXCL_LINE
     546             :          vlvl        ! level ice mean thickness 
     547             :          
     548             :       real (kind=dbl_kind), intent(in) :: &
     549             :          aice     , & ! concentration of ice   ! LCOV_EXCL_LINE
     550             :          vice     , & ! volume per unit area of ice   ! LCOV_EXCL_LINE
     551             :          vsno         ! volume per unit area of snow 
     552             :          
     553             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     554             :          aicen    , & ! concentration of ice   ! LCOV_EXCL_LINE
     555             :          vicen        ! volume per unit area of ice (m)
     556             :      
     557             :       real (kind=dbl_kind), &
     558             :          intent(out) :: &   ! LCOV_EXCL_LINE
     559             :          hfreebd      , & ! freeboard (m)   ! LCOV_EXCL_LINE
     560             :          hdraft       , & ! draught of ice + snow column (Stoessel1993)   ! LCOV_EXCL_LINE
     561             :          hridge       , & ! ridge height   ! LCOV_EXCL_LINE
     562             :          distrdg      , & ! distance between ridges   ! LCOV_EXCL_LINE
     563             :          hkeel        , & ! keel depth   ! LCOV_EXCL_LINE
     564             :          dkeel        , & ! distance between keels   ! LCOV_EXCL_LINE
     565             :          lfloe        , & ! floe length (m)   ! LCOV_EXCL_LINE
     566             :          dfloe        , & ! distance between floes   ! LCOV_EXCL_LINE
     567             :          Cdn_ocn      , & ! ocean-ice neutral drag coefficient    ! LCOV_EXCL_LINE
     568             :          Cdn_ocn_skin , & ! drag coefficient due to skin drag    ! LCOV_EXCL_LINE
     569             :          Cdn_ocn_floe , & ! drag coefficient due to floe edges    ! LCOV_EXCL_LINE
     570             :          Cdn_ocn_keel , & ! drag coefficient due to keels    ! LCOV_EXCL_LINE
     571             :          Cdn_atm      , & ! ice-atmosphere drag coefficient    ! LCOV_EXCL_LINE
     572             :          Cdn_atm_skin , & ! drag coefficient due to skin drag    ! LCOV_EXCL_LINE
     573             :          Cdn_atm_floe , & ! drag coefficient due to floe edges    ! LCOV_EXCL_LINE
     574             :          Cdn_atm_pond , & ! drag coefficient due to ponds    ! LCOV_EXCL_LINE
     575             :          Cdn_atm_rdg      ! drag coefficient due to ridges 
     576             : 
     577             :       real (kind=dbl_kind), parameter :: & 
     578             :                                       ! [,] = range of values that can be tested 
     579             :          csw       = 0.002_dbl_kind ,&! ice-ocn drag coefficient [0.0005,0.005]
     580             :          csa       = 0.0005_dbl_kind,&! ice-air drag coefficient [0.0001,0.001]    ! LCOV_EXCL_LINE
     581             :          mrdg      = c20            ,&! screening effect see Lu2011 [5,50]   ! LCOV_EXCL_LINE
     582             :          mrdgo     = c10            ,&! screening effect see Lu2011 [5,50]   ! LCOV_EXCL_LINE
     583             :          beta      = p5             ,&! power exponent appearing in astar and    ! LCOV_EXCL_LINE
     584             :                                       ! L=Lmin(A*/(A*-A))**beta [0,1]
     585             :          Lmin      = c8             ,&! min length of floe (m) [5,100]
     586             :          Lmax      = 300._dbl_kind  ,&! max length of floe (m) [30,3000]   ! LCOV_EXCL_LINE
     587             :          cfa       = p2             ,&! Eq. 12 ratio of local from drag over    ! LCOV_EXCL_LINE
     588             :                                       ! geometrical parameter [0,1] 
     589             :          cfw       = p2             ,&! Eq. 15 ratio of local from drag over 
     590             :                                       ! geometrical parameter [0,1]
     591             :          cpa       = p2             ,&! Eq. 16 ratio of local form drag over 
     592             :                                       ! geometrical parameter [0,1]
     593             :          cra       = p2             ,&! Eq. 10 local form drag coefficient [0,1]
     594             :          crw       = p2             ,&! Eq. 11 local form drag coefficient [0,1]   ! LCOV_EXCL_LINE
     595             :          sl        = 22._dbl_kind   ,&! Sheltering parameter Lupkes2012 [10,30]   ! LCOV_EXCL_LINE
     596             :          lpmin     = 2.26_dbl_kind  ,&! min pond length (m) see Eq. 17 [1,10]   ! LCOV_EXCL_LINE
     597             :          lpmax     = 24.63_dbl_kind ,&! max pond length (m) see Eq. 17 [10,100]   ! LCOV_EXCL_LINE
     598             :          tanar     = p4             ,&! 0.25 sail slope = 14 deg [0.4,1]   ! LCOV_EXCL_LINE
     599             :          tanak     = p4             ,&! 0.58 keel slope = 30 deg [0.4,1]   ! LCOV_EXCL_LINE
     600             :          phir      = 0.8_dbl_kind   ,&! porosity of ridges [0.4,1]   ! LCOV_EXCL_LINE
     601             :          phik      = 0.8_dbl_kind   ,&! porosity of keels  [0.4,1]   ! LCOV_EXCL_LINE
     602             :          hkoverhr  = c4             ,&! hkeel/hridge ratio [4,8]   ! LCOV_EXCL_LINE
     603             :          dkoverdr  = c1             ,&! dkeel/distrdg ratio [1,5]   ! LCOV_EXCL_LINE
     604             :          sHGB      = 0.18_dbl_kind  ,&! Lupkes2012 Eq. 28, Hanssen1988,    ! LCOV_EXCL_LINE
     605             :                                       ! Steele1989 suggest instead 0.18
     606             :          alpha2    = c0             ,&! weight functions for area of 
     607             :          beta2     = p75              ! ridged ice [0,1]
     608             : 
     609             :        integer (kind=int_kind) :: &
     610             :          n            ! category index       
     611             : 
     612             :       real (kind=dbl_kind) :: &
     613             :          astar,     & ! new constant for form drag   ! LCOV_EXCL_LINE
     614             :          ctecaf,    & ! constante   ! LCOV_EXCL_LINE
     615             :          ctecwf,    & ! constante   ! LCOV_EXCL_LINE
     616             :          sca,       & ! wind attenuation function   ! LCOV_EXCL_LINE
     617             :          scw,       & ! ocean attenuation function       ! LCOV_EXCL_LINE
     618             :          lp,        & ! pond length (m)   ! LCOV_EXCL_LINE
     619             :          ctecar,    &   ! LCOV_EXCL_LINE
     620             :          ctecwk,    &   ! LCOV_EXCL_LINE
     621             :          ai, aii,   & ! ice area and its inverse   ! LCOV_EXCL_LINE
     622             :          ocnrufi,   & ! inverse ocean roughness   ! LCOV_EXCL_LINE
     623             :          icerufi,   & ! inverse ice roughness   ! LCOV_EXCL_LINE
     624      192144 :          tmp1         ! temporary
     625             : 
     626             :       real (kind=dbl_kind) :: &
     627             :          apond    , & ! melt pond fraction of grid cell   ! LCOV_EXCL_LINE
     628             :          ardg     , & ! ridged ice area fraction of grid cell   ! LCOV_EXCL_LINE
     629      192144 :          vrdg         ! ridged ice mean thickness  
     630             : 
     631             :       real (kind=dbl_kind), parameter :: &
     632             :          ocnruf   = 0.000327_dbl_kind ! ocean surface roughness (m)
     633             : 
     634             :       real (kind=dbl_kind), parameter :: &
     635             :          camax    = 0.02_dbl_kind , & ! Maximum for atmospheric drag   ! LCOV_EXCL_LINE
     636             :          cwmax    = 0.06_dbl_kind     ! Maximum for ocean drag
     637             : 
     638             :       character(len=*),parameter :: subname='(neutral_drag_coeffs)'
     639             : 
     640    20367264 :       astar = c1/(c1-(Lmin/Lmax)**(c1/beta))
     641             : 
     642             :       !-----------------------------------------------------------------
     643             :       ! Initialize across entire grid
     644             :       !-----------------------------------------------------------------
     645             : 
     646    20367264 :       ocnrufi  = c1/ocnruf    ! inverse ocean roughness
     647    20367264 :       icerufi  = c1/iceruf    ! inverse ice roughness
     648    20367264 :       hfreebd=c0
     649    20367264 :       hdraft =c0       
     650    20367264 :       hridge =c0       
     651    20367264 :       distrdg=c0    
     652    20367264 :       hkeel  =c0 
     653    20367264 :       dkeel  =c0 
     654    20367264 :       lfloe  =c0
     655    20367264 :       dfloe  =c0 
     656    20367264 :       Cdn_ocn=dragio
     657    20367264 :       Cdn_ocn_skin=c0 
     658    20367264 :       Cdn_ocn_floe=c0 
     659    20367264 :       Cdn_ocn_keel=c0 
     660    20367264 :       Cdn_atm = (vonkar/log(zref/iceruf)) * (vonkar/log(zref/iceruf))
     661    20367264 :       Cdn_atm_skin=c0 
     662    20367264 :       Cdn_atm_floe=c0 
     663    20367264 :       Cdn_atm_pond=c0 
     664    20367264 :       Cdn_atm_rdg =c0
     665             : 
     666    20367264 :       if (aice > p001) then 
     667             :       
     668     4625699 :          Cdn_atm_skin = csa
     669     4625699 :          Cdn_ocn_skin = csw
     670             : 
     671     4625699 :          ai  = aice
     672     4625699 :          aii = c1/ai
     673             :          
     674             :       !------------------------------------------------------------
     675             :       ! Compute average quantities
     676             :       !------------------------------------------------------------
     677             : 
     678             :          ! ponds
     679     4625699 :          apond = c0
     680     4625699 :          if (tr_pond) then
     681    27754194 :             do n = 1,ncat
     682             :                ! area of pond per unit area of grid cell 
     683    27754194 :                apond = apond+apnd(n)*aicen(n)  
     684             :             enddo
     685             :          endif
     686             :          
     687             :          ! draft and freeboard (see Eq. 27)
     688     4625699 :          hdraft = (rhoi*vice+rhos*vsno)*aii/rhow ! without ponds
     689     4625699 :          hfreebd = (vice+vsno)*aii-hdraft
     690             :          
     691             :          ! Do not allow draft larger than ice thickness (see Eq. 28)  
     692     4625699 :          if (hdraft >= vice*aii) then
     693             :             ! replace excess snow with ice so hi~=hdraft
     694             :             hfreebd = (hdraft*ai*(c1-rhoi/rhow) + &
     695             :                       (vsno-(vice-hdraft*ai)*rhoi/rhos) * &   ! LCOV_EXCL_LINE
     696         486 :                       (c1-rhos/rhow))*aii ! Stoessel1993  
     697             :          endif
     698             :          
     699             :          ! floe size parameterization see Eq. 13
     700     4625699 :          lfloe = Lmin * (astar / (astar - ai))**beta
     701             :          
     702             :          ! distance between floes parameterization see Eq. 14
     703     4625699 :          dfloe = lfloe * (c1/sqrt(ai) - c1)
     704             :          
     705             :          ! Relate ridge height and distance between ridges to 
     706             :          ! ridged ice area fraction and ridged ice mean thickness
     707             :          ! Assumes total volume of ridged ice is split into ridges and keels.
     708             :          ! Then assume total ridges volume over total area of ridges = 
     709             :          ! volume of one average ridge / area of one average ridge
     710             :          ! Same for keels.
     711             :          
     712     4625699 :          ardg=c0
     713     4625699 :          vrdg=c0
     714    27754194 :          do n=1,ncat
     715             :             ! ridged ice area fraction over grid cell 
     716    23128495 :             ardg=ardg+(c1-alvl(n))*aicen(n)
     717             :             ! total ridged ice volume per unit grid cell area
     718    27754194 :             vrdg=vrdg+(c1-vlvl(n))*vicen(n)
     719             :          enddo
     720             :          
     721             :          ! hridge, hkeel, distrdg and dkeel estimates from CICE for 
     722             :          ! simple triangular geometry
     723     4625699 :          if (ardg > p001) then 
     724             :             ! see Eq. 25 and Eq. 26
     725             :             hridge = vrdg/ardg*c2 &
     726             :                    * (alpha2+beta2*hkoverhr/dkoverdr*tanar/tanak) &   ! LCOV_EXCL_LINE
     727      683592 :                    / (phir*c1+phik*tanar/tanak*hkoverhr**c2/dkoverdr)
     728             :             distrdg = c2*hridge*ai/ardg &
     729      683592 :                     * (alpha2/tanar+beta2/tanak*hkoverhr/dkoverdr)
     730      683592 :             hkeel = hkoverhr * hridge
     731      683592 :             dkeel = dkoverdr * distrdg
     732             :             
     733             :           ! Use the height of ridges relative to the mean freeboard of
     734             :           ! the pack.  Therefore skin drag and ridge drag differ in
     735             :           ! this code as compared to  Tsamados et al. (2014) equations
     736             :           ! 10 and 18, which reference both to sea level. 
     737      683592 :           tmp1 = max(c0,hridge - hfreebd)
     738             : 
     739             :       !------------------------------------------------------------
     740             :       ! Skin drag (atmo)
     741             :       !------------------------------------------------------------ 
     742             : 
     743      683592 :           Cdn_atm_skin = csa*(c1 - mrdg*tmp1/distrdg)
     744      683592 :           Cdn_atm_skin = max(min(Cdn_atm_skin,camax),c0)
     745             : 
     746             :       !------------------------------------------------------------
     747             :       ! Ridge effect (atmo)
     748             :       !------------------------------------------------------------
     749             : 
     750      683592 :           if (tmp1 > puny) then
     751      683592 :             sca = c1 - exp(-sHGB*distrdg/tmp1) ! see Eq. 9
     752      683592 :             ctecar = cra*p5
     753             :             Cdn_atm_rdg = ctecar*tmp1/distrdg*sca* &
     754      683592 :                        (log(tmp1*icerufi)/log(zref*icerufi))**c2
     755      683592 :             Cdn_atm_rdg = min(Cdn_atm_rdg,camax)
     756             :           endif
     757             : 
     758             :           ! Use the depth of keels relative to the mean draft of
     759             :           ! the pack.  Therefore skin drag and keel drag differ in
     760             :           ! this code as compared to  Tsamados et al. (2014) equations
     761             :           ! 11 and 19, which reference both to  sea level. In some
     762             :           ! circumstances, hkeel can be less than hdraft because hkoverhr
     763             :           ! is constant, and max(c0,...) temporarily addresses this.
     764      683592 :           tmp1 = max(c0,hkeel - hdraft)
     765             : 
     766             :       !------------------------------------------------------------
     767             :       ! Skin drag bottom ice (ocean)
     768             :       !------------------------------------------------------------ 
     769             :   
     770      683592 :           Cdn_ocn_skin = csw * (c1 - mrdgo*tmp1/dkeel)
     771      683592 :           Cdn_ocn_skin = max(min(Cdn_ocn_skin,cwmax), c0)
     772             :   
     773             :       !------------------------------------------------------------
     774             :       ! Keel effect (ocean)
     775             :       !------------------------------------------------------------
     776             : 
     777      683592 :           if (tmp1 > puny) then
     778      683592 :             scw = c1 - exp(-sHGB*dkeel/tmp1) 
     779      683592 :             ctecwk = crw*p5
     780             :             Cdn_ocn_keel = ctecwk*tmp1/dkeel*scw* &
     781      683592 :                         (log(tmp1*icerufi)/log(zref*icerufi))**c2  
     782      683592 :             Cdn_ocn_keel = max(min(Cdn_ocn_keel,cwmax),c0)
     783             :           endif
     784             :   
     785             :          endif ! ardg > 0.001
     786             : 
     787             :       !------------------------------------------------------------
     788             :       ! Floe edge drag effect (atmo)
     789             :       !------------------------------------------------------------
     790             : 
     791     4625699 :         if (hfreebd > puny) then
     792     4625699 :           sca = c1 - exp(-sl*beta*(c1-ai))
     793     4625699 :           ctecaf = cfa*p5*(log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2*sca
     794     4625699 :           Cdn_atm_floe = ctecaf * hfreebd / lfloe
     795     4625699 :           Cdn_atm_floe = max(min(Cdn_atm_floe,camax),c0)
     796             :         endif
     797             : 
     798             :       !------------------------------------------------------------
     799             :       ! Pond edge effect (atmo)
     800             :       !------------------------------------------------------------
     801             : 
     802     4625699 :         if (hfreebd > puny) then
     803     4625699 :           sca = (apond)**(c1/(zref*beta))
     804     4625699 :           lp  = lpmin*(1-apond)+lpmax*apond
     805             :           Cdn_atm_pond = cpa*p5*sca*apond*hfreebd/lp &
     806     4625699 :                    * (log(hfreebd*ocnrufi)/log(zref*ocnrufi))**c2
     807     4625699 :           Cdn_atm_pond = min(Cdn_atm_pond,camax)
     808             :         endif
     809             :   
     810             :       !------------------------------------------------------------
     811             :       ! Floe edge drag effect (ocean)
     812             :       !------------------------------------------------------------
     813             : 
     814     4625699 :         if (hdraft > puny) then
     815     4625699 :           scw = c1 - exp(-sl*beta*(c1-ai))
     816     4625699 :           ctecwf = cfw*p5*(log(hdraft*ocnrufi)/log(zref*ocnrufi))**c2*scw
     817     4625699 :           Cdn_ocn_floe = ctecwf * hdraft / lfloe
     818     4625699 :           Cdn_ocn_floe = max(min(Cdn_ocn_floe,cwmax),c0)
     819             :         endif
     820             : 
     821             :       !------------------------------------------------------------
     822             :       ! Total drag coefficient (atmo)
     823             :       !------------------------------------------------------------
     824             : 
     825     4625699 :          Cdn_atm = Cdn_atm_skin + Cdn_atm_floe + Cdn_atm_pond + Cdn_atm_rdg
     826     4625699 :          Cdn_atm = min(Cdn_atm,camax)
     827             : 
     828             :       !------------------------------------------------------------
     829             :       ! Total drag coefficient (ocean)
     830             :       !------------------------------------------------------------
     831             : 
     832     4625699 :          Cdn_ocn = Cdn_ocn_skin + Cdn_ocn_floe + Cdn_ocn_keel
     833     4625699 :          Cdn_ocn = min(Cdn_ocn,cwmax)  
     834             : 
     835             :       endif
     836             : 
     837    20367264 :       end subroutine neutral_drag_coeffs
     838             : 
     839             : !=======================================================================
     840             : !autodocument_start icepack_atm_boundary
     841             : ! 
     842             : 
     843  1923648360 :       subroutine icepack_atm_boundary(sfctype,                   &
     844             :                                      Tsf,         potT,          &   ! LCOV_EXCL_LINE
     845             :                                      uatm,        vatm,          &   ! LCOV_EXCL_LINE
     846             :                                      wind,        zlvl,          &   ! LCOV_EXCL_LINE
     847             :                                      Qa,          rhoa,          &   ! LCOV_EXCL_LINE
     848             :                                      strx,        stry,          &   ! LCOV_EXCL_LINE
     849             :                                      Tref,        Qref,          &   ! LCOV_EXCL_LINE
     850             :                                      delt,        delq,          &   ! LCOV_EXCL_LINE
     851             :                                      lhcoef,      shcoef,        &   ! LCOV_EXCL_LINE
     852             :                                      Cdn_atm,                    &   ! LCOV_EXCL_LINE
     853             :                                      Cdn_atm_ratio_n,            &   ! LCOV_EXCL_LINE
     854             :                                      Qa_iso,      Qref_iso,      &   ! LCOV_EXCL_LINE
     855             :                                      uvel,        vvel,          &   ! LCOV_EXCL_LINE
     856             :                                      Uref)
     857             : 
     858             :       character (len=3), intent(in) :: &
     859             :          sfctype      ! ice or ocean
     860             : 
     861             :       real (kind=dbl_kind), intent(in) :: &
     862             :          Tsf      , & ! surface temperature of ice or ocean   ! LCOV_EXCL_LINE
     863             :          potT     , & ! air potential temperature  (K)   ! LCOV_EXCL_LINE
     864             :          uatm     , & ! x-direction wind speed (m/s)   ! LCOV_EXCL_LINE
     865             :          vatm     , & ! y-direction wind speed (m/s)   ! LCOV_EXCL_LINE
     866             :          wind     , & ! wind speed (m/s)   ! LCOV_EXCL_LINE
     867             :          zlvl     , & ! atm level height (m)   ! LCOV_EXCL_LINE
     868             :          Qa       , & ! specific humidity (kg/kg)   ! LCOV_EXCL_LINE
     869             :          rhoa         ! air density (kg/m^3)
     870             : 
     871             :       real (kind=dbl_kind), intent(inout) :: &
     872             :          Cdn_atm  , &    ! neutral drag coefficient   ! LCOV_EXCL_LINE
     873             :          Cdn_atm_ratio_n ! ratio drag coeff / neutral drag coeff
     874             : 
     875             :       real (kind=dbl_kind), &
     876             :          intent(inout) :: &   ! LCOV_EXCL_LINE
     877             :          strx     , & ! x surface stress (N)   ! LCOV_EXCL_LINE
     878             :          stry         ! y surface stress (N)
     879             : 
     880             :       real (kind=dbl_kind), intent(inout) :: &
     881             :          Tref     , & ! reference height temperature  (K)   ! LCOV_EXCL_LINE
     882             :          Qref     , & ! reference height specific humidity (kg/kg)   ! LCOV_EXCL_LINE
     883             :          delt     , & ! potential T difference   (K)   ! LCOV_EXCL_LINE
     884             :          delq     , & ! humidity difference      (kg/kg)   ! LCOV_EXCL_LINE
     885             :          shcoef   , & ! transfer coefficient for sensible heat   ! LCOV_EXCL_LINE
     886             :          lhcoef       ! transfer coefficient for latent heat
     887             : 
     888             :       real (kind=dbl_kind), intent(in), optional, dimension(:) :: &
     889             :          Qa_iso       ! specific isotopic humidity (kg/kg)
     890             : 
     891             :       real (kind=dbl_kind), intent(inout), optional, dimension(:) :: &
     892             :          Qref_iso     ! reference specific isotopic humidity (kg/kg)
     893             : 
     894             :       real (kind=dbl_kind), optional, intent(in) :: &
     895             :          uvel     , & ! x-direction ice speed (m/s)   ! LCOV_EXCL_LINE
     896             :          vvel         ! y-direction ice speed (m/s)
     897             : 
     898             :       real (kind=dbl_kind), optional, intent(out) :: &
     899             :          Uref         ! reference height wind speed (m/s)
     900             : 
     901             : !autodocument_end
     902             : 
     903             :       ! local variables
     904             : 
     905             :       real (kind=dbl_kind) :: &
     906   216108905 :          l_uvel, l_vvel, l_Uref
     907             : 
     908             :       real (kind=dbl_kind), dimension(:), allocatable :: &
     909  2139757265 :          l_Qa_iso, l_Qref_iso   ! local copies of Qa_iso, Qref_iso
     910             : 
     911             :       logical (kind=log_kind) :: &
     912             :          iso_flag  ! flag to turn on iso calcs in other subroutines
     913             : 
     914             :       character(len=*),parameter :: subname='(icepack_atm_boundary)'
     915             : 
     916  2139757265 :       l_uvel = c0
     917  2139757265 :       l_vvel = c0
     918  2139757265 :       l_Uref = c0
     919  2139757265 :       if (present(uvel)) then
     920   975413849 :          l_uvel = uvel
     921             :       endif
     922  2139757265 :       if (present(vvel)) then
     923   975413849 :          l_vvel = vvel
     924             :       endif
     925  2139757265 :       if (present(Qa_iso) .and. present(Qref_iso)) then
     926   975413849 :          iso_flag = .true.
     927   975413849 :          allocate(l_Qa_iso(size(Qa_iso,dim=1)))
     928   975413849 :          allocate(l_Qref_iso(size(Qref_iso,dim=1)))
     929  3901655396 :          l_Qa_iso = Qa_iso
     930  1041260231 :          l_Qref_iso = Qref_iso
     931             :       else
     932  1164343416 :          iso_flag = .false.
     933  1164343416 :          allocate(l_Qa_iso(1))
     934  1164343416 :          allocate(l_Qref_iso(1))
     935  2328686832 :          l_Qa_iso = c0
     936  2328686832 :          l_Qref_iso = c0
     937             :       endif
     938             : 
     939  2139757265 :       Cdn_atm_ratio_n = c1
     940             : 
     941  2139757265 :       if (trim(atmbndy) == 'constant') then
     942             :          call atmo_boundary_const (sfctype,  calc_strair, &
     943             :                                    uatm,     vatm,     &   ! LCOV_EXCL_LINE
     944             :                                    wind,     rhoa,     &   ! LCOV_EXCL_LINE
     945             :                                    strx,     stry,     &   ! LCOV_EXCL_LINE
     946             :                                    Tsf,      potT,     &   ! LCOV_EXCL_LINE
     947             :                                    Qa,                 &   ! LCOV_EXCL_LINE
     948             :                                    delt,     delq,     &   ! LCOV_EXCL_LINE
     949     2079360 :                                    lhcoef,   shcoef    )
     950     2079360 :          if (icepack_warnings_aborted(subname)) return
     951             :       else ! default
     952             :          call atmo_boundary_layer (sfctype,                 &
     953             :                                    calc_strair, formdrag,   &   ! LCOV_EXCL_LINE
     954             :                                    Tsf,      potT,          &   ! LCOV_EXCL_LINE
     955             :                                    uatm,     vatm,          &   ! LCOV_EXCL_LINE
     956             :                                    wind,     zlvl,          &   ! LCOV_EXCL_LINE
     957             :                                    Qa,       rhoa,          &   ! LCOV_EXCL_LINE
     958             :                                    strx,     stry,          &   ! LCOV_EXCL_LINE
     959             :                                    Tref,     Qref,          &   ! LCOV_EXCL_LINE
     960             :                                    delt,     delq,          &   ! LCOV_EXCL_LINE
     961             :                                    lhcoef,   shcoef,        &   ! LCOV_EXCL_LINE
     962             :                                    Cdn_atm,                 &   ! LCOV_EXCL_LINE
     963             :                                    Cdn_atm_ratio_n,         &   ! LCOV_EXCL_LINE
     964             :                                    iso_flag = iso_flag,     &   ! LCOV_EXCL_LINE
     965             :                                    Qa_iso=l_Qa_iso,         &   ! LCOV_EXCL_LINE
     966             :                                    Qref_iso=l_Qref_iso,     &   ! LCOV_EXCL_LINE
     967             :                                    uvel=l_uvel, vvel=l_vvel,  &   ! LCOV_EXCL_LINE
     968  2137677905 :                                    Uref=l_Uref)
     969  2137677905 :          if (icepack_warnings_aborted(subname)) return
     970             :       endif ! atmbndy
     971             : 
     972  2139757265 :       if (present(Uref)) then
     973   975413849 :          Uref = l_Uref
     974             :       endif
     975             : 
     976  2139757265 :       if (present(Qref_iso)) then
     977  1041260231 :          Qref_iso = l_Qref_iso
     978             :       endif
     979             : 
     980  2139757265 :       deallocate(l_Qa_iso,l_Qref_iso)
     981             : 
     982  2139757265 :       end subroutine icepack_atm_boundary
     983             : 
     984             : !------------------------------------------------------------
     985             : ! Define functions
     986             : !------------------------------------------------------------
     987             : 
     988             : !=======================================================================
     989             : 
     990 10912790654 :       real(kind=dbl_kind) function psimhu(xd)
     991             : 
     992             :       real(kind=dbl_kind), intent(in) :: xd
     993             : 
     994             :       psimhu = log((c1+xd*(c2+xd))*(c1+xd*xd)/c8) &
     995 10912790654 :              - c2*atan(xd) + pih
     996             : !ech         - c2*atan(xd) + 1.571_dbl_kind
     997             : 
     998 10912790654 :       end function psimhu
     999             : 
    1000             : !=======================================================================
    1001             : 
    1002 13050468559 :       real(kind=dbl_kind) function psixhu(xd)
    1003             : 
    1004             :       real(kind=dbl_kind), intent(in) :: xd
    1005             : 
    1006 13050468559 :       psixhu =  c2 * log((c1 + xd*xd)/c2)
    1007             : 
    1008 13050468559 :       end function psixhu
    1009             : 
    1010             : !=======================================================================
    1011             : 
    1012             :       end module icepack_atmo
    1013             : 
    1014             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd