LCOV - code coverage report
Current view: top level - columnphysics - icepack_therm_shared.F90 (source / functions) Hit Total Coverage
Test: 200704-015006:b1e41d9f12:3:base,travis,quick Lines: 87 102 85.29 %
Date: 2020-07-03 20:11:40 Functions: 7 10 70.00 %

          Line data    Source code
       1             : !=========================================================================
       2             : !
       3             : ! Shared thermo variables, subroutines
       4             : !
       5             : ! authors: Elizabeth C. Hunke, LANL
       6             : 
       7             :       module icepack_therm_shared
       8             : 
       9             :       use icepack_kinds
      10             : 
      11             :       use icepack_parameters, only: c0, c1, c2, c4, p5, pi
      12             :       use icepack_parameters, only: cp_ocn, cp_ice, rhoi, rhos, Tffresh, TTTice, qqqice
      13             :       use icepack_parameters, only: stefan_boltzmann, emissivity, Lfresh, Tsmelt
      14             :       use icepack_parameters, only: saltmax, min_salin, depressT
      15             :       use icepack_parameters, only: ktherm, heat_capacity, tfrz_option
      16             :       use icepack_parameters, only: calc_Tsfc
      17             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      18             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      19             : 
      20             :       use icepack_mushy_physics, only: enthalpy_mush
      21             :       use icepack_mushy_physics, only: temperature_snow
      22             :       use icepack_mushy_physics, only: enthalpy_snow
      23             :       use icepack_mushy_physics, only: icepack_mushy_temperature_mush
      24             :       use icepack_mushy_physics, only: liquidus_temperature_mush
      25             :     
      26             :       implicit none
      27             : 
      28             :       private
      29             :       public :: calculate_Tin_from_qin, &
      30             :                 surface_heat_flux, &
      31             :                 dsurface_heat_flux_dTsf, &
      32             :                 icepack_init_thermo, &
      33             :                 icepack_init_trcr, &
      34             :                 icepack_ice_temperature, &
      35             :                 icepack_snow_temperature, &
      36             :                 icepack_liquidus_temperature, &
      37             :                 icepack_sea_freezing_temperature, &
      38             :                 icepack_enthalpy_snow
      39             : 
      40             :       real (kind=dbl_kind), parameter, public :: &
      41             :          ferrmax = 1.0e-3_dbl_kind    ! max allowed energy flux error (W m-2)
      42             :                                       ! recommend ferrmax < 0.01 W m-2
      43             : 
      44             :       real (kind=dbl_kind), parameter, public :: &
      45             :          Tmin = -100.0_dbl_kind ! min allowed internal temperature (deg C)
      46             : 
      47             :       logical (kind=log_kind), public :: &
      48             :          l_brine         ! if true, treat brine pocket effects
      49             : 
      50             :       real (kind=dbl_kind), parameter, public :: &
      51             :          hfrazilmin = 0.05_dbl_kind ! min thickness of new frazil ice (m)
      52             : 
      53             :       real (kind=dbl_kind), public :: &
      54             :          hi_min          ! minimum ice thickness allowed (m)
      55             : 
      56             : !=======================================================================
      57             : 
      58             :       contains
      59             : 
      60             : !=======================================================================
      61             : !
      62             : !  Compute the internal ice temperatures from enthalpy using
      63             : !  quadratic formula
      64             : 
      65     9054323 :       function calculate_Tin_from_qin (qin, Tmltk) &
      66     3232100 :                result(Tin)
      67             : 
      68             :       real (kind=dbl_kind), intent(in) :: &
      69             :          qin   , &              ! enthalpy
      70             :          Tmltk                  ! melting temperature at one level
      71             : 
      72             :       real (kind=dbl_kind) :: &
      73             :          Tin                 ! internal temperature
      74             : 
      75             :       ! local variables
      76             : 
      77             :       real (kind=dbl_kind) :: &
      78     3232100 :          aa1,bb1,cc1         ! quadratic solvers
      79             : 
      80             :       character(len=*),parameter :: subname='(calculate_Tin_from_qin)'
      81             : 
      82     9054323 :       if (l_brine) then
      83     8751050 :          aa1 = cp_ice
      84     8751050 :          bb1 = (cp_ocn-cp_ice)*Tmltk - qin/rhoi - Lfresh 
      85     8751050 :          cc1 = Lfresh * Tmltk
      86             :          Tin =  min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) /  &
      87     8751050 :                          (c2*aa1),Tmltk)
      88             : 
      89             :       else                ! fresh ice
      90      303273 :          Tin = (Lfresh + qin/rhoi) / cp_ice
      91             :       endif
      92             :  
      93     9054323 :       end function calculate_Tin_from_qin
      94             : 
      95             : !=======================================================================
      96             : ! Surface heat flux
      97             : !=======================================================================
      98             : 
      99             : ! heat flux into ice
     100             :     
     101    12748521 :       subroutine surface_heat_flux(Tsf,     fswsfc, &
     102             :                                    rhoa,    flw,    &
     103             :                                    potT,    Qa,     &
     104             :                                    shcoef,  lhcoef, &
     105             :                                    flwoutn, fsensn, &
     106             :                                    flatn,   fsurfn)
     107             : 
     108             :       ! input surface temperature
     109             :       real(kind=dbl_kind), intent(in) :: &
     110             :          Tsf             ! ice/snow surface temperature (C)
     111             :     
     112             :       ! input variables
     113             :       real(kind=dbl_kind), intent(in) :: &
     114             :          fswsfc      , & ! SW absorbed at ice/snow surface (W m-2)
     115             :          rhoa        , & ! air density (kg/m^3)
     116             :          flw         , & ! incoming longwave radiation (W/m^2)
     117             :          potT        , & ! air potential temperature  (K)
     118             :          Qa          , & ! specific humidity (kg/kg)
     119             :          shcoef      , & ! transfer coefficient for sensible heat
     120             :          lhcoef          ! transfer coefficient for latent heat
     121             :     
     122             :       ! output
     123             :       real(kind=dbl_kind), intent(out) :: &
     124             :          fsensn      , & ! surface downward sensible heat (W m-2)
     125             :          flatn       , & ! surface downward latent heat (W m-2)
     126             :          flwoutn     , & ! upward LW at surface (W m-2)
     127             :          fsurfn          ! net flux to top surface, excluding fcondtopn
     128             :     
     129             :       ! local variables
     130             :       real(kind=dbl_kind) :: &
     131     4359646 :          TsfK        , & ! ice/snow surface temperature (K)
     132     4359646 :          Qsfc        , & ! saturated surface specific humidity (kg/kg)
     133     4359646 :          qsat        , & ! the saturation humidity of air (kg/m^3)
     134     4359646 :          flwdabs     , & ! downward longwave absorbed heat flx (W/m^2)
     135     4359646 :          tmpvar          ! 1/TsfK
     136             :     
     137             :       character(len=*),parameter :: subname='(surface_heat_flux)'
     138             : 
     139             :       ! ice surface temperature in Kelvin
     140    12748521 :       TsfK = Tsf + Tffresh
     141             : !      TsfK = max(Tsf + Tffresh, c1)
     142    12748521 :       tmpvar = c1/TsfK
     143             :     
     144             :       ! saturation humidity
     145    12748521 :       qsat    = qqqice * exp(-TTTice*tmpvar)
     146    12748521 :       Qsfc    = qsat / rhoa
     147             :     
     148             :       ! longwave radiative flux
     149    12748521 :       flwdabs =  emissivity * flw
     150    12748521 :       flwoutn = -emissivity * stefan_boltzmann * TsfK**4
     151             :     
     152             :       ! downward latent and sensible heat fluxes
     153    12748521 :       fsensn = shcoef * (potT - TsfK)
     154    12748521 :       flatn  = lhcoef * (Qa - Qsfc)
     155             :     
     156             :       ! combine fluxes
     157    12748521 :       fsurfn = fswsfc + flwdabs + flwoutn + fsensn + flatn
     158             : 
     159    12748521 :       end subroutine surface_heat_flux
     160             : 
     161             :   !=======================================================================
     162             :   
     163     9388321 :       subroutine dsurface_heat_flux_dTsf(Tsf,  rhoa,      &
     164             :                                          shcoef,  lhcoef, &
     165             :                                          dfsurfn_dTsf, dflwoutn_dTsf, &
     166             :                                          dfsensn_dTsf, dflatn_dTsf)
     167             :     
     168             :       ! input surface temperature
     169             :       real(kind=dbl_kind), intent(in) :: &
     170             :          Tsf               ! ice/snow surface temperature (C)
     171             :     
     172             :       ! input variables
     173             :       real(kind=dbl_kind), intent(in) :: &
     174             :          rhoa          , & ! air density (kg/m^3)
     175             :          shcoef        , & ! transfer coefficient for sensible heat
     176             :          lhcoef            ! transfer coefficient for latent heat
     177             :     
     178             :       ! output
     179             :       real(kind=dbl_kind), intent(out) :: &
     180             :          dfsurfn_dTsf      ! derivative of net flux to top surface, excluding fcondtopn
     181             :     
     182             :       real(kind=dbl_kind), intent(out) :: &
     183             :          dflwoutn_dTsf , & ! derivative of longwave flux wrt surface temperature
     184             :          dfsensn_dTsf  , & ! derivative of sensible heat flux wrt surface temperature
     185             :          dflatn_dTsf       ! derivative of latent heat flux wrt surface temperature
     186             :     
     187             :       ! local variables
     188             :       real(kind=dbl_kind) :: &
     189     3220984 :          TsfK          , & ! ice/snow surface temperature (K)
     190     3220984 :          dQsfc_dTsf    , & ! saturated surface specific humidity (kg/kg)
     191     3220984 :          qsat          , & ! the saturation humidity of air (kg/m^3)
     192     3220984 :          tmpvar            ! 1/TsfK
     193             :     
     194             :       character(len=*),parameter :: subname='(dsurface_heat_flux_dTsf)'
     195             : 
     196             :       ! ice surface temperature in Kelvin
     197             : !      TsfK = max(Tsf + Tffresh, c1)
     198     9388321 :       TsfK = Tsf + Tffresh
     199     9388321 :       tmpvar = c1/TsfK
     200             :     
     201             :       ! saturation humidity
     202     9388321 :       qsat          = qqqice * exp(-TTTice*tmpvar)
     203     9388321 :       dQsfc_dTsf    = TTTice * tmpvar * tmpvar * (qsat / rhoa)
     204             :     
     205             :       ! longwave radiative flux
     206     9388321 :       dflwoutn_dTsf = -emissivity * stefan_boltzmann * c4*TsfK**3
     207             :     
     208             :       ! downward latent and sensible heat fluxes
     209     9388321 :       dfsensn_dTsf = -shcoef
     210     9388321 :       dflatn_dTsf  = -lhcoef * dQsfc_dTsf
     211             :     
     212             :       ! combine fluxes
     213     9388321 :       dfsurfn_dTsf = dflwoutn_dTsf + dfsensn_dTsf + dflatn_dTsf
     214             :     
     215     9388321 :       end subroutine dsurface_heat_flux_dTsf
     216             : 
     217             : !=======================================================================
     218             : !autodocument_start icepack_init_thermo
     219             : ! Initialize the vertical profile of ice salinity and melting temperature.
     220             : !
     221             : ! authors: C. M. Bitz, UW
     222             : !          William H. Lipscomb, LANL
     223             : 
     224          45 :       subroutine icepack_init_thermo(nilyr, sprofile)
     225             : 
     226             :       integer (kind=int_kind), intent(in) :: &
     227             :          nilyr                            ! number of ice layers
     228             : 
     229             :       real (kind=dbl_kind), dimension(:), intent(out) :: &
     230             :          sprofile                         ! vertical salinity profile
     231             : 
     232             : !autodocument_end
     233             : 
     234             :       real (kind=dbl_kind), parameter :: &
     235             :          nsal    = 0.407_dbl_kind, &
     236             :          msal    = 0.573_dbl_kind
     237             : 
     238             :       integer (kind=int_kind) :: k        ! ice layer index
     239          16 :       real (kind=dbl_kind)    :: zn       ! normalized ice thickness
     240             : 
     241             :       character(len=*),parameter :: subname='(icepack_init_thermo)'
     242             : 
     243             :       !-----------------------------------------------------------------
     244             :       ! Determine l_brine based on saltmax.
     245             :       ! Thermodynamic solver will not converge if l_brine is true and
     246             :       !  saltmax is close to zero.
     247             :       ! Set l_brine to false for zero layer thermodynamics
     248             :       !-----------------------------------------------------------------
     249             : 
     250          45 :       heat_capacity = .true.      
     251          45 :       if (ktherm == 0) heat_capacity = .false. ! 0-layer thermodynamics
     252             : 
     253          45 :       l_brine = .false.
     254          45 :       if (saltmax > min_salin .and. heat_capacity) l_brine = .true.
     255             : 
     256             :       !-----------------------------------------------------------------
     257             :       ! Prescibe vertical profile of salinity and melting temperature.
     258             :       ! Note this profile is only used for BL99 thermodynamics.
     259             :       !-----------------------------------------------------------------
     260             : 
     261          45 :       if (l_brine) then
     262         312 :          do k = 1, nilyr
     263             :             zn = (real(k,kind=dbl_kind)-p5) /  &
     264         273 :                   real(nilyr,kind=dbl_kind)
     265         273 :             sprofile(k)=(saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
     266         312 :             sprofile(k) = max(sprofile(k), min_salin)
     267             :          enddo ! k
     268          39 :          sprofile(nilyr+1) = saltmax
     269             : 
     270             :       else ! .not. l_brine
     271          18 :          do k = 1, nilyr+1
     272          18 :             sprofile(k) = c0
     273             :          enddo
     274             :       endif ! l_brine
     275             : 
     276          45 :       end subroutine icepack_init_thermo
     277             : 
     278             : !=======================================================================
     279             : !autodocument_start icepack_init_trcr
     280             : !
     281         426 :       subroutine icepack_init_trcr(Tair,     Tf,       &
     282         426 :                                   Sprofile, Tprofile, &
     283             :                                   Tsfc,               &
     284             :                                   nilyr,    nslyr,    &
     285         852 :                                   qin,      qsn)
     286             : 
     287             :       integer (kind=int_kind), intent(in) :: &
     288             :          nilyr, &    ! number of ice layers
     289             :          nslyr       ! number of snow layers
     290             : 
     291             :       real (kind=dbl_kind), intent(in) :: &
     292             :          Tair, &     ! air temperature (C)
     293             :          Tf          ! freezing temperature (C)
     294             : 
     295             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
     296             :          Sprofile, & ! vertical salinity profile (ppt)
     297             :          Tprofile    ! vertical temperature profile (C)
     298             : 
     299             :       real (kind=dbl_kind), intent(out) :: &
     300             :          Tsfc        ! surface temperature (C)
     301             : 
     302             :       real (kind=dbl_kind), dimension(:), intent(out) :: &
     303             :          qin, &      ! ice enthalpy profile (J/m3)
     304             :          qsn         ! snow enthalpy profile (J/m3)
     305             : 
     306             : !autodocument_end
     307             : 
     308             :       ! local variables
     309             : 
     310             :       integer (kind=int_kind) :: k
     311             : 
     312             :       real (kind=dbl_kind) :: &
     313         152 :          slope, Ti
     314             : 
     315             :       character(len=*),parameter :: subname='(icepack_init_trcr)'
     316             : 
     317             :       ! surface temperature
     318         426 :       Tsfc = Tf ! default
     319         426 :       if (calc_Tsfc) Tsfc = min(Tsmelt, Tair - Tffresh) ! deg C
     320             :       
     321         426 :       if (heat_capacity) then
     322             :         
     323             :         ! ice enthalpy
     324        3120 :         do k = 1, nilyr
     325             :           ! assume linear temp profile and compute enthalpy
     326        2730 :           slope = Tf - Tsfc
     327             :           Ti = Tsfc + slope*(real(k,kind=dbl_kind)-p5) &
     328        2730 :               /real(nilyr,kind=dbl_kind)
     329        3120 :           if (ktherm == 2) then
     330        2100 :             qin(k) = enthalpy_mush(Ti, Sprofile(k))
     331             :           else
     332         210 :             qin(k) = -(rhoi * (cp_ice*(Tprofile(k)-Ti) &
     333         630 :                 + Lfresh*(c1-Tprofile(k)/Ti) - cp_ocn*Tprofile(k)))
     334             :           endif
     335             :         enddo               ! nilyr
     336             :         
     337             :         ! snow enthalpy
     338         900 :         do k = 1, nslyr
     339         510 :           Ti = min(c0, Tsfc)
     340         900 :           qsn(k) = -rhos*(Lfresh - cp_ice*Ti)
     341             :         enddo               ! nslyr
     342             :         
     343             :       else  ! one layer with zero heat capacity
     344             :         
     345             :         ! ice energy
     346          36 :         qin(1) = -rhoi * Lfresh 
     347             :         
     348             :         ! snow energy
     349          36 :         qsn(1) = -rhos * Lfresh 
     350             :         
     351             :       endif               ! heat_capacity
     352             :       
     353         426 :     end subroutine icepack_init_trcr
     354             : 
     355             : !=======================================================================
     356             : !autodocument_start icepack_liquidus_temperature
     357             : ! compute liquidus temperature
     358             : 
     359      887700 :       function icepack_liquidus_temperature(Sin) result(Tmlt)
     360             : 
     361             :         real(dbl_kind), intent(in) :: Sin
     362             :         real(dbl_kind) :: Tmlt
     363             : 
     364             : !autodocument_end
     365             : 
     366             :         character(len=*),parameter :: subname='(icepack_liquidus_temperature)'
     367             : 
     368      887700 :         if (ktherm == 2) then
     369             : 
     370      887640 :            Tmlt = liquidus_temperature_mush(Sin)
     371             : 
     372             :         else
     373             : 
     374          60 :            Tmlt = -depressT * Sin
     375             : 
     376             :         endif
     377             : 
     378      887700 :       end function icepack_liquidus_temperature
     379             : 
     380             : !=======================================================================
     381             : !autodocument_start icepack_sea_freezing_temperature
     382             : ! compute ocean freezing temperature
     383             : 
     384     1370160 :       function icepack_sea_freezing_temperature(sss) result(Tf)
     385             : 
     386             :         real(dbl_kind), intent(in) :: sss
     387             :         real(dbl_kind) :: Tf
     388             : 
     389             : !autodocument_end
     390             : 
     391             :         character(len=*),parameter :: subname='(icepack_sea_freezing_temperature)'
     392             : 
     393     1370160 :         if (trim(tfrz_option) == 'mushy') then
     394             : 
     395      887520 :            Tf = icepack_liquidus_temperature(sss) ! deg C
     396             :            
     397      482640 :         elseif (trim(tfrz_option) == 'linear_salt') then
     398             : 
     399      386112 :            Tf = -depressT * sss ! deg C
     400             : 
     401             :         else
     402             : 
     403       96528 :            Tf = -1.8_dbl_kind
     404             : 
     405             :         endif
     406             : 
     407     1370160 :       end function icepack_sea_freezing_temperature
     408             : 
     409             : !=======================================================================
     410             : !autodocument_start icepack_ice_temperature
     411             : ! compute ice temperature
     412             : 
     413           0 :       function icepack_ice_temperature(qin, Sin) result(Tin)
     414             : 
     415             :         real(kind=dbl_kind), intent(in) :: qin, Sin
     416             :         real(kind=dbl_kind) :: Tin
     417             : 
     418             : !autodocument_end
     419             : 
     420           0 :         real(kind=dbl_kind) :: Tmlts
     421             : 
     422             :         character(len=*),parameter :: subname='(icepack_ice_temperature)'
     423             : 
     424           0 :         if (ktherm == 2) then
     425             : 
     426           0 :            Tin = icepack_mushy_temperature_mush(qin, Sin)
     427             : 
     428             :         else
     429             : 
     430           0 :            Tmlts = -depressT * Sin
     431           0 :            Tin = calculate_Tin_from_qin(qin,Tmlts)
     432             : 
     433             :         endif
     434             : 
     435           0 :       end function icepack_ice_temperature
     436             : 
     437             : !=======================================================================
     438             : !autodocument_start icepack_snow_temperature
     439             : ! compute snow temperature
     440             : 
     441           0 :       function icepack_snow_temperature(qin) result(Tsn)
     442             : 
     443             :         real(kind=dbl_kind), intent(in) :: qin
     444             :         real(kind=dbl_kind) :: Tsn
     445             : 
     446             : !autodocument_end
     447             : 
     448             :         character(len=*),parameter :: subname='(icepack_snow_temperature)'
     449             : 
     450           0 :         if (ktherm == 2) then
     451             : 
     452           0 :            Tsn = temperature_snow(qin)
     453             : 
     454             :         else
     455             : 
     456           0 :            Tsn = (Lfresh + qin/rhos)/cp_ice
     457             : 
     458             :         endif
     459             : 
     460           0 :       end function icepack_snow_temperature
     461             : 
     462             : !=======================================================================
     463             : !autodocument_start icepack_enthalpy_snow
     464             : ! compute snow enthalpy
     465             : 
     466           0 :       function icepack_enthalpy_snow(zTsn) result(qsn)
     467             : 
     468             :         real(kind=dbl_kind), intent(in) :: zTsn
     469             :         real(kind=dbl_kind) :: qsn
     470             : 
     471             : !autodocument_end
     472             : 
     473             :         character(len=*),parameter :: subname='(icepack_enthalpy_snow)'
     474             : 
     475           0 :         qsn = enthalpy_snow(zTsn)
     476             : 
     477           0 :       end function icepack_enthalpy_snow
     478             : 
     479             : !=======================================================================
     480             : 
     481             :       end module icepack_therm_shared
     482             : 
     483             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd