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

Generated by: LCOV version 1.14-6-g40580cd