LCOV - code coverage report
Current view: top level - icepack/columnphysics - icepack_therm_vertical.F90 (source / functions) Hit Total Coverage
Test: 200617-180449:aec9683041:7:first,base,travis,decomp,reprosum,io,quick Lines: 707 931 75.94 %
Date: 2020-06-17 18:05:09 Functions: 9 9 100.00 %

          Line data    Source code
       1             : !=========================================================================
       2             : !
       3             : ! Update ice and snow internal temperatures and compute
       4             : ! thermodynamic growth rates and atmospheric fluxes.
       5             : !
       6             : ! NOTE: The thermodynamic calculation is split in two for load balancing.
       7             : !       First icepack_therm_vertical computes vertical growth rates and coupler
       8             : !       fluxes.  Then icepack_therm_itd does thermodynamic calculations not
       9             : !       needed for coupling.
      10             : !
      11             : ! authors: William H. Lipscomb, LANL
      12             : !          C. M. Bitz, UW
      13             : !          Elizabeth C. Hunke, LANL
      14             : !
      15             : ! 2003: Vectorized by Clifford Chen (Fujitsu) and William Lipscomb
      16             : ! 2004: Block structure added by William Lipscomb
      17             : ! 2006: Streamlined for efficiency by Elizabeth Hunke
      18             : !       Converted to free source form (F90)
      19             : 
      20             :       module icepack_therm_vertical
      21             : 
      22             :       use icepack_kinds
      23             :       use icepack_parameters, only: c0, c1, p001, p5, puny
      24             :       use icepack_parameters, only: pi, depressT, Lvap, hs_min, cp_ice, min_salin
      25             :       use icepack_parameters, only: cp_ocn, rhow, rhoi, rhos, Lfresh, rhofresh, ice_ref_salinity
      26             :       use icepack_parameters, only: ktherm, heat_capacity, calc_Tsfc
      27             :       use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair
      28             :       use icepack_parameters, only: rfracmin, rfracmax, pndaspect, dpscale, frzpnd
      29             :       use icepack_parameters, only: phi_i_mushy
      30             : 
      31             :       use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso
      32             :       use icepack_tracers, only: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
      33             :       use icepack_tracers, only: n_aero, n_iso
      34             : 
      35             :       use icepack_therm_shared, only: ferrmax, l_brine
      36             :       use icepack_therm_shared, only: calculate_tin_from_qin, Tmin
      37             :       use icepack_therm_shared, only: hi_min
      38             :       use icepack_therm_bl99,   only: temperature_changes
      39             :       use icepack_therm_0layer, only: zerolayer_temperature
      40             :       use icepack_therm_mushy,  only: temperature_changes_salinity
      41             : 
      42             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      43             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      44             : 
      45             :       use icepack_mushy_physics, only: icepack_mushy_temperature_mush
      46             :       use icepack_mushy_physics, only: liquidus_temperature_mush
      47             :       use icepack_mushy_physics, only: enthalpy_mush, enthalpy_of_melting
      48             : 
      49             :       use icepack_aerosol, only: update_aerosol
      50             :       use icepack_isotope, only: update_isotope
      51             :       use icepack_atmo, only: neutral_drag_coeffs, icepack_atm_boundary
      52             :       use icepack_age, only: increment_age
      53             :       use icepack_firstyear, only: update_FYarea
      54             :       use icepack_flux, only: set_sfcflux, merge_fluxes
      55             :       use icepack_meltpond_cesm, only: compute_ponds_cesm
      56             :       use icepack_meltpond_lvl, only: compute_ponds_lvl
      57             :       use icepack_meltpond_topo, only: compute_ponds_topo
      58             : 
      59             :       implicit none
      60             : 
      61             :       private
      62             :       public :: frzmlt_bottom_lateral, &
      63             :                 thermo_vertical, &
      64             :                 icepack_step_therm1
      65             : 
      66             : !=======================================================================
      67             : 
      68             :       contains
      69             : 
      70             : !=======================================================================
      71             : !
      72             : ! Driver for updating ice and snow internal temperatures and
      73             : ! computing thermodynamic growth rates and atmospheric fluxes.
      74             : !
      75             : ! authors: William H. Lipscomb, LANL
      76             : !          C. M. Bitz, UW
      77             : 
      78   774714937 :       subroutine thermo_vertical (nilyr,       nslyr,     &
      79             :                                   dt,          aicen,     &
      80             :                                   vicen,       vsnon,     &
      81   774714937 :                                   Tsf,         zSin,      &
      82   774714937 :                                   zqin,        zqsn,      &
      83             :                                   apond,       hpond,     &
      84             :                                   tr_pond_topo,&
      85             :                                   flw,         potT,      &
      86             :                                   Qa,          rhoa,      &
      87             :                                   fsnow,       fpond,     &
      88             :                                   fbot,        Tbot,      &
      89             :                                   Tsnice,       sss,       &
      90             :                                   lhcoef,      shcoef,    &
      91             :                                   fswsfc,      fswint,    &
      92   774714937 :                                   Sswabs,      Iswabs,    &
      93             :                                   fsurfn,      fcondtopn, &
      94             :                                   fcondbotn,              &
      95             :                                   fsensn,      flatn,     &
      96             :                                   flwoutn,     evapn,     &
      97             :                                   evapsn,      evapin,    &
      98             :                                   freshn,      fsaltn,    &
      99             :                                   fhocnn,      meltt,     &
     100             :                                   melts,       meltb,     &
     101             :                                   congel,      snoice,    &
     102             :                                   mlt_onset,   frz_onset, &
     103             :                                   yday,        dsnow,     &
     104             :                                   prescribed_ice)
     105             : 
     106             :       integer (kind=int_kind), intent(in) :: &
     107             :          nilyr   , & ! number of ice layers
     108             :          nslyr       ! number of snow layers
     109             : 
     110             :       real (kind=dbl_kind), intent(in) :: &
     111             :          dt          ! time step
     112             : 
     113             :       ! ice state variables
     114             :       real (kind=dbl_kind), intent(inout) :: &
     115             :          aicen   , & ! concentration of ice
     116             :          vicen   , & ! volume per unit area of ice          (m)
     117             :          vsnon       ! volume per unit area of snow         (m)
     118             : 
     119             :       ! tracers
     120             :       real (kind=dbl_kind), intent(inout) :: &
     121             :          Tsf     , & ! ice/snow top surface temp, same as Tsfcn (deg C)
     122             :          apond   , & ! melt pond area fraction
     123             :          hpond       ! melt pond depth (m)
     124             : !        iage        ! ice age (s)
     125             : 
     126             :       logical (kind=log_kind), intent(in) :: &
     127             :          tr_pond_topo    ! if .true., use melt pond tracer
     128             : 
     129             :       logical (kind=log_kind), intent(in), optional :: &
     130             :          prescribed_ice  ! if .true., use prescribed ice instead of computed
     131             : 
     132             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     133             :          zqsn    , & ! snow layer enthalpy, zqsn < 0 (J m-3)
     134             :          zqin    , & ! ice layer enthalpy, zqin < 0 (J m-3)
     135             :          zSin        ! internal ice layer salinities
     136             : 
     137             :       ! input from atmosphere
     138             :       real (kind=dbl_kind), &
     139             :          intent(in) :: &
     140             :          flw     , & ! incoming longwave radiation (W/m^2)
     141             :          potT    , & ! air potential temperature  (K) 
     142             :          Qa      , & ! specific humidity (kg/kg) 
     143             :          rhoa    , & ! air density (kg/m^3) 
     144             :          fsnow   , & ! snowfall rate (kg m-2 s-1)
     145             :          shcoef  , & ! transfer coefficient for sensible heat
     146             :          lhcoef      ! transfer coefficient for latent heat
     147             : 
     148             :       real (kind=dbl_kind), &
     149             :          intent(inout) :: &
     150             :          fswsfc  , & ! SW absorbed at ice/snow surface (W m-2)
     151             :          fswint  , & ! SW absorbed in ice interior, below surface (W m-2)
     152             :          fpond       ! fresh water flux to ponds (kg/m^2/s)
     153             : 
     154             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     155             :          Sswabs  , & ! SW radiation absorbed in snow layers (W m-2)
     156             :          Iswabs      ! SW radiation absorbed in ice layers (W m-2)
     157             : 
     158             :       ! input from ocean
     159             :       real (kind=dbl_kind), intent(in) :: &
     160             :          fbot    , & ! ice-ocean heat flux at bottom surface (W/m^2)
     161             :          Tbot    , & ! ice bottom surface temperature (deg C)
     162             :          sss         ! ocean salinity
     163             : 
     164             :       ! coupler fluxes to atmosphere
     165             :       real (kind=dbl_kind), intent(out):: &
     166             :          flwoutn , & ! outgoing longwave radiation (W/m^2) 
     167             :          evapn   , & ! evaporative water flux (kg/m^2/s) 
     168             :          evapsn  , & ! evaporative water flux over snow (kg/m^2/s) 
     169             :          evapin      ! evaporative water flux over ice (kg/m^2/s) 
     170             : 
     171             :       ! Note: these are intent out if calc_Tsfc = T, otherwise intent in
     172             :       real (kind=dbl_kind), intent(inout):: &
     173             :          fsensn   , & ! sensible heat flux (W/m^2) 
     174             :          flatn    , & ! latent heat flux   (W/m^2) 
     175             :          fsurfn   , & ! net flux to top surface, excluding fcondtopn
     176             :          fcondtopn, & ! downward cond flux at top surface (W m-2)
     177             :          fcondbotn    ! downward cond flux at bottom surface (W m-2)
     178             : 
     179             :       ! coupler fluxes to ocean
     180             :       real (kind=dbl_kind), intent(out):: &
     181             :          freshn  , & ! fresh water flux to ocean (kg/m^2/s)
     182             :          fsaltn  , & ! salt flux to ocean (kg/m^2/s)
     183             :          fhocnn      ! net heat flux to ocean (W/m^2) 
     184             : 
     185             :       ! diagnostic fields
     186             :       real (kind=dbl_kind), &
     187             :          intent(inout):: &
     188             :          Tsnice    , & ! snow ice interface temperature (deg C)
     189             :          meltt    , & ! top ice melt             (m/step-->cm/day) 
     190             :          melts    , & ! snow melt                (m/step-->cm/day) 
     191             :          meltb    , & ! basal ice melt           (m/step-->cm/day) 
     192             :          congel   , & ! basal ice growth         (m/step-->cm/day) 
     193             :          snoice   , & ! snow-ice formation       (m/step-->cm/day) 
     194             :          dsnow    , & ! change in snow thickness (m/step-->cm/day) 
     195             :          mlt_onset, & ! day of year that sfc melting begins 
     196             :          frz_onset    ! day of year that freezing begins (congel or frazil) 
     197             : 
     198             :       real (kind=dbl_kind), intent(in) :: &
     199             :          yday         ! day of year
     200             : 
     201             :       ! local variables
     202             : 
     203             :       integer (kind=int_kind) :: &
     204             :          k               ! ice layer index
     205             : 
     206             :       real (kind=dbl_kind) :: &
     207    54309860 :          dhi         , & ! change in ice thickness
     208    54309860 :          dhs             ! change in snow thickness
     209             : 
     210             : ! 2D state variables (thickness, temperature)
     211             : 
     212             :       real (kind=dbl_kind) :: &
     213    54309860 :          hilyr       , & ! ice layer thickness
     214    54309860 :          hslyr       , & ! snow layer thickness
     215    54309860 :          hin         , & ! ice thickness (m)
     216    54309860 :          hsn         , & ! snow thickness (m)
     217    54309860 :          hsn_new     , & ! thickness of new snow (m)
     218    54309860 :          worki       , & ! local work array
     219    54309860 :          works           ! local work array
     220             : 
     221             :       real (kind=dbl_kind), dimension (nilyr) :: &
     222  1876471506 :          zTin            ! internal ice layer temperatures
     223             : 
     224             :       real (kind=dbl_kind), dimension (nslyr) :: &
     225   883334657 :          zTsn            ! internal snow layer temperatures
     226             : 
     227             : ! other 2D flux and energy variables
     228             : 
     229             :       real (kind=dbl_kind) :: &
     230    54309860 :          einit       , & ! initial energy of melting (J m-2)
     231    54309860 :          efinal      , & ! final energy of melting (J m-2)
     232    54309860 :          einter          ! intermediate energy
     233             : 
     234             :       real (kind=dbl_kind) :: &
     235    54309860 :          fadvocn ! advective heat flux to ocean
     236             : 
     237             :       character(len=*),parameter :: subname='(thermo_vertical)'
     238             : 
     239             :       !-----------------------------------------------------------------
     240             :       ! Initialize
     241             :       !-----------------------------------------------------------------
     242             : 
     243   774714937 :       flwoutn = c0
     244   774714937 :       evapn   = c0
     245   774714937 :       evapsn   = c0
     246   774714937 :       evapin   = c0
     247   774714937 :       freshn  = c0
     248   774714937 :       fsaltn  = c0
     249   774714937 :       fhocnn  = c0
     250   774714937 :       fadvocn = c0
     251             : 
     252   774714937 :       meltt   = c0
     253   774714937 :       meltb   = c0
     254   774714937 :       melts   = c0
     255   774714937 :       congel  = c0
     256   774714937 :       snoice  = c0
     257   774714937 :       dsnow   = c0
     258  1549429874 :       zTsn(:) = c0
     259  5517993338 :       zTin(:) = c0
     260             : 
     261   774714937 :       if (calc_Tsfc) then
     262   744838517 :          fsensn  = c0
     263   744838517 :          flatn     = c0
     264   744838517 :          fsurfn    = c0
     265   744838517 :          fcondtopn = c0
     266             :       endif
     267             : 
     268             :       !-----------------------------------------------------------------
     269             :       ! Compute variables needed for vertical thermo calculation
     270             :       !-----------------------------------------------------------------
     271             : 
     272             :       call init_vertical_profile (nilyr,    nslyr,   &
     273             :                                   aicen,             &
     274             :                                   vicen,    vsnon,   &
     275             :                                   hin,      hilyr,   &
     276             :                                   hsn,      hslyr,   &
     277           0 :                                   zqin,     zTin,    &
     278           0 :                                   zqsn,     zTsn,    &
     279           0 :                                   zSin,              &
     280   774714937 :                                   einit )
     281   774714937 :       if (icepack_warnings_aborted(subname)) return
     282             : 
     283             :       ! Save initial ice and snow thickness (for fresh and fsalt)
     284   774714937 :       worki = hin
     285   774714937 :       works = hsn
     286             : 
     287             :       !-----------------------------------------------------------------
     288             :       ! Compute new surface temperature and internal ice and snow
     289             :       !  temperatures.
     290             :       !-----------------------------------------------------------------
     291             : 
     292   774714937 :       if (heat_capacity) then   ! usual case
     293             : 
     294   661427244 :          if (ktherm == 2) then
     295             : 
     296             :             call temperature_changes_salinity(dt,                   & 
     297             :                                               nilyr,     nslyr,     &
     298             :                                               rhoa,      flw,       &
     299             :                                               potT,      Qa,        &
     300             :                                               shcoef,    lhcoef,    &
     301             :                                               fswsfc,    fswint,    &
     302           0 :                                               Sswabs,    Iswabs,    &
     303             :                                               hilyr,     hslyr,     &
     304             :                                               apond,     hpond,     &
     305           0 :                                               zqin,      zTin,      &
     306           0 :                                               zqsn,      zTsn,      &
     307           0 :                                               zSin,                 &
     308             :                                               Tsf,       Tbot,      &
     309             :                                               sss,                  &
     310             :                                               fsensn,    flatn,     &
     311             :                                               flwoutn,   fsurfn,    &
     312             :                                               fcondtopn, fcondbotn,  &
     313   609317819 :                                               fadvocn,   snoice)
     314   609317819 :             if (icepack_warnings_aborted(subname)) return
     315             : 
     316             :          else ! ktherm
     317             : 
     318             :             call temperature_changes(dt,                   &  
     319             :                                      nilyr,     nslyr,     &
     320             :                                      rhoa,      flw,       &
     321             :                                      potT,      Qa,        &
     322             :                                      shcoef,    lhcoef,    &
     323             :                                      fswsfc,    fswint,    &
     324           0 :                                      Sswabs,    Iswabs,    &
     325             :                                      hilyr,     hslyr,     &
     326           0 :                                      zqin,      zTin,      &
     327           0 :                                      zqsn,      zTsn,      &
     328           0 :                                      zSin,                 &
     329             :                                      Tsf,       Tbot,      &
     330             :                                      fsensn,    flatn,     &
     331             :                                      flwoutn,   fsurfn,    &
     332             :                                      fcondtopn, fcondbotn,  &
     333    52109425 :                                      einit                 )
     334    52109425 :             if (icepack_warnings_aborted(subname)) return
     335             : 
     336             :          endif ! ktherm
     337             :             
     338             :       else
     339             : 
     340   113287693 :          if (calc_Tsfc) then       
     341             : 
     342             :             call zerolayer_temperature(nilyr,     nslyr,    &
     343             :                                        rhoa,      flw,      &
     344             :                                        potT,      Qa,       &
     345             :                                        shcoef,    lhcoef,   &
     346             :                                        fswsfc,              &
     347             :                                        hilyr,     hslyr,    &
     348             :                                        Tsf,       Tbot,     &
     349             :                                        fsensn,    flatn,    &
     350             :                                        flwoutn,   fsurfn,   &
     351   113287693 :                                        fcondtopn, fcondbotn  )
     352   113287693 :             if (icepack_warnings_aborted(subname)) return
     353             : 
     354             :          else
     355             : 
     356             :             !------------------------------------------------------------
     357             :             ! Set fcondbot = fcondtop for zero layer thermodynamics
     358             :             ! fcondtop is set in call to set_sfcflux in step_therm1
     359             :             !------------------------------------------------------------
     360             : 
     361           0 :             fcondbotn  = fcondtopn   ! zero layer         
     362             :       
     363             :          endif      ! calc_Tsfc
     364             : 
     365             :       endif         ! heat_capacity
     366             : 
     367             :       ! intermediate energy for error check
     368             :       
     369   774714937 :       einter = c0
     370  1549429874 :       do k = 1, nslyr
     371  1549429874 :          einter = einter + hslyr * zqsn(k)
     372             :       enddo ! k
     373  5517993338 :       do k = 1, nilyr
     374  5517993338 :          einter = einter + hilyr * zqin(k)
     375             :       enddo ! k
     376             : 
     377   774714937 :       Tsnice = c0
     378   774714937 :       if ((hslyr+hilyr) > puny) then
     379   774714937 :          if (hslyr > puny) then
     380   643151396 :             Tsnice = (hslyr*zTsn(nslyr) + hilyr*zTin(1)) / (hslyr+hilyr)
     381             :          else
     382   131563541 :             Tsnice = Tsf
     383             :          endif
     384             :       endif
     385             : 
     386   774714937 :       if (icepack_warnings_aborted(subname)) return
     387             : 
     388             :       !-----------------------------------------------------------------
     389             :       ! Compute growth and/or melting at the top and bottom surfaces.
     390             :       ! Add new snowfall.
     391             :       ! Repartition ice into equal-thickness layers, conserving energy.
     392             :       !----------------------------------------------------------------- 
     393             : 
     394             :       call thickness_changes(nilyr,       nslyr,     &
     395             :                              dt,          yday,      &
     396             :                              efinal,                 &
     397             :                              hin,         hilyr,     &
     398             :                              hsn,         hslyr,     &
     399           0 :                              zqin,        zqsn,      &
     400             :                              fbot,        Tbot,      &
     401             :                              flatn,       fsurfn,    &
     402             :                              fcondtopn,   fcondbotn,  &
     403             :                              fsnow,       hsn_new,   &
     404             :                              fhocnn,      evapn,     &
     405             :                              evapsn,      evapin,    &
     406             :                              meltt,       melts,     &
     407             :                              meltb,       &
     408             :                              congel,      snoice,    &
     409             :                              mlt_onset,   frz_onset, &
     410           0 :                              zSin,        sss,       &
     411   774714937 :                              dsnow)
     412   774714937 :       if (icepack_warnings_aborted(subname)) return
     413             : 
     414             :       !-----------------------------------------------------------------
     415             :       ! Check for energy conservation by comparing the change in energy
     416             :       ! to the net energy input
     417             :       !-----------------------------------------------------------------
     418             : 
     419             :       call conservation_check_vthermo(dt,                  &
     420             :                                       fsurfn,    flatn,    &
     421             :                                       fhocnn,    fswint,   &
     422             :                                       fsnow,     einit,    &
     423             :                                       einter,    efinal,   &
     424             :                                       fcondtopn, fcondbotn, &
     425   774714937 :                                       fadvocn,   fbot      )
     426   774714937 :       if (icepack_warnings_aborted(subname)) return
     427             : 
     428             :       !-----------------------------------------------------------------
     429             :       ! If prescribed ice, set hi back to old values
     430             :       !-----------------------------------------------------------------
     431             : 
     432             : #ifdef CESMCOUPLED
     433             :       if (present(prescribed_ice)) then
     434             :           if (prescribed_ice) then
     435             :             hin    = worki
     436             :             fhocnn = c0             ! for diagnostics
     437             :           endif
     438             :       endif
     439             : #endif
     440             : 
     441             :       !-----------------------------------------------------------------
     442             :       ! Compute fluxes of water and salt from ice to ocean.
     443             :       ! evapn < 0 => sublimation, evapn > 0 => condensation
     444             :       ! aerosol flux is accounted for in icepack_aerosol.F90
     445             :       !-----------------------------------------------------------------
     446             :             
     447   774714937 :       dhi = hin - worki
     448   774714937 :       dhs = hsn - works - hsn_new
     449             :       
     450   774714937 :       freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt
     451   774714937 :       fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
     452   774714937 :       fhocnn = fhocnn + fadvocn ! for ktherm=2 
     453             : 
     454   774714937 :       if (hin == c0) then
     455      105574 :          if (tr_pond_topo) fpond = fpond - aicen * apond * hpond
     456             :       endif
     457             : 
     458             :       !-----------------------------------------------------------------
     459             :       !  Given the vertical thermo state variables, compute the new ice 
     460             :       !   state variables.
     461             :       !-----------------------------------------------------------------
     462             : 
     463             :       call update_state_vthermo(nilyr,   nslyr,   &
     464             :                                 Tbot,    Tsf,     &     
     465             :                                 hin,     hsn,     &
     466           0 :                                 zqin,    zSin,    &
     467           0 :                                 zqsn,             &
     468             :                                 aicen,            &
     469   774714937 :                                 vicen,   vsnon)
     470   774714937 :       if (icepack_warnings_aborted(subname)) return
     471             : 
     472             :       !-----------------------------------------------------------------
     473             :       ! Reload passive tracer array
     474             :       !-----------------------------------------------------------------
     475             : 
     476             :     end subroutine thermo_vertical
     477             : 
     478             : !=======================================================================
     479             : !
     480             : ! Compute heat flux to bottom surface.
     481             : ! Compute fraction of ice that melts laterally.
     482             : !
     483             : ! authors C. M. Bitz, UW
     484             : !         William H. Lipscomb, LANL
     485             : !         Elizabeth C. Hunke, LANL
     486             : 
     487   722474472 :       subroutine frzmlt_bottom_lateral (dt,       ncat,     &
     488             :                                         nilyr,    nslyr,    &
     489             :                                         aice,     frzmlt,   &
     490  1444948944 :                                         vicen,    vsnon,    &
     491  1444948944 :                                         qicen,    qsnon,    &
     492             :                                         sst,      Tf,       & 
     493             :                                         ustar_min,          &
     494    45106800 :                                         fbot_xfer_type,     &
     495             :                                         strocnxT, strocnyT, &
     496             :                                         Tbot,     fbot,     &
     497             :                                         rside,    Cdn_ocn,  &
     498             :                                         fside)
     499             : 
     500             :       integer (kind=int_kind), intent(in) :: &
     501             :          ncat  , & ! number of thickness categories
     502             :          nilyr , & ! number of ice layers
     503             :          nslyr     ! number of snow layers
     504             : 
     505             :       real (kind=dbl_kind), intent(in) :: &
     506             :          dt                  ! time step
     507             : 
     508             :       real (kind=dbl_kind), intent(in) :: &
     509             :          aice    , & ! ice concentration
     510             :          frzmlt  , & ! freezing/melting potential (W/m^2)
     511             :          sst     , & ! sea surface temperature (C)
     512             :          Tf      , & ! freezing temperature (C)
     513             :          ustar_min,& ! minimum friction velocity for ice-ocean heat flux
     514             :          Cdn_ocn , & ! ocean-ice neutral drag coefficient
     515             :          strocnxT, & ! ice-ocean stress, x-direction
     516             :          strocnyT    ! ice-ocean stress, y-direction
     517             : 
     518             :       character (char_len), intent(in) :: &
     519             :          fbot_xfer_type  ! transfer coefficient type for ice-ocean heat flux
     520             : 
     521             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
     522             :          vicen   , & ! ice volume (m)
     523             :          vsnon       ! snow volume (m)
     524             : 
     525             :       real (kind=dbl_kind), dimension(:,:), intent(in) :: &
     526             :          qicen   , & ! ice layer enthalpy (J m-3)
     527             :          qsnon       ! snow layer enthalpy (J m-3)
     528             : 
     529             :       real (kind=dbl_kind), intent(out) :: &
     530             :          Tbot    , & ! ice bottom surface temperature (deg C)
     531             :          fbot    , & ! heat flux to ice bottom  (W/m^2)
     532             :          rside   , & ! fraction of ice that melts laterally
     533             :          fside       ! lateral heat flux (W/m^2)
     534             : 
     535             :       ! local variables
     536             : 
     537             :       integer (kind=int_kind) :: &
     538             :          n              , & ! thickness category index
     539             :          k                  ! layer index
     540             : 
     541             :       real (kind=dbl_kind) :: &
     542    45106800 :          etot    , & ! total energy in column
     543    45106800 :          qavg        ! average enthalpy in column (approximate)
     544             : 
     545             :       real (kind=dbl_kind) :: &
     546    45106800 :          deltaT    , & ! SST - Tbot >= 0
     547    45106800 :          ustar     , & ! skin friction velocity for fbot (m/s)
     548    45106800 :          wlat      , & ! lateral melt rate (m/s)
     549    45106800 :          xtmp          ! temporary variable
     550             : 
     551             :       ! Parameters for bottom melting
     552             : 
     553             :       real (kind=dbl_kind) :: &
     554    45106800 :          cpchr         ! -cp_ocn*rhow*exchange coefficient
     555             : 
     556             :       ! Parameters for lateral melting
     557             : 
     558             :       real (kind=dbl_kind), parameter :: &
     559             :          floediam = 300.0_dbl_kind, & ! effective floe diameter (m)
     560             :          floeshape = 0.66_dbl_kind , & ! constant from Steele (unitless)
     561             :          m1 = 1.6e-6_dbl_kind     , & ! constant from Maykut & Perovich
     562             :                                       ! (m/s/deg^(-m2))
     563             :          m2 = 1.36_dbl_kind           ! constant from Maykut & Perovich
     564             :                                       ! (unitless)
     565             : 
     566             :       character(len=*),parameter :: subname='(frzmlt_bottom_lateral)'
     567             : 
     568             :       !-----------------------------------------------------------------
     569             :       ! Identify grid cells where ice can melt.
     570             :       !-----------------------------------------------------------------
     571             : 
     572   722474472 :       rside = c0
     573   722474472 :       fside = c0
     574   722474472 :       Tbot  = Tf
     575   722474472 :       fbot  = c0
     576   722474472 :       wlat  = c0
     577             : 
     578   722474472 :       if (aice > puny .and. frzmlt < c0) then ! ice can melt
     579             :          
     580             :       !-----------------------------------------------------------------
     581             :       ! Use boundary layer theory for fbot.
     582             :       ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703.
     583             :       !-----------------------------------------------------------------
     584             : 
     585    66606681 :          deltaT = max((sst-Tbot),c0)
     586             :          
     587             :          ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
     588    66606681 :          ustar = sqrt (sqrt(strocnxT**2+strocnyT**2)/rhow)
     589    66606681 :          ustar = max (ustar,ustar_min)
     590             : 
     591    66606681 :          if (trim(fbot_xfer_type) == 'Cdn_ocn') then
     592             :             ! Note: Cdn_ocn has already been used for calculating ustar 
     593             :             ! (formdrag only) --- David Schroeder (CPOM)
     594     2897616 :             cpchr = -cp_ocn*rhow*Cdn_ocn
     595             :          else ! fbot_xfer_type == 'constant'
     596             :             ! 0.006 = unitless param for basal heat flx ala McPhee and Maykut
     597    63709065 :             cpchr = -cp_ocn*rhow*0.006_dbl_kind
     598             :          endif
     599             : 
     600    66606681 :          fbot = cpchr * deltaT * ustar ! < 0
     601    66606681 :          fbot = max (fbot, frzmlt) ! frzmlt < fbot < 0
     602             :             
     603             : !!! uncomment to use all frzmlt for standalone runs
     604             :    !     fbot = min (c0, frzmlt)
     605             : 
     606             :       !-----------------------------------------------------------------
     607             :       ! Compute rside.  See these references:
     608             :       !    Maykut and Perovich (1987): JGR, 92, 7032-7044
     609             :       !    Steele (1992): JGR, 97, 17,729-17,738
     610             :       !-----------------------------------------------------------------
     611             : 
     612    66606681 :          wlat = m1 * deltaT**m2 ! Maykut & Perovich
     613    66606681 :          rside = wlat*dt*pi/(floeshape*floediam) ! Steele
     614    66606681 :          rside = max(c0,min(rside,c1))
     615             : 
     616             :       !-----------------------------------------------------------------
     617             :       ! Compute heat flux associated with this value of rside.
     618             :       !-----------------------------------------------------------------
     619             : 
     620   397671208 :          do n = 1, ncat
     621             :             
     622   331064527 :             etot = c0
     623   331064527 :             qavg = c0
     624             :             
     625             :             ! melting energy/unit area in each column, etot < 0
     626             :             
     627   662129054 :             do k = 1, nslyr
     628   331064527 :                etot = etot + qsnon(k,n) * vsnon(n)/real(nslyr,kind=dbl_kind)
     629   662129054 :                qavg = qavg + qsnon(k,n)
     630             :             enddo
     631             :             
     632  2648516216 :             do k = 1, nilyr
     633  2317451689 :                etot = etot + qicen(k,n) * vicen(n)/real(nilyr,kind=dbl_kind)
     634  2648516216 :                qavg = qavg + qicen(k,n)
     635             :             enddo                  ! nilyr
     636             :             
     637             :             ! lateral heat flux, fside < 0
     638   397671208 :             if (tr_fsd) then ! floe size distribution
     639    13734570 :                fside = fside + wlat*qavg
     640             :             else             ! default floe size
     641   317329957 :                fside = fside + rside*etot/dt
     642             :             endif
     643             :             
     644             :          enddo                     ! n
     645             :          
     646             :       !-----------------------------------------------------------------
     647             :       ! Limit bottom and lateral heat fluxes if necessary.
     648             :       !-----------------------------------------------------------------
     649             :             
     650    66606681 :          xtmp = frzmlt/(fbot + fside + puny) 
     651    66606681 :          xtmp = min(xtmp, c1)
     652    66606681 :          fbot  = fbot  * xtmp
     653    66606681 :          rside = rside * xtmp
     654    66606681 :          fside = fside * xtmp
     655             :          
     656             :       endif
     657             : 
     658   722474472 :       end subroutine frzmlt_bottom_lateral
     659             : 
     660             : !=======================================================================
     661             : !
     662             : ! Given the state variables (vicen, vsnon, zqin, etc.),
     663             : ! compute variables needed for the vertical thermodynamics
     664             : ! (hin, hsn, zTin, etc.)
     665             : !
     666             : ! authors William H. Lipscomb, LANL
     667             : !         C. M. Bitz, UW
     668             : 
     669   774714937 :       subroutine init_vertical_profile(nilyr,    nslyr,    &
     670             :                                        aicen,    vicen,    &
     671             :                                        vsnon,              &
     672             :                                        hin,      hilyr,    &
     673             :                                        hsn,      hslyr,    &
     674   774714937 :                                        zqin,     zTin,     &
     675   774714937 :                                        zqsn,     zTsn,     &
     676   774714937 :                                        zSin,               &
     677             :                                        einit )
     678             : 
     679             :       integer (kind=int_kind), intent(in) :: &
     680             :          nilyr , & ! number of ice layers
     681             :          nslyr     ! number of snow layers
     682             : 
     683             :       real (kind=dbl_kind), intent(in) :: &
     684             :          aicen , & ! concentration of ice
     685             :          vicen , & ! volume per unit area of ice          (m)
     686             :          vsnon     ! volume per unit area of snow         (m)
     687             :  
     688             :       real (kind=dbl_kind), intent(out):: &
     689             :          hilyr       , & ! ice layer thickness
     690             :          hslyr       , & ! snow layer thickness
     691             :          einit           ! initial energy of melting (J m-2)
     692             :  
     693             :       real (kind=dbl_kind), intent(out):: &
     694             :          hin         , & ! ice thickness (m)
     695             :          hsn             ! snow thickness (m)
     696             : 
     697             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
     698             :          zqin        , & ! ice layer enthalpy (J m-3)
     699             :          zTin            ! internal ice layer temperatures
     700             : 
     701             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
     702             :          zSin            ! internal ice layer salinities
     703             :         
     704             :       real (kind=dbl_kind), dimension (:), &
     705             :          intent(inout) :: &
     706             :          zqsn        , & ! snow enthalpy
     707             :          zTsn            ! snow temperature
     708             : 
     709             :       ! local variables
     710             :       real (kind=dbl_kind), dimension(nilyr) :: &
     711  1101756569 :          Tmlts           ! melting temperature
     712             : 
     713             :       integer (kind=int_kind) :: &
     714             :          k               ! ice layer index
     715             : 
     716             :       real (kind=dbl_kind) :: &
     717    54309860 :          rnslyr,       & ! real(nslyr)
     718    54309860 :          Tmax            ! maximum allowed snow/ice temperature (deg C)
     719             : 
     720             :       logical (kind=log_kind) :: &   ! for vector-friendly error checks
     721             :          tsno_high   , & ! flag for zTsn > Tmax
     722             :          tice_high   , & ! flag for zTin > Tmlt
     723             :          tsno_low    , & ! flag for zTsn < Tmin
     724             :          tice_low        ! flag for zTin < Tmin
     725             : 
     726             :       character(len=*),parameter :: subname='(init_vertical_profile)'
     727             : 
     728             :       !-----------------------------------------------------------------
     729             :       ! Initialize
     730             :       !-----------------------------------------------------------------
     731             : 
     732   774714937 :       rnslyr = real(nslyr,kind=dbl_kind)
     733             : 
     734   774714937 :       tsno_high = .false.
     735   774714937 :       tice_high = .false.
     736   774714937 :       tsno_low  = .false.
     737   774714937 :       tice_low  = .false.
     738             : 
     739   774714937 :       einit = c0
     740             :  
     741             :       !-----------------------------------------------------------------
     742             :       ! Surface temperature, ice and snow thickness
     743             :       ! Initialize internal energy
     744             :       !-----------------------------------------------------------------
     745             : 
     746   774714937 :       hin    = vicen / aicen
     747   774714937 :       hsn    = vsnon / aicen
     748   774714937 :       hilyr    = hin / real(nilyr,kind=dbl_kind)
     749   774714937 :       hslyr    = hsn / rnslyr
     750             : 
     751             :       !-----------------------------------------------------------------
     752             :       ! Snow enthalpy and maximum allowed snow temperature
     753             :       ! If heat_capacity = F, zqsn and zTsn are used only for checking
     754             :       ! conservation.
     755             :       !-----------------------------------------------------------------
     756             : 
     757  1549429874 :       do k = 1, nslyr
     758             : 
     759             :       !-----------------------------------------------------------------
     760             :       ! Tmax based on the idea that dT ~ dq / (rhos*cp_ice)
     761             :       !                             dq ~ q dv / v
     762             :       !                             dv ~ puny = eps11
     763             :       ! where 'd' denotes an error due to roundoff.
     764             :       !-----------------------------------------------------------------
     765             : 
     766   774714937 :          if (hslyr > hs_min/rnslyr .and. heat_capacity) then
     767             :             ! zqsn < 0              
     768    34110740 :             Tmax = -zqsn(k)*puny*rnslyr / &
     769   582101107 :                  (rhos*cp_ice*vsnon)
     770             :          else
     771   192613830 :             zqsn  (k) = -rhos * Lfresh
     772   192613830 :             Tmax = puny
     773             :          endif
     774             : 
     775             :       !-----------------------------------------------------------------
     776             :       ! Compute snow temperatures from enthalpies.
     777             :       ! Note: zqsn <= -rhos*Lfresh, so zTsn <= 0.
     778             :       !-----------------------------------------------------------------
     779   774714937 :          zTsn(k) = (Lfresh + zqsn(k)/rhos)/cp_ice
     780             : 
     781             :       !-----------------------------------------------------------------
     782             :       ! Check for zTsn > Tmax (allowing for roundoff error) and zTsn < Tmin.
     783             :       !-----------------------------------------------------------------
     784  1549429874 :          if (zTsn(k) > Tmax) then
     785           0 :             tsno_high = .true.
     786   774714937 :          elseif (zTsn(k) < Tmin) then
     787           0 :             tsno_low  = .true.
     788             :          endif
     789             : 
     790             :       enddo                     ! nslyr
     791             : 
     792             :       !-----------------------------------------------------------------
     793             :       ! If zTsn is out of bounds, print diagnostics and exit.
     794             :       !-----------------------------------------------------------------
     795             : 
     796   774714937 :       if (tsno_high .and. heat_capacity) then
     797           0 :          do k = 1, nslyr
     798             : 
     799           0 :             if (hslyr > hs_min/rnslyr) then
     800           0 :                Tmax = -zqsn(k)*puny*rnslyr / &
     801           0 :                     (rhos*cp_ice*vsnon)
     802             :             else
     803           0 :                Tmax = puny
     804             :             endif
     805             : 
     806           0 :             if (zTsn(k) > Tmax) then
     807           0 :                write(warnstr,*) ' '
     808           0 :                call icepack_warnings_add(warnstr)
     809           0 :                write(warnstr,*) subname, 'Starting thermo, zTsn > Tmax'
     810           0 :                call icepack_warnings_add(warnstr)
     811           0 :                write(warnstr,*) subname, 'zTsn=',zTsn(k)
     812           0 :                call icepack_warnings_add(warnstr)
     813           0 :                write(warnstr,*) subname, 'Tmax=',Tmax
     814           0 :                call icepack_warnings_add(warnstr)
     815           0 :                write(warnstr,*) subname, 'zqsn',zqsn(k),-Lfresh*rhos,zqsn(k)+Lfresh*rhos
     816           0 :                call icepack_warnings_add(warnstr)
     817           0 :                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     818           0 :                call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTsn > Tmax" ) 
     819           0 :                return
     820             :             endif
     821             : 
     822             :          enddo                  ! nslyr
     823             :       endif                     ! tsno_high
     824             : 
     825   774714937 :       if (tsno_low .and. heat_capacity) then
     826           0 :          do k = 1, nslyr
     827             : 
     828           0 :             if (zTsn(k) < Tmin) then ! allowing for roundoff error
     829           0 :                write(warnstr,*) ' '
     830           0 :                call icepack_warnings_add(warnstr)
     831           0 :                write(warnstr,*) subname, 'Starting thermo, zTsn < Tmin'
     832           0 :                call icepack_warnings_add(warnstr)
     833           0 :                write(warnstr,*) subname, 'zTsn=', zTsn(k)
     834           0 :                call icepack_warnings_add(warnstr)
     835           0 :                write(warnstr,*) subname, 'Tmin=', Tmin
     836           0 :                call icepack_warnings_add(warnstr)
     837           0 :                write(warnstr,*) subname, 'zqsn', zqsn(k)
     838           0 :                call icepack_warnings_add(warnstr)
     839           0 :                write(warnstr,*) subname, hin
     840           0 :                call icepack_warnings_add(warnstr)
     841           0 :                write(warnstr,*) subname, hsn
     842           0 :                call icepack_warnings_add(warnstr)
     843           0 :                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     844           0 :                call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, zTsn < Tmin" ) 
     845           0 :                return
     846             :             endif
     847             : 
     848             :          enddo                  ! nslyr
     849             :       endif                     ! tsno_low
     850             : 
     851  1549429874 :       do k = 1, nslyr
     852             : 
     853   774714937 :          if (zTsn(k) > c0) then   ! correct roundoff error
     854     4009830 :             zTsn(k) = c0
     855     4009830 :             zqsn(k) = -rhos*Lfresh
     856             :          endif
     857             : 
     858             :       !-----------------------------------------------------------------
     859             :       ! initial energy per unit area of ice/snow, relative to 0 C
     860             :       !-----------------------------------------------------------------
     861  1549429874 :          einit = einit + hslyr*zqsn(k)
     862             : 
     863             :       enddo                     ! nslyr
     864             : 
     865  5517993338 :       do k = 1, nilyr
     866             : 
     867             :       !---------------------------------------------------------------------
     868             :       !  Use initial salinity profile for thin ice
     869             :       !---------------------------------------------------------------------
     870             : 
     871  4743278401 :          if (ktherm == 1 .and. zSin(k) < min_salin-puny) then
     872           0 :             write(warnstr,*) ' '
     873           0 :             call icepack_warnings_add(warnstr)
     874           0 :             write(warnstr,*) subname, 'Starting zSin < min_salin, layer', k
     875           0 :             call icepack_warnings_add(warnstr)
     876           0 :             write(warnstr,*) subname, 'zSin =', zSin(k)
     877           0 :             call icepack_warnings_add(warnstr)
     878           0 :             write(warnstr,*) subname, 'min_salin =', min_salin
     879           0 :             call icepack_warnings_add(warnstr)
     880           0 :             call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     881           0 :             call icepack_warnings_add(subname//" init_vertical_profile: Starting zSin < min_salin, layer" ) 
     882           0 :             return
     883             :          endif
     884             :          
     885  4743278401 :          if (ktherm == 2) then
     886  4265224733 :             Tmlts(k) = liquidus_temperature_mush(zSin(k))
     887             :          else
     888   478053668 :             Tmlts(k) = -zSin(k) * depressT
     889             :          endif
     890             : 
     891             :       !-----------------------------------------------------------------
     892             :       ! Compute ice enthalpy
     893             :       ! If heat_capacity = F, zqin and zTin are used only for checking
     894             :       ! conservation.
     895             :       !-----------------------------------------------------------------
     896             : 
     897             :       !-----------------------------------------------------------------
     898             :       ! Compute ice temperatures from enthalpies using quadratic formula
     899             :       !-----------------------------------------------------------------
     900             :          
     901  4743278401 :          if (ktherm == 2) then
     902  4265224733 :             zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k))
     903             :          else
     904   478053668 :             zTin(k) = calculate_Tin_from_qin(zqin(k),Tmlts(k))
     905             :          endif
     906             : 
     907  4743278401 :          if (l_brine) then
     908  4629990708 :             Tmax = Tmlts(k)
     909             :          else                ! fresh ice
     910   113287693 :             Tmax = -zqin(k)*puny/(rhos*cp_ice*vicen)
     911             :          endif
     912             : 
     913             :       !-----------------------------------------------------------------
     914             :       ! Check for zTin > Tmax and zTin < Tmin
     915             :       !-----------------------------------------------------------------
     916  4743278401 :          if (zTin(k) > Tmax) then
     917           0 :             tice_high = .true.
     918  4743278401 :          elseif (zTin(k) < Tmin) then
     919           0 :             tice_low  = .true.
     920             :          endif
     921             : 
     922             :       !-----------------------------------------------------------------
     923             :       ! If zTin is out of bounds, print diagnostics and exit.
     924             :       !-----------------------------------------------------------------
     925             : 
     926  4743278401 :          if (tice_high .and. heat_capacity) then
     927             : 
     928           0 :             if (l_brine) then
     929           0 :                Tmax = Tmlts(k)
     930             :             else             ! fresh ice
     931           0 :                Tmax = -zqin(k)*puny/(rhos*cp_ice*vicen)
     932             :             endif
     933             : 
     934           0 :             if (zTin(k) > Tmax) then
     935           0 :                write(warnstr,*) ' '
     936           0 :                call icepack_warnings_add(warnstr)
     937           0 :                write(warnstr,*) subname, 'Starting thermo, T > Tmax, layer', k
     938           0 :                call icepack_warnings_add(warnstr)
     939           0 :                write(warnstr,*) subname, 'k:', k
     940           0 :                call icepack_warnings_add(warnstr)
     941           0 :                write(warnstr,*) subname, 'zTin =',zTin(k),', Tmax=',Tmax
     942           0 :                call icepack_warnings_add(warnstr)
     943           0 :                write(warnstr,*) subname, 'zSin =',zSin(k)
     944           0 :                call icepack_warnings_add(warnstr)
     945           0 :                write(warnstr,*) subname, 'hin =',hin
     946           0 :                call icepack_warnings_add(warnstr)
     947           0 :                write(warnstr,*) subname, 'zqin =',zqin(k)
     948           0 :                call icepack_warnings_add(warnstr)
     949           0 :                write(warnstr,*) subname, 'qmlt=',enthalpy_of_melting(zSin(k))
     950           0 :                call icepack_warnings_add(warnstr)
     951           0 :                write(warnstr,*) subname, 'Tmlt=',Tmlts(k)
     952           0 :                call icepack_warnings_add(warnstr)
     953             :                
     954           0 :                if (ktherm == 2) then
     955           0 :                   zqin(k) = enthalpy_of_melting(zSin(k)) - c1
     956           0 :                   zTin(k) = icepack_mushy_temperature_mush(zqin(k),zSin(k))
     957           0 :                   write(warnstr,*) subname, 'Corrected quantities'
     958           0 :                   call icepack_warnings_add(warnstr)
     959           0 :                   write(warnstr,*) subname, 'zqin=',zqin(k)
     960           0 :                   call icepack_warnings_add(warnstr)
     961           0 :                   write(warnstr,*) subname, 'zTin=',zTin(k)
     962           0 :                   call icepack_warnings_add(warnstr)
     963             :                else
     964           0 :                   call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     965           0 :                   call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, T > Tmax, layer" ) 
     966           0 :                   return
     967             :                endif
     968             :             endif
     969             :          endif                  ! tice_high
     970             : 
     971  4743278401 :          if (tice_low .and. heat_capacity) then
     972             : 
     973           0 :             if (zTin(k) < Tmin) then
     974           0 :                write(warnstr,*) ' '
     975           0 :                call icepack_warnings_add(warnstr)
     976           0 :                write(warnstr,*) subname, 'Starting thermo T < Tmin, layer', k
     977           0 :                call icepack_warnings_add(warnstr)
     978           0 :                write(warnstr,*) subname, 'zTin =', zTin(k)
     979           0 :                call icepack_warnings_add(warnstr)
     980           0 :                write(warnstr,*) subname, 'Tmin =', Tmin
     981           0 :                call icepack_warnings_add(warnstr)
     982           0 :                call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
     983           0 :                call icepack_warnings_add(subname//" init_vertical_profile: Starting thermo, T < Tmin, layer" ) 
     984           0 :                return
     985             :             endif
     986             :          endif                  ! tice_low
     987             : 
     988             :       !-----------------------------------------------------------------
     989             :       ! correct roundoff error
     990             :       !-----------------------------------------------------------------
     991             : 
     992  4743278401 :          if (ktherm /= 2) then
     993             : 
     994   478053668 :             if (zTin(k) > c0) then
     995    20630681 :                zTin(k) = c0
     996    20630681 :                zqin(k) = -rhoi*Lfresh
     997             :             endif
     998             : 
     999             :          endif
    1000             : 
    1001             : ! echmod: is this necessary?
    1002             : !         if (ktherm == 1) then
    1003             : !               if (zTin(k)>= -zSin(k)*depressT) then
    1004             : !                   zTin(k) = -zSin(k)*depressT - puny
    1005             : !                   zqin(k) = -rhoi*cp_ocn*zSin(k)*depressT
    1006             : !               endif
    1007             : !         endif
    1008             : 
    1009             :       !-----------------------------------------------------------------
    1010             :       ! initial energy per unit area of ice/snow, relative to 0 C
    1011             :       !-----------------------------------------------------------------
    1012             : 
    1013  5517993338 :          einit = einit + hilyr*zqin(k) 
    1014             : 
    1015             :       enddo                     ! nilyr
    1016             : 
    1017             :       end subroutine init_vertical_profile
    1018             : 
    1019             : !=======================================================================
    1020             : !
    1021             : ! Compute growth and/or melting at the top and bottom surfaces.
    1022             : ! Convert snow to ice if necessary.
    1023             : !
    1024             : ! authors William H. Lipscomb, LANL
    1025             : !         C. M. Bitz, UW
    1026             : 
    1027   774714937 :       subroutine thickness_changes (nilyr,     nslyr,    &
    1028             :                                     dt,        yday,     &
    1029             :                                     efinal,              & 
    1030             :                                     hin,       hilyr,    &
    1031             :                                     hsn,       hslyr,    &
    1032   774714937 :                                     zqin,      zqsn,     &
    1033             :                                     fbot,      Tbot,     &
    1034             :                                     flatn,     fsurfn,   &
    1035             :                                     fcondtopn, fcondbotn, &
    1036             :                                     fsnow,     hsn_new,  &
    1037             :                                     fhocnn,    evapn,    &
    1038             :                                     evapsn,    evapin,   &
    1039             :                                     meltt,     melts,    &
    1040             :                                     meltb,     &
    1041             :                                     congel,    snoice,   &  
    1042             :                                     mlt_onset, frz_onset,&
    1043   774714937 :                                     zSin,      sss,      &
    1044             :                                     dsnow)
    1045             : 
    1046             :       integer (kind=int_kind), intent(in) :: &
    1047             :          nilyr , & ! number of ice layers
    1048             :          nslyr     ! number of snow layers
    1049             : 
    1050             :       real (kind=dbl_kind), intent(in) :: &
    1051             :          dt          , & ! time step
    1052             :          yday            ! day of the year
    1053             : 
    1054             :       real (kind=dbl_kind), intent(in) :: &
    1055             :          fbot        , & ! ice-ocean heat flux at bottom surface (W/m^2)
    1056             :          Tbot        , & ! ice bottom surface temperature (deg C)
    1057             :          fsnow       , & ! snowfall rate (kg m-2 s-1)
    1058             :          flatn       , & ! surface downward latent heat (W m-2)
    1059             :          fsurfn      , & ! net flux to top surface, excluding fcondtopn
    1060             :          fcondtopn       ! downward cond flux at top surface (W m-2)
    1061             : 
    1062             :       real (kind=dbl_kind), intent(inout) :: &
    1063             :          fcondbotn       ! downward cond flux at bottom surface (W m-2)
    1064             : 
    1065             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1066             :          zqin        , & ! ice layer enthalpy (J m-3)
    1067             :          zqsn            ! snow layer enthalpy (J m-3)
    1068             : 
    1069             :       real (kind=dbl_kind), intent(inout) :: &
    1070             :          hilyr       , & ! ice layer thickness (m)
    1071             :          hslyr           ! snow layer thickness (m)
    1072             : 
    1073             :       real (kind=dbl_kind), intent(inout) :: &
    1074             :          meltt       , & ! top ice melt             (m/step-->cm/day)
    1075             :          melts       , & ! snow melt                (m/step-->cm/day)
    1076             :          meltb       , & ! basal ice melt           (m/step-->cm/day)
    1077             :          congel      , & ! basal ice growth         (m/step-->cm/day)
    1078             :          snoice      , & ! snow-ice formation       (m/step-->cm/day)
    1079             :          dsnow       , & ! snow  formation          (m/step-->cm/day)
    1080             : !        iage        , & ! ice age (s)
    1081             :          mlt_onset   , & ! day of year that sfc melting begins
    1082             :          frz_onset       ! day of year that freezing begins (congel or frazil)
    1083             : 
    1084             :       real (kind=dbl_kind), intent(inout) :: &
    1085             :          hin         , & ! total ice thickness (m)
    1086             :          hsn             ! total snow thickness (m)
    1087             : 
    1088             :       real (kind=dbl_kind), intent(out):: &
    1089             :          efinal          ! final energy of melting (J m-2)
    1090             : 
    1091             :       real (kind=dbl_kind), intent(out):: &
    1092             :          fhocnn      , & ! fbot, corrected for any surplus energy (W m-2)
    1093             :          evapn       , & ! ice/snow mass sublimated/condensed (kg m-2 s-1)
    1094             :          evapsn      , & ! ice/snow mass sublimated/condensed over snow (kg m-2 s-1)
    1095             :          evapin          ! ice/snow mass sublimated/condensed over ice (kg m-2 s-1)
    1096             : 
    1097             :       real (kind=dbl_kind), intent(out):: &
    1098             :          hsn_new         ! thickness of new snow (m)
    1099             : 
    1100             :       ! changes to zSin in this subroutine are not reloaded into the
    1101             :       ! trcrn array for ktherm /= 2, so we could remove ktherm=2 conditionals
    1102             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1103             :          zSin            ! ice layer salinity (ppt)
    1104             : 
    1105             :       real (kind=dbl_kind), intent(in) :: &
    1106             :          sss             ! ocean salinity (PSU) 
    1107             : 
    1108             :       ! local variables
    1109             : 
    1110             :       integer (kind=int_kind) :: &
    1111             :          k               ! vertical index
    1112             : 
    1113             :       real (kind=dbl_kind) :: &
    1114    54309860 :          esub        , & ! energy for sublimation, > 0    (J m-2)
    1115    54309860 :          econ        , & ! energy for condensation, < 0   (J m-2)
    1116    54309860 :          etop_mlt    , & ! energy for top melting, > 0    (J m-2)
    1117    54309860 :          ebot_mlt    , & ! energy for bottom melting, > 0 (J m-2)
    1118    54309860 :          ebot_gro    , & ! energy for bottom growth, < 0  (J m-2)
    1119    54309860 :          emlt_atm    , & ! total energy of brine, from atmosphere (J m-2)
    1120    54309860 :          emlt_ocn        ! total energy of brine, to ocean        (J m-2)
    1121             : 
    1122             :       real (kind=dbl_kind) :: &
    1123    54309860 :          qbotmax     , & ! max enthalpy of ice growing at bottom
    1124    54309860 :          dhi         , & ! change in ice thickness
    1125    54309860 :          dhs         , & ! change in snow thickness
    1126    54309860 :          Ti          , & ! ice temperature
    1127    54309860 :          Ts          , & ! snow temperature
    1128    54309860 :          qbot        , & ! enthalpy of ice growing at bottom surface (J m-3)
    1129    54309860 :          qsub        , & ! energy/unit volume to sublimate ice/snow (J m-3)
    1130    54309860 :          hqtot       , & ! sum of h*q for two layers
    1131    54309860 :          wk1         , & ! temporary variable
    1132    54309860 :          zqsnew      , & ! enthalpy of new snow (J m-3)
    1133    54309860 :          hstot       , & ! snow thickness including new snow (m)
    1134    54309860 :          Tmlts           ! melting temperature
    1135             : 
    1136             :       real (kind=dbl_kind), dimension (nilyr+1) :: &
    1137  1822161646 :          zi1         , & ! depth of ice layer boundaries (m)
    1138  1822161646 :          zi2             ! adjusted depths, with equal hilyr (m)
    1139             : 
    1140             :       real (kind=dbl_kind), dimension (nslyr+1) :: &
    1141  1603739734 :          zs1         , & ! depth of snow layer boundaries (m)
    1142  1603739734 :          zs2             ! adjusted depths, with equal hslyr (m)
    1143             : 
    1144             :       real (kind=dbl_kind), dimension (nilyr) :: &
    1145  1101756569 :          dzi             ! ice layer thickness after growth/melting
    1146             : 
    1147             :       real (kind=dbl_kind), dimension (nslyr) :: &
    1148  1549429874 :          dzs             ! snow layer thickness after growth/melting
    1149             : 
    1150             :       real (kind=dbl_kind), dimension (nilyr) :: &
    1151  1767851786 :          qm          , & ! energy of melting (J m-3) = zqin in BL99 formulation
    1152  1156066429 :          qmlt            ! enthalpy of melted ice (J m-3) = zero in BL99 formulation
    1153             : 
    1154             :       real (kind=dbl_kind) :: &
    1155    54309860 :          qbotm       , &
    1156    54309860 :          qbotp       , &
    1157    54309860 :          qbot0
    1158             : 
    1159             :       character(len=*),parameter :: subname='(thickness_changes)'
    1160             : 
    1161             :       !-----------------------------------------------------------------
    1162             :       ! Initialize
    1163             :       !-----------------------------------------------------------------
    1164             : 
    1165   774714937 :       dhi = c0
    1166   774714937 :       dhs = c0
    1167   774714937 :       hsn_new  = c0
    1168             : 
    1169  5517993338 :       do k = 1, nilyr
    1170  5517993338 :          dzi(k) = hilyr
    1171             :       enddo
    1172             : 
    1173  1549429874 :       do k = 1, nslyr
    1174  1549429874 :          dzs(k) = hslyr
    1175             :       enddo
    1176             : 
    1177  5517993338 :       do k = 1, nilyr
    1178  4743278401 :          if (ktherm == 2) then
    1179  4265224733 :             qmlt(k) = enthalpy_of_melting(zSin(k))
    1180             :          else
    1181   478053668 :             qmlt(k) = c0
    1182             :          endif
    1183  4743278401 :          qm(k) = zqin(k) - qmlt(k)
    1184  4743278401 :          emlt_atm = c0
    1185  5517993338 :          emlt_ocn = c0
    1186             :       enddo
    1187             : 
    1188             :       !-----------------------------------------------------------------
    1189             :       ! For l_brine = false (fresh ice), check for temperatures > 0.
    1190             :       !  Melt ice or snow as needed to bring temperatures back to 0.
    1191             :       ! For l_brine = true, this should not be necessary.
    1192             :       !-----------------------------------------------------------------
    1193             : 
    1194   774714937 :       if (.not. l_brine) then 
    1195             : 
    1196   226575386 :          do k = 1, nslyr
    1197   113287693 :             Ts = (Lfresh + zqsn(k)/rhos) / cp_ice
    1198   226575386 :             if (Ts > c0) then
    1199           0 :                dhs = cp_ice*Ts*dzs(k) / Lfresh
    1200           0 :                dzs(k) = dzs(k) - dhs
    1201           0 :                zqsn(k) = -rhos*Lfresh
    1202             :             endif
    1203             :          enddo
    1204             : 
    1205   226575386 :          do k = 1, nilyr
    1206   113287693 :             Ti = (Lfresh + zqin(k)/rhoi) / cp_ice
    1207   226575386 :             if (Ti > c0) then
    1208           0 :                dhi = cp_ice*Ti*dzi(k) / Lfresh
    1209           0 :                dzi(k) = dzi(k) - dhi
    1210           0 :                zqin(k) = -rhoi*Lfresh
    1211             :             endif
    1212             :          enddo                  ! k
    1213             : 
    1214             :       endif                     ! .not. l_brine
    1215             : 
    1216             :       !-----------------------------------------------------------------
    1217             :       ! Compute energy available for sublimation/condensation, top melt,
    1218             :       ! and bottom growth/melt.
    1219             :       !-----------------------------------------------------------------
    1220             : 
    1221   774714937 :       wk1 = -flatn * dt
    1222   774714937 :       esub = max(wk1, c0)     ! energy for sublimation, > 0
    1223   774714937 :       econ = min(wk1, c0)     ! energy for condensation, < 0
    1224             : 
    1225   774714937 :       wk1 = (fsurfn - fcondtopn) * dt
    1226   774714937 :       etop_mlt = max(wk1, c0)           ! etop_mlt > 0
    1227             : 
    1228   774714937 :       wk1 = (fcondbotn - fbot) * dt
    1229   774714937 :       ebot_mlt = max(wk1, c0)           ! ebot_mlt > 0
    1230   774714937 :       ebot_gro = min(wk1, c0)           ! ebot_gro < 0
    1231             : 
    1232             :       !--------------------------------------------------------------
    1233             :       ! Condensation (evapn > 0)
    1234             :       ! Note: evapn here has unit of kg/m^2.  Divide by dt later.
    1235             :       ! This is the only case in which energy from the atmosphere
    1236             :       ! is used for changes in the brine energy (emlt_atm).
    1237             :       !--------------------------------------------------------------
    1238             : 
    1239   774714937 :       evapn = c0          ! initialize
    1240   774714937 :       evapsn = c0          ! initialize
    1241   774714937 :       evapin = c0          ! initialize
    1242             : 
    1243   774714937 :       if (hsn > puny) then    ! add snow with enthalpy zqsn(1)
    1244   643151396 :          dhs = econ / (zqsn(1) - rhos*Lvap) ! econ < 0, dhs > 0
    1245   643151396 :          dzs(1) = dzs(1) + dhs
    1246   643151396 :          evapn = evapn + dhs*rhos
    1247   643151396 :          evapsn = evapsn + dhs*rhos
    1248             :       else                        ! add ice with enthalpy zqin(1)
    1249   131563541 :          dhi = econ / (qm(1) - rhoi*Lvap) ! econ < 0, dhi > 0
    1250   131563541 :          dzi(1) = dzi(1) + dhi
    1251   131563541 :          evapn = evapn + dhi*rhoi
    1252   131563541 :          evapin = evapin + dhi*rhoi
    1253             :          ! enthalpy of melt water
    1254   131563541 :          emlt_atm = emlt_atm - qmlt(1) * dhi 
    1255             :       endif
    1256             : 
    1257             :       !--------------------------------------------------------------
    1258             :       ! Grow ice (bottom)
    1259             :       !--------------------------------------------------------------
    1260             : 
    1261   774714937 :       if (ktherm == 2) then
    1262             : 
    1263   609317819 :          qbotm = enthalpy_mush(Tbot, sss)
    1264   609317819 :          qbotp = -Lfresh * rhoi * (c1 - phi_i_mushy)
    1265   609317819 :          qbot0 = qbotm - qbotp
    1266             : 
    1267   609317819 :          dhi = ebot_gro / qbotp     ! dhi > 0
    1268             : 
    1269   609317819 :          hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbotm
    1270   609317819 :          hstot = dzi(nilyr)*zSin(nilyr) + dhi*sss
    1271   609317819 :          emlt_ocn = emlt_ocn - qbot0 * dhi
    1272             : 
    1273             :       else
    1274             : 
    1275   165397118 :          Tmlts = -zSin(nilyr) * depressT 
    1276             : 
    1277             :          ! enthalpy of new ice growing at bottom surface
    1278   165397118 :          if (heat_capacity) then
    1279    52109425 :             if (l_brine) then
    1280    52109425 :                qbotmax = -p5*rhoi*Lfresh  ! max enthalpy of ice growing at bottom
    1281             :                qbot = -rhoi * (cp_ice * (Tmlts-Tbot) &
    1282             :                     + Lfresh * (c1-Tmlts/Tbot) &
    1283    52109425 :                     - cp_ocn * Tmlts)
    1284    52109425 :                qbot = min (qbot, qbotmax)      ! in case Tbot is close to Tmlt
    1285             :             else
    1286           0 :                qbot = -rhoi * (-cp_ice * Tbot + Lfresh)
    1287             :             endif
    1288             :          else   ! zero layer
    1289   113287693 :             qbot = -rhoi * Lfresh
    1290             :          endif
    1291             : 
    1292   165397118 :          dhi = ebot_gro / qbot     ! dhi > 0
    1293             : 
    1294   165397118 :          hqtot = dzi(nilyr)*zqin(nilyr) + dhi*qbot
    1295   165397118 :          hstot = c0
    1296             :       endif ! ktherm
    1297             : 
    1298   774714937 :       dzi(nilyr) = dzi(nilyr) + dhi
    1299   774714937 :       if (dzi(nilyr) > puny) then
    1300   774714937 :          zqin(nilyr) = hqtot / dzi(nilyr)
    1301   774714937 :          if (ktherm == 2) then
    1302   609317819 :             zSin(nilyr) = hstot / dzi(nilyr)
    1303   609317819 :             qmlt(nilyr) = enthalpy_of_melting(zSin(nilyr))
    1304             :          else
    1305   165397118 :             qmlt(nilyr) = c0
    1306             :          endif
    1307   774714937 :          qm(nilyr) = zqin(nilyr) - qmlt(nilyr)
    1308             :       endif
    1309             : 
    1310             :       ! update ice age due to freezing (new ice age = dt)
    1311             :       !         if (tr_iage) &
    1312             :       !            iage = (iage*hin + dt*dhi) / (hin + dhi)
    1313             : 
    1314             :       ! history diagnostics
    1315   774714937 :       congel = congel + dhi
    1316   774714937 :       if (dhi > puny .and. frz_onset < puny) &
    1317      408439 :            frz_onset = yday
    1318             : 
    1319  1549429874 :       do k = 1, nslyr
    1320             : 
    1321             :          !--------------------------------------------------------------
    1322             :          ! Remove internal snow melt 
    1323             :          !--------------------------------------------------------------
    1324             :          
    1325   774714937 :          if (ktherm == 2 .and. zqsn(k) > -rhos * Lfresh) then
    1326             : 
    1327     6403933 :             dhs = max(-dzs(k), &
    1328    48471814 :                 -((zqsn(k) + rhos*Lfresh) / (rhos*Lfresh)) * dzs(k))
    1329    48471814 :             dzs(k) = dzs(k) + dhs
    1330    48471814 :             zqsn(k) = -rhos * Lfresh
    1331    48471814 :             melts = melts - dhs
    1332             :             ! delta E = zqsn(k) + rhos * Lfresh
    1333             : 
    1334             :          endif
    1335             : 
    1336             :          !--------------------------------------------------------------
    1337             :          ! Sublimation of snow (evapn < 0)
    1338             :          !--------------------------------------------------------------
    1339             : 
    1340   774714937 :          qsub = zqsn(k) - rhos*Lvap ! qsub < 0
    1341   774714937 :          dhs  = max (-dzs(k), esub/qsub)  ! esub > 0, dhs < 0
    1342   774714937 :          dzs(k) = dzs(k) + dhs
    1343   774714937 :          esub = esub - dhs*qsub
    1344   774714937 :          esub = max(esub, c0)   ! in case of roundoff error
    1345   774714937 :          evapn = evapn + dhs*rhos
    1346   774714937 :          evapsn = evapsn + dhs*rhos
    1347             : 
    1348             :          !--------------------------------------------------------------
    1349             :          ! Melt snow (top)
    1350             :          !--------------------------------------------------------------
    1351             : 
    1352   774714937 :          dhs = max(-dzs(k), etop_mlt/zqsn(k))
    1353   774714937 :          dzs(k) = dzs(k) + dhs         ! zqsn < 0, dhs < 0
    1354   774714937 :          etop_mlt = etop_mlt - dhs*zqsn(k)
    1355   774714937 :          etop_mlt = max(etop_mlt, c0) ! in case of roundoff error
    1356             : 
    1357             :          ! history diagnostics
    1358   774714937 :          if (dhs < -puny .and. mlt_onset < puny) &
    1359      155883 :               mlt_onset = yday
    1360  1549429874 :          melts = melts - dhs
    1361             : 
    1362             :       enddo                     ! nslyr
    1363             : 
    1364  5517993338 :       do k = 1, nilyr
    1365             : 
    1366             :          !--------------------------------------------------------------
    1367             :          ! Sublimation of ice (evapn < 0)
    1368             :          !--------------------------------------------------------------
    1369             : 
    1370  4743278401 :          qsub = qm(k) - rhoi*Lvap              ! qsub < 0
    1371  4743278401 :          dhi  = max (-dzi(k), esub/qsub) ! esub < 0, dhi < 0
    1372  4743278401 :          dzi(k) = dzi(k) + dhi
    1373  4743278401 :          esub = esub - dhi*qsub
    1374  4743278401 :          esub = max(esub, c0)
    1375  4743278401 :          evapn = evapn + dhi*rhoi
    1376  4743278401 :          evapin = evapin + dhi*rhoi
    1377  4743278401 :          emlt_ocn = emlt_ocn - qmlt(k) * dhi 
    1378             : 
    1379             :          !--------------------------------------------------------------
    1380             :          ! Melt ice (top)
    1381             :          !--------------------------------------------------------------
    1382             :    
    1383  4743278401 :          if (qm(k) < c0) then
    1384  4743239690 :             dhi = max(-dzi(k), etop_mlt/qm(k))
    1385             :          else
    1386       38711 :             qm(k) = c0
    1387       38711 :             dhi = -dzi(k)
    1388             :          endif
    1389  4743278401 :          emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
    1390             : 
    1391  4743278401 :          dzi(k) = dzi(k) + dhi         ! zqin < 0, dhi < 0
    1392  4743278401 :          etop_mlt = max(etop_mlt - dhi*qm(k), c0)
    1393             : 
    1394             :          ! history diagnostics
    1395  4743278401 :          if (dhi < -puny .and. mlt_onset < puny) &
    1396       42495 :               mlt_onset = yday
    1397  5517993338 :          meltt = meltt - dhi
    1398             : 
    1399             :       enddo                     ! nilyr
    1400             : 
    1401  5517993338 :       do k = nilyr, 1, -1
    1402             : 
    1403             :          !--------------------------------------------------------------
    1404             :          ! Melt ice (bottom)
    1405             :          !--------------------------------------------------------------
    1406             : 
    1407  4743278401 :          if (qm(k) < c0) then
    1408  4743239690 :             dhi = max(-dzi(k), ebot_mlt/qm(k))
    1409             :          else
    1410       38711 :             qm(k) = c0
    1411       38711 :             dhi = -dzi(k)
    1412             :          endif
    1413  4743278401 :          emlt_ocn = emlt_ocn - max(zqin(k),qmlt(k)) * dhi
    1414             : 
    1415  4743278401 :          dzi(k) = dzi(k) + dhi         ! zqin < 0, dhi < 0
    1416  4743278401 :          ebot_mlt = max(ebot_mlt - dhi*qm(k), c0)
    1417             : 
    1418             :          ! history diagnostics 
    1419  5517993338 :          meltb = meltb -dhi
    1420             : 
    1421             :       enddo                     ! nilyr
    1422             : 
    1423  1549429874 :       do k = nslyr, 1, -1
    1424             : 
    1425             :          !--------------------------------------------------------------
    1426             :          ! Melt snow (only if all the ice has melted)
    1427             :          !--------------------------------------------------------------
    1428             :          
    1429   774714937 :          dhs = max(-dzs(k), ebot_mlt/zqsn(k))
    1430   774714937 :          dzs(k) = dzs(k) + dhs         ! zqsn < 0, dhs < 0
    1431   774714937 :          ebot_mlt = ebot_mlt - dhs*zqsn(k)
    1432  1549429874 :          ebot_mlt = max(ebot_mlt, c0)
    1433             : 
    1434             :       enddo                     ! nslyr
    1435             : 
    1436             :       !-----------------------------------------------------------------
    1437             :       ! Compute heat flux used by the ice (<=0).
    1438             :       ! fhocn is the available ocean heat that is left after use by ice
    1439             :       !-----------------------------------------------------------------
    1440             : 
    1441             :       fhocnn = fbot &
    1442   774714937 :              + (esub + etop_mlt + ebot_mlt)/dt
    1443             : 
    1444             : !---!-----------------------------------------------------------------
    1445             : !---! Add new snowfall at top surface.
    1446             : !---!-----------------------------------------------------------------
    1447             : 
    1448             :       !----------------------------------------------------------------
    1449             :       ! NOTE: If heat flux diagnostics are to work, new snow should
    1450             :       !       have T = 0 (i.e. q = -rhos*Lfresh) and should not be
    1451             :       !       converted to rain.
    1452             :       !----------------------------------------------------------------
    1453             : 
    1454   774714937 :       if (fsnow > c0) then
    1455             : 
    1456   418869473 :          hsn_new = fsnow/rhos * dt
    1457   418869473 :          zqsnew = -rhos*Lfresh
    1458   418869473 :          hstot = dzs(1) + hsn_new
    1459             : 
    1460   418869473 :          if (hstot > c0) then
    1461    24474633 :             zqsn(1) =  (dzs(1) * zqsn(1) &
    1462   418869473 :                     + hsn_new * zqsnew) / hstot
    1463             :             ! avoid roundoff errors
    1464   418869473 :             zqsn(1) = min(zqsn(1), -rhos*Lfresh)
    1465             : 
    1466   418869473 :             dzs(1) = hstot
    1467             :          endif
    1468             :       endif
    1469             : 
    1470             :     !-----------------------------------------------------------------
    1471             :     ! Find the new ice and snow thicknesses.
    1472             :     !-----------------------------------------------------------------
    1473             : 
    1474   774714937 :       hin = c0
    1475   774714937 :       hsn = c0
    1476             : 
    1477  5517993338 :       do k = 1, nilyr
    1478  5517993338 :          hin = hin + dzi(k)
    1479             :       enddo                     ! k
    1480             : 
    1481  1549429874 :       do k = 1, nslyr
    1482   774714937 :          hsn = hsn + dzs(k)
    1483  1549429874 :          dsnow = dsnow + dzs(k) - hslyr  
    1484             :       enddo                     ! k
    1485             : 
    1486             :     !-------------------------------------------------------------------
    1487             :     ! Convert snow to ice if snow lies below freeboard.
    1488             :     !-------------------------------------------------------------------
    1489             : 
    1490   774714937 :       if (ktherm /= 2) &
    1491             :          call freeboard (nslyr, &
    1492             :                          snoice, &
    1493             :                          hin,      hsn,      &
    1494           0 :                          zqin,     zqsn,     &
    1495           0 :                          dzi,      dzs,      &
    1496   165397118 :                          dsnow)
    1497   774714937 :          if (icepack_warnings_aborted(subname)) return
    1498             : 
    1499             : !---!-------------------------------------------------------------------
    1500             : !---! Repartition the ice and snow into equal-thickness layers,
    1501             : !---! conserving energy.
    1502             : !---!-------------------------------------------------------------------
    1503             : 
    1504             :       !-----------------------------------------------------------------
    1505             :       ! Compute desired layer thicknesses.
    1506             :       !-----------------------------------------------------------------
    1507             :  
    1508   774714937 :       if (hin > c0) then
    1509   774609363 :          hilyr = hin / real(nilyr,kind=dbl_kind)
    1510             :       else
    1511      105574 :          hin = c0
    1512      105574 :          hilyr = c0
    1513             :       endif
    1514   774714937 :       if (hsn > c0) then
    1515   631947991 :          hslyr = hsn / real(nslyr,kind=dbl_kind)
    1516             :       else
    1517   142766946 :          hsn = c0
    1518   142766946 :          hslyr = c0
    1519             :       endif
    1520             : 
    1521             :       !-----------------------------------------------------------------
    1522             :       ! Compute depths zi1 of old layers (unequal thickness).
    1523             :       ! Compute depths zi2 of new layers (equal thickness).
    1524             :       !-----------------------------------------------------------------
    1525             : 
    1526   774714937 :       zi1(1) = c0
    1527   774714937 :       zi1(1+nilyr) = hin
    1528             : 
    1529   774714937 :       zi2(1) = c0
    1530   774714937 :       zi2(1+nilyr) = hin
    1531             : 
    1532   774714937 :       if (heat_capacity) then
    1533             : 
    1534  4629990708 :          do k = 1, nilyr-1
    1535  3968563464 :             zi1(k+1) = zi1(k) + dzi(k)
    1536  4629990708 :             zi2(k+1) = zi2(k) + hilyr
    1537             :          enddo
    1538             : 
    1539             :         !-----------------------------------------------------------------
    1540             :         ! Conserving energy, compute the enthalpy of the new equal layers.
    1541             :         !-----------------------------------------------------------------
    1542             :             
    1543             :          call adjust_enthalpy (nilyr,              &
    1544           0 :                                zi1,      zi2,      &
    1545             :                                hilyr,    hin,      &
    1546   661427244 :                                zqin)   
    1547   661427244 :          if (icepack_warnings_aborted(subname)) return
    1548             : 
    1549   661427244 :          if (ktherm == 2) &
    1550             :               call adjust_enthalpy (nilyr,              &
    1551           0 :                                     zi1,      zi2,      &
    1552             :                                     hilyr,    hin,      &
    1553   609317819 :                                     zSin)   
    1554   661427244 :          if (icepack_warnings_aborted(subname)) return
    1555             : 
    1556             :       else ! zero layer (nilyr=1)
    1557             : 
    1558   113287693 :          zqin(1) = -rhoi * Lfresh
    1559   113287693 :          zqsn(1) = -rhos * Lfresh
    1560             :        
    1561             :       endif
    1562             : 
    1563   774714937 :       if (nslyr > 1) then
    1564             : 
    1565             :       !-----------------------------------------------------------------
    1566             :       ! Compute depths zs1 of old layers (unequal thickness).
    1567             :       ! Compute depths zs2 of new layers (equal thickness).
    1568             :       !-----------------------------------------------------------------
    1569             : 
    1570           0 :          zs1(1) = c0
    1571           0 :          zs1(1+nslyr) = hsn
    1572             :          
    1573           0 :          zs2(1) = c0
    1574           0 :          zs2(1+nslyr) = hsn
    1575             :          
    1576           0 :          do k = 1, nslyr-1
    1577           0 :             zs1(k+1) = zs1(k) + dzs(k)
    1578           0 :             zs2(k+1) = zs2(k) + hslyr
    1579             :          enddo
    1580             : 
    1581             :       !-----------------------------------------------------------------
    1582             :       ! Conserving energy, compute the enthalpy of the new equal layers.
    1583             :       !-----------------------------------------------------------------
    1584             : 
    1585             :          call adjust_enthalpy (nslyr,              &
    1586           0 :                                zs1,      zs2,      &
    1587             :                                hslyr,    hsn,      &
    1588           0 :                                zqsn)   
    1589           0 :          if (icepack_warnings_aborted(subname)) return
    1590             : 
    1591             :       endif   ! nslyr > 1
    1592             : 
    1593             :       !-----------------------------------------------------------------
    1594             :       ! Remove very thin snow layers (ktherm = 2)
    1595             :       !-----------------------------------------------------------------
    1596             : 
    1597   774714937 :       if (ktherm == 2) then
    1598  1218635638 :          do k = 1, nslyr
    1599  1218635638 :             if (hsn <= puny) then
    1600             :                fhocnn = fhocnn &
    1601    70334279 :                       + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
    1602    70334279 :                zqsn(k) = -rhos*Lfresh
    1603    70334279 :                hslyr = c0
    1604             :             endif
    1605             :          enddo
    1606             :       endif
    1607             : 
    1608             :       !-----------------------------------------------------------------
    1609             :       ! Compute final ice-snow energy, including the energy of
    1610             :       !  sublimated/condensed ice.
    1611             :       !-----------------------------------------------------------------
    1612             : 
    1613   774714937 :       efinal = -evapn*Lvap
    1614   774714937 :       evapn =  evapn/dt
    1615   774714937 :       evapsn =  evapsn/dt
    1616   774714937 :       evapin =  evapin/dt
    1617             : 
    1618  1549429874 :       do k = 1, nslyr
    1619  1549429874 :          efinal = efinal + hslyr*zqsn(k)
    1620             :       enddo
    1621             : 
    1622  5517993338 :       do k = 1, nilyr
    1623  5517993338 :          efinal = efinal + hilyr*zqin(k)
    1624             :       enddo                     ! k
    1625             : 
    1626   774714937 :       if (ktherm < 2) then
    1627   165397118 :          emlt_atm = c0
    1628   165397118 :          emlt_ocn = c0
    1629             :       endif
    1630             : 
    1631             :       ! melt water is no longer zero enthalpy with ktherm=2
    1632   774714937 :       fhocnn = fhocnn + emlt_ocn/dt
    1633   774714937 :       efinal = efinal + emlt_atm ! for conservation check
    1634             : 
    1635             :       end subroutine thickness_changes
    1636             : 
    1637             : !=======================================================================
    1638             : !
    1639             : ! If there is enough snow to lower the ice/snow interface below
    1640             : ! sea level, convert enough snow to ice to bring the interface back
    1641             : ! to sea level.
    1642             : !
    1643             : ! authors William H. Lipscomb, LANL
    1644             : !         Elizabeth C. Hunke, LANL
    1645             : 
    1646   165397118 :       subroutine freeboard (nslyr, &
    1647             :                             snoice,             &
    1648             :                             hin,      hsn,      &
    1649   330794236 :                             zqin,     zqsn,     &
    1650   330794236 :                             dzi,      dzs,      &
    1651             :                             dsnow)
    1652             : 
    1653             :       integer (kind=int_kind), intent(in) :: &
    1654             :          nslyr     ! number of snow layers
    1655             : 
    1656             : !     real (kind=dbl_kind), intent(in) :: &
    1657             : !        dt      ! time step
    1658             : 
    1659             :       real (kind=dbl_kind), &
    1660             :          intent(inout) :: &
    1661             :          snoice  , & ! snow-ice formation       (m/step-->cm/day)
    1662             :          dsnow       ! change in snow thickness after snow-ice formation (m)
    1663             : !        iage        ! ice age (s)
    1664             : 
    1665             :       real (kind=dbl_kind), &
    1666             :          intent(inout) :: &
    1667             :          hin     , & ! ice thickness (m)
    1668             :          hsn         ! snow thickness (m)
    1669             : 
    1670             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1671             :          zqsn        ! snow layer enthalpy (J m-3)
    1672             : 
    1673             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1674             :          zqin    , & ! ice layer enthalpy (J m-3)
    1675             :          dzi     , & ! ice layer thicknesses (m)
    1676             :          dzs         ! snow layer thicknesses (m)
    1677             : 
    1678             :       ! local variables
    1679             : 
    1680             :       integer (kind=int_kind) :: &
    1681             :          k               ! vertical index
    1682             : 
    1683             :       real (kind=dbl_kind) :: &
    1684    21858357 :          dhin        , & ! change in ice thickness (m)
    1685    21858357 :          dhsn        , & ! change in snow thickness (m)
    1686    21858357 :          hqs             ! sum of h*q for snow (J m-2)
    1687             : 
    1688             :       real (kind=dbl_kind) :: &
    1689    21858357 :          wk1         , & ! temporary variable
    1690    21858357 :          dhs             ! snow to remove from layer (m)
    1691             : 
    1692             :       character(len=*),parameter :: subname='(freeboard)'
    1693             : 
    1694             :       !-----------------------------------------------------------------
    1695             :       ! Determine whether snow lies below freeboard.
    1696             :       !-----------------------------------------------------------------
    1697             :       
    1698   165397118 :       dhin = c0
    1699   165397118 :       dhsn = c0
    1700   165397118 :       hqs  = c0
    1701             : 
    1702   165397118 :       wk1 = hsn - hin*(rhow-rhoi)/rhos
    1703             : 
    1704   165397118 :       if (wk1 > puny .and. hsn > puny) then  ! snow below freeboard
    1705     1431289 :          dhsn = min(wk1*rhoi/rhow, hsn) ! snow to remove
    1706     1431289 :          dhin = dhsn * rhos/rhoi        ! ice to add
    1707             :       endif
    1708             : 
    1709             :       !-----------------------------------------------------------------
    1710             :       ! Adjust snow layer thickness.
    1711             :       ! Compute energy to transfer from snow to ice.
    1712             :       !-----------------------------------------------------------------
    1713             : 
    1714   330794236 :       do k = nslyr, 1, -1
    1715   330794236 :          if (dhin > puny) then
    1716     1431289 :             dhs = min(dhsn, dzs(k)) ! snow to remove from layer
    1717     1431289 :             hsn = hsn - dhs
    1718     1431289 :             dsnow = dsnow -dhs   !new snow addition term 
    1719     1431289 :             dzs(k) = dzs(k) - dhs
    1720     1431289 :             dhsn = dhsn - dhs
    1721     1431289 :             dhsn = max(dhsn,c0)
    1722     1431289 :             hqs = hqs + dhs * zqsn(k)
    1723             :          endif               ! dhin > puny
    1724             :       enddo
    1725             : 
    1726             :       !-----------------------------------------------------------------
    1727             :       ! Transfer volume and energy from snow to top ice layer.
    1728             :       !-----------------------------------------------------------------
    1729             : 
    1730   165397118 :       if (dhin > puny) then
    1731             :          ! update ice age due to freezing (new ice age = dt)
    1732             :          !            if (tr_iage) &
    1733             :          !               iage = (iage*hin+dt*dhin)/(hin+dhin)
    1734             :          
    1735     1431289 :          wk1 = dzi(1) + dhin
    1736     1431289 :          hin = hin + dhin
    1737     1431289 :          zqin(1) = (dzi(1)*zqin(1) + hqs) / wk1
    1738     1431289 :          dzi(1) = wk1
    1739             : 
    1740             :          ! history diagnostic
    1741     1431289 :          snoice = snoice + dhin
    1742             :       endif               ! dhin > puny
    1743             : 
    1744   165397118 :       end subroutine freeboard
    1745             : 
    1746             : !=======================================================================
    1747             : !
    1748             : ! Conserving energy, compute the new enthalpy of equal-thickness ice
    1749             : ! or snow layers.
    1750             : !
    1751             : ! authors William H. Lipscomb, LANL
    1752             : !         C. M. Bitz, UW
    1753             : 
    1754  1270745063 :       subroutine adjust_enthalpy (nlyr,               &
    1755  1270745063 :                                   z1,       z2,       &
    1756             :                                   hlyr,     hn,       &
    1757  1270745063 :                                   qn)
    1758             : 
    1759             :       integer (kind=int_kind), intent(in) :: &
    1760             :          nlyr            ! number of layers (nilyr or nslyr)
    1761             : 
    1762             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
    1763             :          z1          , & ! interface depth for old, unequal layers (m)
    1764             :          z2              ! interface depth for new, equal layers (m)
    1765             : 
    1766             :       real (kind=dbl_kind), intent(in) :: &
    1767             :          hlyr            ! new layer thickness (m)
    1768             : 
    1769             :       real (kind=dbl_kind), intent(in) :: &
    1770             :          hn              ! total thickness (m)
    1771             : 
    1772             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1773             :          qn              ! layer quantity (enthalpy, salinity...)
    1774             : 
    1775             :       ! local variables
    1776             : 
    1777             :       integer (kind=int_kind) :: &
    1778             :          k, k1, k2       ! vertical indices
    1779             : 
    1780             :       real (kind=dbl_kind) :: &
    1781    68855155 :          hovlp           ! overlap between old and new layers (m)
    1782             : 
    1783             :       real (kind=dbl_kind) :: &
    1784    68855155 :          rhlyr           ! 1./hlyr
    1785             : 
    1786             :       real (kind=dbl_kind), dimension (nlyr) :: &
    1787  2954621056 :          hq              ! h * q for a layer
    1788             : 
    1789             :       character(len=*),parameter :: subname='(adjust_enthalpy)'
    1790             : 
    1791             :       !-----------------------------------------------------------------
    1792             :       ! Compute reciprocal layer thickness.
    1793             :       !-----------------------------------------------------------------
    1794             : 
    1795  1270745063 :       rhlyr = c0
    1796  1270745063 :       if (hn > puny) rhlyr = c1 / hlyr
    1797             : 
    1798             :       !-----------------------------------------------------------------
    1799             :       ! Compute h*q for new layers (k2) given overlap with old layers (k1)
    1800             :       !-----------------------------------------------------------------
    1801             : 
    1802 10165960504 :       do k2 = 1, nlyr
    1803 10165960504 :          hq(k2) = c0
    1804             :       enddo                     ! k
    1805  1270745063 :       k1 = 1
    1806  1270745063 :       k2 = 1
    1807 17789151197 :       do while (k1 <= nlyr .and. k2 <= nlyr)
    1808  1790006594 :          hovlp = min (z1(k1+1), z2(k2+1)) &
    1809 18308412728 :                - max (z1(k1),   z2(k2))
    1810 16518406134 :          hovlp = max (hovlp, c0)
    1811 16518406134 :          hq(k2) = hq(k2) + hovlp*qn(k1)
    1812 17789151197 :          if (z1(k1+1) > z2(k2+1)) then
    1813  7623211275 :             k2 = k2 + 1
    1814             :          else
    1815  8895194859 :             k1 = k1 + 1
    1816             :          endif
    1817             :       enddo                  ! while
    1818             : 
    1819             :       !-----------------------------------------------------------------
    1820             :       ! Compute new enthalpies.
    1821             :       !-----------------------------------------------------------------
    1822             : 
    1823 10165960504 :       do k = 1, nlyr
    1824 10165960504 :          qn(k) = hq(k) * rhlyr
    1825             :       enddo                     ! k
    1826             : 
    1827  1270745063 :       end subroutine adjust_enthalpy
    1828             : 
    1829             : !=======================================================================
    1830             : !
    1831             : ! Check for energy conservation by comparing the change in energy
    1832             : ! to the net energy input.
    1833             : !
    1834             : ! authors William H. Lipscomb, LANL
    1835             : !         C. M. Bitz, UW
    1836             : !         Adrian K. Turner, LANL
    1837             : 
    1838   774714937 :       subroutine conservation_check_vthermo(dt,                 &
    1839             :                                             fsurfn,   flatn,    &
    1840             :                                             fhocnn,   fswint,   &
    1841             :                                             fsnow,              &
    1842             :                                             einit,    einter,   &
    1843             :                                             efinal,             &
    1844             :                                             fcondtopn,fcondbotn, &
    1845             :                                             fadvocn,  fbot      )
    1846             : 
    1847             :       real (kind=dbl_kind), intent(in) :: &
    1848             :          dt              ! time step
    1849             : 
    1850             :       real (kind=dbl_kind), intent(in) :: &
    1851             :          fsurfn      , & ! net flux to top surface, excluding fcondtopn
    1852             :          flatn       , & ! surface downward latent heat (W m-2)
    1853             :          fhocnn      , & ! fbot, corrected for any surplus energy
    1854             :          fswint      , & ! SW absorbed in ice interior, below surface (W m-2)
    1855             :          fsnow       , & ! snowfall rate (kg m-2 s-1)
    1856             :          fcondtopn   , &
    1857             :          fadvocn     , &
    1858             :          fbot           
    1859             : 
    1860             :       real (kind=dbl_kind), intent(in) :: &
    1861             :          einit       , & ! initial energy of melting (J m-2)
    1862             :          einter      , & ! intermediate energy of melting (J m-2)
    1863             :          efinal      , & ! final energy of melting (J m-2)
    1864             :          fcondbotn
    1865             : 
    1866             :       ! local variables
    1867             : 
    1868             :       real (kind=dbl_kind) :: &
    1869    54309860 :          einp        , & ! energy input during timestep (J m-2)
    1870    54309860 :          ferr            ! energy conservation error (W m-2)
    1871             : 
    1872             :       character(len=*),parameter :: subname='(conservation_check_vthermo)'
    1873             : 
    1874             :       !----------------------------------------------------------------
    1875             :       ! If energy is not conserved, print diagnostics and exit.
    1876             :       !----------------------------------------------------------------
    1877             : 
    1878             :       !-----------------------------------------------------------------
    1879             :       ! Note that fsurf - flat = fsw + flw + fsens; i.e., the latent
    1880             :       ! heat is not included in the energy input, since (efinal - einit)
    1881             :       ! is the energy change in the system ice + vapor, and the latent
    1882             :       ! heat lost by the ice is equal to that gained by the vapor.
    1883             :       !-----------------------------------------------------------------
    1884             :        
    1885             :       einp = (fsurfn - flatn + fswint - fhocnn &
    1886   774714937 :            - fsnow*Lfresh - fadvocn) * dt
    1887   774714937 :       ferr = abs(efinal-einit-einp) / dt
    1888             : 
    1889   774714937 :       if (ferr > ferrmax) then
    1890           0 :          call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
    1891           0 :          call icepack_warnings_add(subname//" conservation_check_vthermo: Thermo energy conservation error" ) 
    1892             : 
    1893           0 :          write(warnstr,*) subname, 'Thermo energy conservation error'
    1894           0 :          call icepack_warnings_add(warnstr)
    1895           0 :          write(warnstr,*) subname, 'Flux error (W/m^2) =', ferr
    1896           0 :          call icepack_warnings_add(warnstr)
    1897           0 :          write(warnstr,*) subname, 'Energy error (J) =', ferr*dt
    1898           0 :          call icepack_warnings_add(warnstr)
    1899           0 :          write(warnstr,*) subname, 'Initial energy =', einit
    1900           0 :          call icepack_warnings_add(warnstr)
    1901           0 :          write(warnstr,*) subname, 'Final energy   =', efinal
    1902           0 :          call icepack_warnings_add(warnstr)
    1903           0 :          write(warnstr,*) subname, 'efinal - einit  =', efinal-einit
    1904           0 :          call icepack_warnings_add(warnstr)
    1905           0 :          write(warnstr,*) subname, 'fsurfn,flatn,fswint,fhocn, fsnow*Lfresh:'
    1906           0 :          call icepack_warnings_add(warnstr)
    1907           0 :          write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn, fsnow*Lfresh
    1908           0 :          call icepack_warnings_add(warnstr)
    1909           0 :          write(warnstr,*) subname, 'Input energy =', einp
    1910           0 :          call icepack_warnings_add(warnstr)
    1911           0 :          write(warnstr,*) subname, 'fbot,fcondbot:'
    1912           0 :          call icepack_warnings_add(warnstr)
    1913           0 :          write(warnstr,*) subname, fbot,fcondbotn
    1914           0 :          call icepack_warnings_add(warnstr)
    1915             : 
    1916             :          !         if (ktherm == 2) then
    1917           0 :          write(warnstr,*) subname, 'Intermediate energy =', einter
    1918           0 :          call icepack_warnings_add(warnstr)
    1919           0 :          write(warnstr,*) subname, 'efinal - einter =', &
    1920           0 :               efinal-einter
    1921           0 :          call icepack_warnings_add(warnstr)
    1922           0 :          write(warnstr,*) subname, 'einter - einit  =', &
    1923           0 :               einter-einit
    1924           0 :          call icepack_warnings_add(warnstr)
    1925           0 :          write(warnstr,*) subname, 'Conduction Error =', (einter-einit) &
    1926           0 :               - (fcondtopn*dt - fcondbotn*dt + fswint*dt)
    1927           0 :          call icepack_warnings_add(warnstr)
    1928           0 :          write(warnstr,*) subname, 'Melt/Growth Error =', (einter-einit) &
    1929           0 :               + ferr*dt - (fcondtopn*dt - fcondbotn*dt + fswint*dt)
    1930           0 :          call icepack_warnings_add(warnstr)
    1931           0 :          write(warnstr,*) subname, 'Advection Error =', fadvocn*dt
    1932           0 :          call icepack_warnings_add(warnstr)
    1933             :          !         endif
    1934             : 
    1935             :          !         write(warnstr,*) subname, fsurfn,flatn,fswint,fhocnn
    1936             :          !         call icepack_warnings_add(warnstr)
    1937             :          
    1938           0 :          write(warnstr,*) subname, 'dt*(fsurfn, flatn, fswint, fhocn, fsnow*Lfresh, fadvocn):'
    1939           0 :          call icepack_warnings_add(warnstr)
    1940           0 :          write(warnstr,*) subname, fsurfn*dt, flatn*dt, &
    1941           0 :               fswint*dt, fhocnn*dt, &
    1942           0 :               fsnow*Lfresh*dt, fadvocn*dt
    1943           0 :          call icepack_warnings_add(warnstr)
    1944           0 :          return
    1945             :       endif
    1946             : 
    1947             :       end subroutine conservation_check_vthermo
    1948             : 
    1949             : !=======================================================================
    1950             : !
    1951             : ! Given the vertical thermo state variables (hin, hsn),
    1952             : ! compute the new ice state variables (vicen, vsnon).
    1953             : ! Zero out state variables if ice has melted entirely.
    1954             : !
    1955             : ! authors William H. Lipscomb, LANL
    1956             : !         C. M. Bitz, UW
    1957             : !         Elizabeth C. Hunke, LANL
    1958             : 
    1959   774714937 :       subroutine update_state_vthermo(nilyr,    nslyr,    &
    1960             :                                       Tf,       Tsf,      &
    1961             :                                       hin,      hsn,      &
    1962   774714937 :                                       zqin,     zSin,     &
    1963   774714937 :                                       zqsn,               &
    1964             :                                       aicen,    vicen,    &
    1965             :                                       vsnon)
    1966             : 
    1967             :       integer (kind=int_kind), intent(in) :: &
    1968             :          nilyr , & ! number of ice layers
    1969             :          nslyr     ! number of snow layers
    1970             : 
    1971             :       real (kind=dbl_kind), intent(in) :: &
    1972             :          Tf              ! freezing temperature (C)
    1973             : 
    1974             :       real (kind=dbl_kind), intent(inout) :: &
    1975             :          Tsf             ! ice/snow surface temperature, Tsfcn
    1976             : 
    1977             :       real (kind=dbl_kind), intent(in) :: &
    1978             :          hin         , & ! ice thickness (m)
    1979             :          hsn             ! snow thickness (m)
    1980             : 
    1981             :       real (kind=dbl_kind), dimension (:), intent(inout) :: &
    1982             :          zqin        , & ! ice layer enthalpy (J m-3)
    1983             :          zSin        , & ! ice salinity    (ppt)
    1984             :          zqsn            ! snow layer enthalpy (J m-3)
    1985             : 
    1986             :       real (kind=dbl_kind), intent(inout) :: &
    1987             :          aicen       , & ! concentration of ice
    1988             :          vicen       , & ! volume per unit area of ice          (m)
    1989             :          vsnon           ! volume per unit area of snow         (m)
    1990             : 
    1991             :       ! local variables
    1992             : 
    1993             :       integer (kind=int_kind) :: &
    1994             :          k               ! ice layer index
    1995             : 
    1996             :       character(len=*),parameter :: subname='(update_state_vthermo)'
    1997             : 
    1998   774714937 :       if (hin <= c0) then
    1999      105574 :          aicen = c0
    2000      105574 :          vicen = c0
    2001      105574 :          vsnon = c0
    2002      105574 :          Tsf = Tf
    2003      844592 :          do k = 1, nilyr
    2004      844592 :             zqin(k) = c0
    2005             :          enddo
    2006      105574 :          if (ktherm == 2) then
    2007      844592 :             do k = 1, nilyr
    2008      844592 :                zSin(k) = c0
    2009             :             enddo
    2010             :          endif
    2011      211148 :          do k = 1, nslyr
    2012      211148 :             zqsn(k) = c0
    2013             :          enddo
    2014             :       else
    2015             :          ! aicen is already up to date
    2016   774609363 :          vicen = aicen * hin
    2017   774609363 :          vsnon = aicen * hsn
    2018             :       endif
    2019             :       
    2020   774714937 :       end subroutine update_state_vthermo
    2021             : 
    2022             : !=======================================================================
    2023             : !autodocument_start icepack_step_therm1
    2024             : ! Driver for thermodynamic changes not needed for coupling:
    2025             : ! transport in thickness space, lateral growth and melting.
    2026             : !
    2027             : ! authors: William H. Lipscomb, LANL
    2028             : !          Elizabeth C. Hunke, LANL
    2029             : 
    2030   722474472 :       subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
    2031   722474472 :                                     aicen_init  ,               &
    2032   722474472 :                                     vicen_init  , vsnon_init  , &
    2033   722474472 :                                     aice        , aicen       , &
    2034   722474472 :                                     vice        , vicen       , &
    2035   722474472 :                                     vsno        , vsnon       , &
    2036             :                                     uvel        , vvel        , &
    2037  1444948944 :                                     Tsfc        , zqsn        , &
    2038  1444948944 :                                     zqin        , zSin        , &
    2039   722474472 :                                     alvl        , vlvl        , &
    2040   722474472 :                                     apnd        , hpnd        , &
    2041   722474472 :                                     ipnd        ,               &
    2042   722474472 :                                     iage        , FY          , &
    2043   722474472 :                                     aerosno     , aeroice     , &
    2044   722474472 :                                     isosno      , isoice      , &
    2045             :                                     uatm        , vatm        , &
    2046             :                                     wind        , zlvl        , &
    2047             :                                     Qa          , rhoa        , &
    2048   722474472 :                                     Qa_iso      , &
    2049             :                                     Tair        , Tref        , &
    2050             :                                     Qref        , Uref        , &
    2051   722474472 :                                     Qref_iso    , &
    2052             :                                     Cdn_atm_ratio,              &
    2053             :                                     Cdn_ocn     , Cdn_ocn_skin, &
    2054             :                                     Cdn_ocn_floe, Cdn_ocn_keel, &
    2055             :                                     Cdn_atm     , Cdn_atm_skin, &
    2056             :                                     Cdn_atm_floe, Cdn_atm_pond, &
    2057             :                                     Cdn_atm_rdg , hfreebd     , &
    2058             :                                     hdraft      , hridge      , &
    2059             :                                     distrdg     , hkeel       , &
    2060             :                                     dkeel       , lfloe       , &
    2061             :                                     dfloe       ,               &
    2062             :                                     strax       , stray       , &
    2063             :                                     strairxT    , strairyT    , &
    2064             :                                     potT        , sst         , &
    2065             :                                     sss         , Tf          , &
    2066             :                                     strocnxT    , strocnyT    , &
    2067             :                                     fbot        ,               &
    2068             :                                     Tbot        , Tsnice      , &
    2069             :                                     frzmlt      , rside       , &
    2070             :                                     fside       ,               &
    2071             :                                     fsnow       , frain       , &
    2072             :                                     fpond       ,               &
    2073   722474472 :                                     fsurf       , fsurfn      , &
    2074   722474472 :                                     fcondtop    , fcondtopn   , &
    2075   722474472 :                                     fcondbot    , fcondbotn   , &
    2076   722474472 :                                     fswsfcn     , fswintn     , &
    2077   722474472 :                                     fswthrun    , fswabs      , &
    2078             :                                     flwout      ,               &
    2079   722474472 :                                     Sswabsn     , Iswabsn     , &
    2080             :                                     flw         , & 
    2081   722474472 :                                     fsens       , fsensn      , &
    2082   722474472 :                                     flat        , flatn       , &
    2083             :                                     evap        ,               &
    2084             :                                     evaps       , evapi       , &
    2085             :                                     fresh       , fsalt       , &
    2086             :                                     fhocn       , fswthru     , &
    2087   722474472 :                                     flatn_f     , fsensn_f    , &
    2088   722474472 :                                     fsurfn_f    , fcondtopn_f , &
    2089   722474472 :                                     faero_atm   , faero_ocn   , &
    2090   722474472 :                                     fiso_atm    , fiso_ocn    , &
    2091   722474472 :                                     fiso_evap   , &
    2092             :                                     HDO_ocn     , H2_16O_ocn  , &
    2093             :                                     H2_18O_ocn  ,  &
    2094   722474472 :                                     dhsn        , ffracn      , &
    2095   722474472 :                                     meltt       , melttn      , &
    2096   722474472 :                                     meltb       , meltbn      , &
    2097   722474472 :                                     melts       , meltsn      , &
    2098   722474472 :                                     congel      , congeln     , &
    2099   722474472 :                                     snoice      , snoicen     , &
    2100   722474472 :                                     dsnown      , &
    2101             :                                     lmask_n     , lmask_s     , &
    2102             :                                     mlt_onset   , frz_onset   , &
    2103             :                                     yday        , prescribed_ice)
    2104             : 
    2105             :       integer (kind=int_kind), intent(in) :: &
    2106             :          ncat    , & ! number of thickness categories
    2107             :          nilyr   , & ! number of ice layers
    2108             :          nslyr       ! number of snow layers
    2109             : 
    2110             :       real (kind=dbl_kind), intent(in) :: &
    2111             :          dt          , & ! time step
    2112             :          uvel        , & ! x-component of velocity (m/s)
    2113             :          vvel        , & ! y-component of velocity (m/s)
    2114             :          strax       , & ! wind stress components (N/m^2)
    2115             :          stray       , & ! 
    2116             :          yday            ! day of year
    2117             : 
    2118             :       logical (kind=log_kind), intent(in) :: &
    2119             :          lmask_n     , & ! northern hemisphere mask
    2120             :          lmask_s         ! southern hemisphere mask
    2121             : 
    2122             :       logical (kind=log_kind), intent(in), optional :: &
    2123             :          prescribed_ice  ! if .true., use prescribed ice instead of computed
    2124             : 
    2125             :       real (kind=dbl_kind), intent(inout) :: &
    2126             :          aice        , & ! sea ice concentration
    2127             :          vice        , & ! volume per unit area of ice          (m)
    2128             :          vsno        , & ! volume per unit area of snow         (m)
    2129             :          zlvl        , & ! atm level height (m)
    2130             :          uatm        , & ! wind velocity components (m/s)
    2131             :          vatm        , &
    2132             :          wind        , & ! wind speed (m/s)
    2133             :          potT        , & ! air potential temperature  (K)
    2134             :          Tair        , & ! air temperature  (K)
    2135             :          Qa          , & ! specific humidity (kg/kg)
    2136             :          rhoa        , & ! air density (kg/m^3)
    2137             :          frain       , & ! rainfall rate (kg/m^2 s)
    2138             :          fsnow       , & ! snowfall rate (kg/m^2 s)
    2139             :          fpond       , & ! fresh water flux to ponds (kg/m^2/s)
    2140             :          fresh       , & ! fresh water flux to ocean (kg/m^2/s)
    2141             :          fsalt       , & ! salt flux to ocean (kg/m^2/s)
    2142             :          fhocn       , & ! net heat flux to ocean (W/m^2)
    2143             :          fswthru     , & ! shortwave penetrating to ocean (W/m^2)
    2144             :          fsurf       , & ! net surface heat flux (excluding fcondtop)(W/m^2)
    2145             :          fcondtop    , & ! top surface conductive flux        (W/m^2)
    2146             :          fcondbot    , & ! bottom surface conductive flux     (W/m^2)
    2147             :          fsens       , & ! sensible heat flux (W/m^2)
    2148             :          flat        , & ! latent heat flux   (W/m^2)
    2149             :          fswabs      , & ! shortwave flux absorbed in ice and ocean (W/m^2)
    2150             :          flw         , & ! incoming longwave radiation (W/m^2)
    2151             :          flwout      , & ! outgoing longwave radiation (W/m^2)
    2152             :          evap        , & ! evaporative water flux (kg/m^2/s)
    2153             :          evaps       , & ! evaporative water flux over snow (kg/m^2/s)
    2154             :          evapi       , & ! evaporative water flux over ice (kg/m^2/s)
    2155             :          congel      , & ! basal ice growth         (m/step-->cm/day)
    2156             :          snoice      , & ! snow-ice formation       (m/step-->cm/day)
    2157             :          Tref        , & ! 2m atm reference temperature (K)
    2158             :          Qref        , & ! 2m atm reference spec humidity (kg/kg)
    2159             :          Uref        , & ! 10m atm reference wind speed (m/s)
    2160             :          Cdn_atm     , & ! atm drag coefficient
    2161             :          Cdn_ocn     , & ! ocn drag coefficient
    2162             :          hfreebd     , & ! freeboard (m)
    2163             :          hdraft      , & ! draft of ice + snow column (Stoessel1993)
    2164             :          hridge      , & ! ridge height
    2165             :          distrdg     , & ! distance between ridges
    2166             :          hkeel       , & ! keel depth
    2167             :          dkeel       , & ! distance between keels
    2168             :          lfloe       , & ! floe length
    2169             :          dfloe       , & ! distance between floes
    2170             :          Cdn_atm_skin, & ! neutral skin drag coefficient
    2171             :          Cdn_atm_floe, & ! neutral floe edge drag coefficient
    2172             :          Cdn_atm_pond, & ! neutral pond edge drag coefficient
    2173             :          Cdn_atm_rdg , & ! neutral ridge drag coefficient
    2174             :          Cdn_ocn_skin, & ! skin drag coefficient
    2175             :          Cdn_ocn_floe, & ! floe edge drag coefficient
    2176             :          Cdn_ocn_keel, & ! keel drag coefficient
    2177             :          Cdn_atm_ratio,& ! ratio drag atm / neutral drag atm
    2178             :          strairxT    , & ! stress on ice by air, x-direction
    2179             :          strairyT    , & ! stress on ice by air, y-direction
    2180             :          strocnxT    , & ! ice-ocean stress, x-direction
    2181             :          strocnyT    , & ! ice-ocean stress, y-direction
    2182             :          fbot        , & ! ice-ocean heat flux at bottom surface (W/m^2)
    2183             :          frzmlt      , & ! freezing/melting potential (W/m^2)
    2184             :          rside       , & ! fraction of ice that melts laterally
    2185             :          fside       , & ! lateral heat flux (W/m^2)
    2186             :          sst         , & ! sea surface temperature (C)
    2187             :          Tf          , & ! freezing temperature (C)
    2188             :          Tbot        , & ! ice bottom surface temperature (deg C)
    2189             :          Tsnice      , & ! snow ice interface temperature (deg C)
    2190             :          sss         , & ! sea surface salinity (ppt)
    2191             :          meltt       , & ! top ice melt             (m/step-->cm/day)
    2192             :          melts       , & ! snow melt                (m/step-->cm/day)
    2193             :          meltb       , & ! basal ice melt           (m/step-->cm/day)
    2194             :          mlt_onset   , & ! day of year that sfc melting begins
    2195             :          frz_onset       ! day of year that freezing begins (congel or frazil)
    2196             : 
    2197             :       real (kind=dbl_kind), dimension(:), optional, intent(inout) :: &
    2198             :          Qa_iso      , & ! isotope specific humidity (kg/kg)
    2199             :          Qref_iso    , & ! isotope 2m atm reference spec humidity (kg/kg)
    2200             :          fiso_atm    , & ! isotope deposition rate (kg/m^2 s)
    2201             :          fiso_ocn    , & ! isotope flux to ocean  (kg/m^2/s)
    2202             :          fiso_evap       ! isotope evaporation (kg/m^2/s)
    2203             : 
    2204             :       real (kind=dbl_kind), optional, intent(in) :: &
    2205             :          HDO_ocn     , & ! ocean concentration of HDO (kg/kg)
    2206             :          H2_16O_ocn  , & ! ocean concentration of H2_16O (kg/kg)
    2207             :          H2_18O_ocn      ! ocean concentration of H2_18O (kg/kg)
    2208             : 
    2209             :       real (kind=dbl_kind), dimension(:), intent(inout) :: &
    2210             :          aicen_init  , & ! fractional area of ice
    2211             :          vicen_init  , & ! volume per unit area of ice (m)
    2212             :          vsnon_init  , & ! volume per unit area of snow (m)
    2213             :          aicen       , & ! concentration of ice
    2214             :          vicen       , & ! volume per unit area of ice          (m)
    2215             :          vsnon       , & ! volume per unit area of snow         (m)
    2216             :          Tsfc        , & ! ice/snow surface temperature, Tsfcn
    2217             :          alvl        , & ! level ice area fraction
    2218             :          vlvl        , & ! level ice volume fraction
    2219             :          apnd        , & ! melt pond area fraction
    2220             :          hpnd        , & ! melt pond depth (m)
    2221             :          ipnd        , & ! melt pond refrozen lid thickness (m)
    2222             :          iage        , & ! volume-weighted ice age
    2223             :          FY          , & ! area-weighted first-year ice area
    2224             :          fsurfn      , & ! net flux to top surface, excluding fcondtop
    2225             :          fcondtopn   , & ! downward cond flux at top surface (W m-2)
    2226             :          fcondbotn   , & ! downward cond flux at bottom surface (W m-2)
    2227             :          flatn       , & ! latent heat flux (W m-2)
    2228             :          fsensn      , & ! sensible heat flux (W m-2)
    2229             :          fsurfn_f    , & ! net flux to top surface, excluding fcondtop
    2230             :          fcondtopn_f , & ! downward cond flux at top surface (W m-2)
    2231             :          flatn_f     , & ! latent heat flux (W m-2)
    2232             :          fsensn_f    , & ! sensible heat flux (W m-2)
    2233             :          fswsfcn     , & ! SW absorbed at ice/snow surface (W m-2)
    2234             :          fswthrun    , & ! SW through ice to ocean            (W/m^2)
    2235             :          fswintn     , & ! SW absorbed in ice interior, below surface (W m-2)
    2236             :          faero_atm   , & ! aerosol deposition rate (kg/m^2 s)
    2237             :          faero_ocn   , & ! aerosol flux to ocean  (kg/m^2/s)
    2238             :          dhsn        , & ! depth difference for snow on sea ice and pond ice
    2239             :          ffracn      , & ! fraction of fsurfn used to melt ipond
    2240             :          meltsn      , & ! snow melt                       (m)
    2241             :          melttn      , & ! top ice melt                    (m)
    2242             :          meltbn      , & ! bottom ice melt                 (m)
    2243             :          congeln     , & ! congelation ice growth          (m)
    2244             :          snoicen     , & ! snow-ice growth                 (m)
    2245             :          dsnown          ! change in snow thickness (m/step-->cm/day)
    2246             : 
    2247             :       real (kind=dbl_kind), dimension(:,:), intent(inout) :: &
    2248             :          zqsn        , & ! snow layer enthalpy (J m-3)
    2249             :          zqin        , & ! ice layer enthalpy (J m-3)
    2250             :          zSin        , & ! internal ice layer salinities
    2251             :          Sswabsn     , & ! SW radiation absorbed in snow layers (W m-2)
    2252             :          Iswabsn         ! SW radiation absorbed in ice layers (W m-2)
    2253             : 
    2254             :       real (kind=dbl_kind), dimension(:,:,:), intent(inout) :: &
    2255             :          aerosno    , &  ! snow aerosol tracer (kg/m^2)
    2256             :          aeroice         ! ice aerosol tracer (kg/m^2)
    2257             : 
    2258             :       real (kind=dbl_kind), dimension(:,:), optional, intent(inout) :: &
    2259             :          isosno     , &  ! snow isotope tracer (kg/m^2)
    2260             :          isoice          ! ice isotope tracer (kg/m^2)
    2261             : !autodocument_end
    2262             : 
    2263             :       ! local variables
    2264             : 
    2265             :       integer (kind=int_kind) :: &
    2266             :          n               ! category index
    2267             : 
    2268             :       real (kind=dbl_kind) :: &
    2269    45106800 :          worka, workb    ! temporary variables
    2270             : 
    2271             :       ! 2D coupler variables (computed for each category, then aggregated)
    2272             :       real (kind=dbl_kind) :: &
    2273    45106800 :          fswabsn     , & ! shortwave absorbed by ice          (W/m^2)
    2274    45106800 :          flwoutn     , & ! upward LW at surface               (W/m^2)
    2275    45106800 :          evapn       , & ! flux of vapor, atmos to ice   (kg m-2 s-1)
    2276    45106800 :          evapsn      , & ! flux of vapor, atmos to ice over snow  (kg m-2 s-1)
    2277    45106800 :          evapin      , & ! flux of vapor, atmos to ice over ice  (kg m-2 s-1)
    2278    45106800 :          freshn      , & ! flux of water, ice to ocean     (kg/m^2/s)
    2279    45106800 :          fsaltn      , & ! flux of salt, ice to ocean      (kg/m^2/s)
    2280    45106800 :          fhocnn      , & ! fbot corrected for leftover energy (W/m^2)
    2281    45106800 :          strairxn    , & ! air/ice zonal  stress,             (N/m^2)
    2282    45106800 :          strairyn    , & ! air/ice meridional stress,         (N/m^2)
    2283    45106800 :          Cdn_atm_ratio_n, & ! drag coefficient ratio
    2284    45106800 :          Trefn       , & ! air tmp reference level                (K)
    2285    45106800 :          Urefn       , & ! air speed reference level            (m/s)
    2286    45106800 :          Qrefn       , & ! air sp hum reference level         (kg/kg)
    2287    45106800 :          shcoef      , & ! transfer coefficient for sensible heat
    2288    45106800 :          lhcoef      , & ! transfer coefficient for latent heat
    2289    45106800 :          rfrac           ! water fraction retained for melt ponds
    2290             : 
    2291             :       real (kind=dbl_kind), dimension(n_iso) :: &
    2292  1400418576 :          Qrefn_iso  , & ! isotope air sp hum reference level         (kg/kg)
    2293  1400418576 :          fiso_ocnn  , & ! isotope flux to ocean  (kg/m^2/s)
    2294  1490632176 :          fiso_evapn     ! isotope evaporation (kg/m^2/s)
    2295             : 
    2296             :       real (kind=dbl_kind), allocatable, dimension(:,:) :: &
    2297   722474472 :          l_isosno   , &  ! local snow isotope tracer (kg/m^2)
    2298   722474472 :          l_isoice        ! local ice isotope tracer (kg/m^2)
    2299             : 
    2300             :       real (kind=dbl_kind), allocatable, dimension(:) :: &
    2301   722474472 :          l_Qa_iso    , & ! local isotope specific humidity (kg/kg)
    2302   722474472 :          l_Qref_iso  , & ! local isotope 2m atm reference spec humidity (kg/kg)
    2303   722474472 :          l_fiso_atm  , & ! local isotope deposition rate (kg/m^2 s)
    2304   722474472 :          l_fiso_ocn  , & ! local isotope flux to ocean  (kg/m^2/s)
    2305   722474472 :          l_fiso_evap     ! local isotope evaporation (kg/m^2/s)
    2306             : 
    2307             :       real (kind=dbl_kind)  :: &
    2308    45106800 :          l_HDO_ocn   , & ! local ocean concentration of HDO (kg/kg)
    2309    45106800 :          l_H2_16O_ocn, & ! local ocean concentration of H2_16O (kg/kg)
    2310    45106800 :          l_H2_18O_ocn    ! local ocean concentration of H2_18O (kg/kg)
    2311             : 
    2312             :       real (kind=dbl_kind) :: &
    2313    45106800 :          pond            ! water retained in ponds (m)
    2314             : 
    2315             :       character(len=*),parameter :: subname='(icepack_step_therm1)'
    2316             : 
    2317             :       !-----------------------------------------------------------------
    2318             :       ! allocate local optional arguments
    2319             :       !-----------------------------------------------------------------
    2320             : 
    2321   722474472 :       if (present(isosno)    ) then
    2322   722474472 :          allocate(l_isosno(size(isosno,dim=1),size(isosno,dim=2)))
    2323  4571183952 :          l_isosno     = isosno
    2324             :       else
    2325           0 :          allocate(l_isosno(1,1))
    2326           0 :          l_isosno     = c0
    2327             :       endif
    2328             : 
    2329   722474472 :       if (present(isoice)    ) then
    2330   722474472 :          allocate(l_isoice(size(isoice,dim=1),size(isoice,dim=2)))
    2331  4571183952 :          l_isoice     = isoice
    2332             :       else
    2333           0 :          allocate(l_isoice(1,1))
    2334           0 :          l_isoice     = c0
    2335             :       endif
    2336             : 
    2337   722474472 :       if (present(Qa_iso)    ) then
    2338   722474472 :          allocate(l_Qa_iso(size(Qa_iso)))
    2339  2889897888 :          l_Qa_iso     = Qa_iso
    2340             :       else
    2341           0 :          allocate(l_Qa_iso(1))
    2342           0 :          l_Qa_iso     = c0
    2343             :       endif
    2344             : 
    2345   722474472 :       if (present(Qref_iso)    ) then
    2346   722474472 :          allocate(l_Qref_iso(size(Qref_iso)))
    2347  2889897888 :          l_Qref_iso     = Qref_iso
    2348             :       else
    2349           0 :          allocate(l_Qref_iso(1))
    2350           0 :          l_Qref_iso     = c0
    2351             :       endif
    2352             : 
    2353   722474472 :       if (present(fiso_atm)  ) then
    2354   722474472 :          allocate(l_fiso_atm(size(fiso_atm)))
    2355  2889897888 :          l_fiso_atm = fiso_atm
    2356             :       else
    2357           0 :          allocate(l_fiso_atm(1))
    2358           0 :          l_fiso_atm   = c0
    2359             :       endif
    2360             : 
    2361   722474472 :       if (present(fiso_ocn)  ) then
    2362   722474472 :          allocate(l_fiso_ocn(size(fiso_ocn)))
    2363  2889897888 :          l_fiso_ocn = fiso_ocn
    2364             :       else
    2365           0 :          allocate(l_fiso_ocn(1))
    2366           0 :          l_fiso_ocn   = c0
    2367             :       endif
    2368             : 
    2369   722474472 :       if (present(fiso_evap)  ) then
    2370   722474472 :          allocate(l_fiso_evap(size(fiso_evap)))
    2371  2889897888 :          l_fiso_evap = fiso_evap
    2372             :       else
    2373           0 :          allocate(l_fiso_evap(1))
    2374           0 :          l_fiso_evap   = c0
    2375             :       endif
    2376             : 
    2377   722474472 :       l_HDO_ocn    = c0
    2378   722474472 :       if (present(HDO_ocn)   ) l_HDO_ocn    = HDO_ocn
    2379             : 
    2380   722474472 :       l_H2_16O_ocn = c0
    2381   722474472 :       if (present(H2_16O_ocn)) l_H2_16O_ocn = H2_16O_ocn
    2382             : 
    2383   722474472 :       l_H2_18O_ocn = c0
    2384   722474472 :       if (present(H2_18O_ocn)) l_H2_18O_ocn = H2_18O_ocn
    2385             : 
    2386             :       !-----------------------------------------------------------------
    2387             :       ! Adjust frzmlt to account for ice-ocean heat fluxes since last
    2388             :       !  call to coupler.
    2389             :       ! Compute lateral and bottom heat fluxes.
    2390             :       !-----------------------------------------------------------------
    2391             : 
    2392             :       call frzmlt_bottom_lateral (dt,        ncat,      &
    2393             :                                   nilyr,     nslyr,     &
    2394             :                                   aice,      frzmlt,    &
    2395           0 :                                   vicen,     vsnon,     &
    2396           0 :                                   zqin,      zqsn,      &
    2397             :                                   sst,       Tf,        &
    2398             :                                   ustar_min,            &
    2399             :                                   fbot_xfer_type,       &
    2400             :                                   strocnxT,  strocnyT,  &
    2401             :                                   Tbot,      fbot,      &
    2402             :                                   rside,     Cdn_ocn,   &
    2403   722474472 :                                   fside)
    2404             : 
    2405   722474472 :       if (icepack_warnings_aborted(subname)) return
    2406             : 
    2407             :       !-----------------------------------------------------------------
    2408             :       ! Update the neutral drag coefficients to account for form drag
    2409             :       ! Oceanic and atmospheric drag coefficients
    2410             :       !-----------------------------------------------------------------
    2411             : 
    2412   722474472 :       if (formdrag) then
    2413           0 :          call neutral_drag_coeffs (apnd         , &
    2414           0 :                                    hpnd        , ipnd         , &
    2415           0 :                                    alvl        , vlvl         , &
    2416             :                                    aice        , vice,          &
    2417           0 :                                    vsno        , aicen        , &
    2418           0 :                                    vicen       , &
    2419             :                                    Cdn_ocn     , Cdn_ocn_skin, &
    2420             :                                    Cdn_ocn_floe, Cdn_ocn_keel, &
    2421             :                                    Cdn_atm     , Cdn_atm_skin, &
    2422             :                                    Cdn_atm_floe, Cdn_atm_pond, &
    2423             :                                    Cdn_atm_rdg , hfreebd     , &
    2424             :                                    hdraft      , hridge      , &
    2425             :                                    distrdg     , hkeel       , &
    2426             :                                    dkeel       , lfloe       , &
    2427    20367264 :                                    dfloe       , ncat)
    2428    20367264 :          if (icepack_warnings_aborted(subname)) return
    2429             :       endif
    2430             : 
    2431  4265674992 :       do n = 1, ncat
    2432             : 
    2433  3543200520 :          meltsn (n) = c0
    2434  3543200520 :          melttn (n) = c0
    2435  3543200520 :          meltbn (n) = c0
    2436  3543200520 :          congeln(n) = c0
    2437  3543200520 :          snoicen(n) = c0
    2438  3543200520 :          dsnown (n) = c0
    2439             : 
    2440  3543200520 :          Trefn  = c0
    2441  3543200520 :          Qrefn  = c0
    2442  3848709480 :          Qrefn_iso(:) = c0
    2443  3848709480 :          fiso_ocnn(:) = c0
    2444  3848709480 :          fiso_evapn(:) = c0
    2445  3543200520 :          Urefn  = c0
    2446  3543200520 :          lhcoef = c0
    2447  3543200520 :          shcoef = c0
    2448  3543200520 :          worka  = c0
    2449  3543200520 :          workb  = c0
    2450             : 
    2451  3543200520 :          if (aicen_init(n) > puny) then
    2452             : 
    2453   774714937 :             if (calc_Tsfc .or. calc_strair) then 
    2454             : 
    2455             :       !-----------------------------------------------------------------
    2456             :       ! Atmosphere boundary layer calculation; compute coefficients
    2457             :       ! for sensible and latent heat fluxes.
    2458             :       !
    2459             :       ! NOTE: The wind stress is computed here for later use if 
    2460             :       !       calc_strair = .true.   Otherwise, the wind stress
    2461             :       !       components are set to the data values.
    2462             :       !-----------------------------------------------------------------
    2463             : 
    2464             :                call icepack_atm_boundary('ice',                  &
    2465    54309860 :                                         Tsfc(n),  potT,          &
    2466             :                                         uatm,     vatm,          &
    2467             :                                         wind,     zlvl,          &
    2468             :                                         Qa,       rhoa,          &
    2469             :                                         strairxn, strairyn,      &
    2470             :                                         Trefn,    Qrefn,         &
    2471             :                                         worka,    workb,         &
    2472             :                                         lhcoef,   shcoef,        &
    2473             :                                         Cdn_atm,                 &
    2474             :                                         Cdn_atm_ratio_n,         &
    2475             :                                         Qa_iso=l_Qa_iso,           &
    2476           0 :                                         Qref_iso=Qrefn_iso,      &
    2477             :                                         uvel=uvel, vvel=vvel,    &
    2478   774714937 :                                         Uref=Urefn)
    2479   774714937 :                if (icepack_warnings_aborted(subname)) return
    2480             : 
    2481             :             endif   ! calc_Tsfc or calc_strair
    2482             : 
    2483   774714937 :             if (.not.(calc_strair)) then
    2484             : #ifndef CICE_IN_NEMO
    2485             :                ! Set to data values (on T points)
    2486           0 :                strairxn = strax
    2487           0 :                strairyn = stray
    2488             : #else
    2489             :                ! NEMO wind stress is supplied on u grid, multipied 
    2490             :                ! by ice concentration and set directly in evp, so
    2491             :                ! strairxT/yT = 0. Zero u-components here for safety.
    2492             :                strairxn = c0
    2493             :                strairyn = c0
    2494             : #endif
    2495             :             endif
    2496             : 
    2497             :       !-----------------------------------------------------------------
    2498             :       ! Update ice age
    2499             :       ! This is further adjusted for freezing in the thermodynamics.
    2500             :       ! Melting does not alter the ice age.
    2501             :       !-----------------------------------------------------------------
    2502             : 
    2503   774714937 :             if (tr_iage) call increment_age (dt, iage(n))
    2504   774714937 :             if (icepack_warnings_aborted(subname)) return
    2505   774714937 :             if (tr_FY)   call update_FYarea (dt,               &
    2506             :                                              lmask_n, lmask_s, &
    2507   634216772 :                                              yday,    FY(n))
    2508   774714937 :             if (icepack_warnings_aborted(subname)) return
    2509             : 
    2510             :       !-----------------------------------------------------------------
    2511             :       ! Vertical thermodynamics: Heat conduction, growth and melting.
    2512             :       !----------------------------------------------------------------- 
    2513             : 
    2514   774714937 :             if (.not.(calc_Tsfc)) then
    2515             : 
    2516             :                ! If not calculating surface temperature and fluxes, set 
    2517             :                ! surface fluxes (flatn, fsurfn, and fcondtopn) to be used 
    2518             :                ! in thickness_changes
    2519             :  
    2520             :                ! hadgem routine sets fluxes to default values in ice-only mode
    2521     3773267 :                call set_sfcflux(aicen      (n),                 &
    2522     3773267 :                                 flatn_f    (n), fsensn_f   (n), &
    2523     3773267 :                                 fcondtopn_f(n),                 &
    2524     3773267 :                                 fsurfn_f   (n),                 &
    2525     3773267 :                                 flatn      (n), fsensn     (n), &
    2526     3773267 :                                 fsurfn     (n),                 &
    2527    29876420 :                                 fcondtopn  (n))
    2528    29876420 :                if (icepack_warnings_aborted(subname)) return
    2529             :             endif
    2530             : 
    2531             :             call thermo_vertical(nilyr,        nslyr,        &
    2532    54309860 :                                  dt,           aicen    (n), &
    2533    54309860 :                                  vicen    (n), vsnon    (n), &
    2534    54309860 :                                  Tsfc     (n), zSin   (:,n), &
    2535           0 :                                  zqin   (:,n), zqsn   (:,n), &
    2536    54309860 :                                  apnd     (n), hpnd     (n), &
    2537             :                                  tr_pond_topo, &
    2538             :                                  flw,          potT,         &
    2539             :                                  Qa,           rhoa,         &
    2540             :                                  fsnow,        fpond,        &
    2541             :                                  fbot,         Tbot,         &
    2542             :                                  Tsnice,        sss,          &
    2543             :                                  lhcoef,       shcoef,       &
    2544    54309860 :                                  fswsfcn  (n), fswintn  (n), &
    2545           0 :                                  Sswabsn(:,n), Iswabsn(:,n), &
    2546    54309860 :                                  fsurfn   (n), fcondtopn(n), &
    2547    54309860 :                                  fcondbotn(n),               &
    2548    54309860 :                                  fsensn   (n), flatn    (n), &
    2549             :                                  flwoutn,      evapn,        &
    2550             :                                  evapsn,       evapin,       &
    2551             :                                  freshn,       fsaltn,       &
    2552             :                                  fhocnn,                     &
    2553    54309860 :                                  melttn   (n), meltsn   (n), &
    2554    54309860 :                                  meltbn   (n),               &
    2555    54309860 :                                  congeln  (n), snoicen  (n), &
    2556             :                                  mlt_onset,    frz_onset,    &
    2557    54309860 :                                  yday,         dsnown   (n), &
    2558   883334657 :                                  prescribed_ice)
    2559             : 
    2560   774714937 :             if (icepack_warnings_aborted(subname)) then
    2561           0 :                call icepack_warnings_add(subname//' ice: Vertical thermo error: ')
    2562           0 :                return
    2563             :             endif
    2564             : 
    2565             :       !-----------------------------------------------------------------
    2566             :       ! Total absorbed shortwave radiation
    2567             :       !-----------------------------------------------------------------
    2568             : 
    2569   774714937 :             fswabsn = fswsfcn(n) + fswintn(n) + fswthrun(n)
    2570             : 
    2571             :       !-----------------------------------------------------------------
    2572             :       ! Aerosol update
    2573             :       !-----------------------------------------------------------------
    2574             : 
    2575   774714937 :             if (tr_aero) then
    2576             :                call update_aerosol (dt,                             &
    2577             :                                     nilyr, nslyr, n_aero,           &
    2578     3952149 :                                     melttn     (n), meltsn     (n), &
    2579     3952149 :                                     meltbn     (n), congeln    (n), &
    2580     3952149 :                                     snoicen    (n), fsnow,          &
    2581           0 :                                     aerosno(:,:,n), aeroice(:,:,n), &
    2582     3952149 :                                     aicen_init (n), vicen_init (n), &
    2583     3952149 :                                     vsnon_init (n),                 &
    2584     3952149 :                                     vicen      (n), vsnon      (n), &
    2585     3952149 :                                     aicen      (n),                 &
    2586    56061574 :                                     faero_atm     ,  faero_ocn)
    2587    52109425 :                if (icepack_warnings_aborted(subname)) return
    2588             :             endif
    2589             : 
    2590   774714937 :             if (tr_iso) then
    2591             :                call update_isotope (dt = dt, &
    2592             :                                     nilyr = nilyr, nslyr = nslyr, &
    2593      205894 :                                     meltt = melttn(n),melts = meltsn(n),     &
    2594      205894 :                                     meltb = meltbn(n),congel=congeln(n),    &
    2595      205894 :                                     snoice=snoicen(n),evap=evapn,         & 
    2596      205894 :                                     fsnow=fsnow,      Tsfc=Tsfc(n),       &
    2597           0 :                                     Qref_iso=Qrefn_iso(:),                 &
    2598           0 :                                     isosno=l_isosno(:,n),isoice=l_isoice(:,n), &
    2599      205894 :                                     aice_old=aicen_init(n),vice_old=vicen_init(n), &
    2600      205894 :                                     vsno_old=vsnon_init(n),                &
    2601      205894 :                                     vicen=vicen(n),vsnon=vsnon(n),      &
    2602      205894 :                                     aicen=aicen(n),                     &
    2603             :                                     fiso_atm=l_fiso_atm(:),                  &
    2604           0 :                                     fiso_evapn=fiso_evapn(:),                &
    2605           0 :                                     fiso_ocnn=fiso_ocnn(:),                 &
    2606             :                                     HDO_ocn=l_HDO_ocn,H2_16O_ocn=l_H2_16O_ocn,    &
    2607    22164656 :                                     H2_18O_ocn=l_H2_18O_ocn)
    2608    21958762 :                if (icepack_warnings_aborted(subname)) return
    2609             :             endif
    2610             :          endif   ! aicen_init
    2611             : 
    2612             :       !-----------------------------------------------------------------
    2613             :       ! Melt ponds
    2614             :       ! If using tr_pond_cesm, the full calculation is performed here.
    2615             :       ! If using tr_pond_topo, the rest of the calculation is done after
    2616             :       ! the surface fluxes are merged, below.
    2617             :       !-----------------------------------------------------------------
    2618             : 
    2619             :          !call ice_timer_start(timer_ponds)
    2620  3543200520 :          if (tr_pond) then
    2621             :                
    2622  3286668360 :             if (tr_pond_cesm) then
    2623   101836320 :                rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n) 
    2624             :                call compute_ponds_cesm(dt,        hi_min,    &
    2625             :                                        pndaspect, rfrac,     &
    2626      960720 :                                        melttn(n), meltsn(n), &
    2627             :                                        frain,                &
    2628      960720 :                                        aicen (n), vicen (n), &
    2629      960720 :                                        Tsfc  (n), &
    2630   101836320 :                                        apnd  (n), hpnd  (n))
    2631   101836320 :                if (icepack_warnings_aborted(subname)) return
    2632             :                   
    2633  3184832040 :             elseif (tr_pond_lvl) then
    2634  3046488360 :                rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
    2635             :                call compute_ponds_lvl(dt,        nilyr,     &
    2636             :                                       ktherm,               &
    2637             :                                       hi_min,               &
    2638             :                                       dpscale,   frzpnd,    &
    2639             :                                       pndaspect, rfrac,     &
    2640   148911600 :                                       melttn(n), meltsn(n), &
    2641             :                                       frain,     Tair,      &
    2642   148911600 :                                       fsurfn(n),            &
    2643   148911600 :                                       dhsn  (n), ffracn(n), &
    2644   148911600 :                                       aicen (n), vicen (n), &
    2645   148911600 :                                       vsnon (n),            &
    2646           0 :                                       zqin(:,n), zSin(:,n), &
    2647   148911600 :                                       Tsfc  (n), alvl  (n), &
    2648   148911600 :                                       apnd  (n), hpnd  (n), &
    2649  3195399960 :                                       ipnd  (n))
    2650  3046488360 :                if (icepack_warnings_aborted(subname)) return
    2651             :                   
    2652   138343680 :             elseif (tr_pond_topo) then
    2653   138343680 :                if (aicen_init(n) > puny) then
    2654             :                      
    2655             :                   ! collect liquid water in ponds
    2656             :                   ! assume salt still runs off
    2657    29876420 :                   rfrac = rfracmin + (rfracmax-rfracmin) * aicen(n)
    2658     3773267 :                   pond = rfrac/rhofresh * (melttn(n)*rhoi &
    2659     3773267 :                        +                   meltsn(n)*rhos &
    2660    29876420 :                        +                   frain *dt)
    2661             : 
    2662             :                   ! if pond does not exist, create new pond over full ice area
    2663             :                   ! otherwise increase pond depth without changing pond area
    2664    29876420 :                   if (apnd(n) < puny) then
    2665    29864249 :                      hpnd(n) = c0
    2666    29864249 :                      apnd(n) = c1
    2667             :                   endif
    2668    29876420 :                   hpnd(n) = (pond + hpnd(n)*apnd(n)) / apnd(n)
    2669    29876420 :                   fpond = fpond + pond * aicen(n) ! m
    2670             :                endif ! aicen_init
    2671             :             endif
    2672             : 
    2673             :          endif ! tr_pond
    2674             :          !call ice_timer_stop(timer_ponds)
    2675             : 
    2676             :       !-----------------------------------------------------------------
    2677             :       ! Increment area-weighted fluxes.
    2678             :       !-----------------------------------------------------------------
    2679             : 
    2680  3543200520 :          if (aicen_init(n) > puny) &
    2681    54309860 :             call merge_fluxes (aicen_init(n),            &
    2682             :                                flw, & 
    2683             :                                strairxn,   strairyn,     &
    2684             :                                Cdn_atm_ratio_n,          &
    2685    54309860 :                                fsurfn(n),  fcondtopn(n), &
    2686    54309860 :                                fcondbotn(n),             &
    2687    54309860 :                                fsensn(n),  flatn(n),     &
    2688             :                                fswabsn,    flwoutn,      &
    2689             :                                evapn,                    &
    2690             :                                evapsn,     evapin,       &
    2691             :                                Trefn,      Qrefn,        &
    2692             :                                freshn,     fsaltn,       &
    2693    54309860 :                                fhocnn,     fswthrun(n),  &
    2694             :                                strairxT,   strairyT,     &
    2695             :                                Cdn_atm_ratio,            &
    2696             :                                fsurf,      fcondtop,     &
    2697             :                                fcondbot,                 &
    2698             :                                fsens,      flat,         &
    2699             :                                fswabs,     flwout,       &
    2700             :                                evap,                     &
    2701             :                                evaps,      evapi,        &
    2702             :                                Tref,       Qref,         &
    2703             :                                fresh,      fsalt,        &
    2704             :                                fhocn,      fswthru,      &
    2705    54309860 :                                melttn (n), meltsn(n),    &
    2706    54309860 :                                meltbn (n), congeln(n),   &
    2707    54309860 :                                snoicen(n),               &
    2708             :                                meltt,      melts,        &
    2709             :                                meltb,      congel,       &
    2710             :                                snoice,                   &
    2711             :                                Uref=Uref,  Urefn=Urefn,  &
    2712             :                                Qref_iso=l_Qref_iso,      &
    2713           0 :                                Qrefn_iso=Qrefn_iso,      &
    2714             :                                fiso_ocn=l_fiso_ocn,      &
    2715           0 :                                fiso_ocnn=fiso_ocnn,      &
    2716             :                                fiso_evap=l_fiso_evap,    &
    2717   774714937 :                                fiso_evapn=fiso_evapn)
    2718             : 
    2719  4265674992 :          if (icepack_warnings_aborted(subname)) return
    2720             : 
    2721             :       enddo                  ! ncat
    2722             : 
    2723  4571183952 :       if (present(isosno)   ) isosno   = l_isosno
    2724  4571183952 :       if (present(isoice)   ) isoice   = l_isoice
    2725  2889897888 :       if (present(Qa_iso)   ) Qa_iso   = l_Qa_iso
    2726  2889897888 :       if (present(Qref_iso) ) Qref_iso = l_Qref_iso
    2727  2889897888 :       if (present(fiso_atm) ) fiso_atm = l_fiso_atm
    2728  2889897888 :       if (present(fiso_ocn) ) fiso_ocn = l_fiso_ocn
    2729  2889897888 :       if (present(fiso_evap)) fiso_evap= l_fiso_evap
    2730   722474472 :       deallocate(l_isosno)
    2731   722474472 :       deallocate(l_isoice)
    2732   722474472 :       deallocate(l_Qa_iso)
    2733   722474472 :       deallocate(l_Qref_iso)
    2734   722474472 :       deallocate(l_fiso_atm)
    2735   722474472 :       deallocate(l_fiso_ocn)
    2736   722474472 :       deallocate(l_fiso_evap)
    2737             : 
    2738             :       !-----------------------------------------------------------------
    2739             :       ! Calculate ponds from the topographic scheme
    2740             :       !-----------------------------------------------------------------
    2741             :       !call ice_timer_start(timer_ponds)
    2742   722474472 :       if (tr_pond_topo) then
    2743             :          call compute_ponds_topo(dt,       ncat,      nilyr,     &
    2744             :                                  ktherm,   heat_capacity,        &
    2745           0 :                                  aice,     aicen,                &
    2746           0 :                                  vice,     vicen,                &
    2747           0 :                                  vsno,     vsnon,                &
    2748             :                                  meltt,                &
    2749             :                                  fsurf,    fpond,                &
    2750           0 :                                  Tsfc,     Tf,                   &
    2751           0 :                                  zqin,     zSin,                 &
    2752    23057280 :                                  apnd,     hpnd,      ipnd       )
    2753    23057280 :          if (icepack_warnings_aborted(subname)) return
    2754             :       endif
    2755             :       !call ice_timer_stop(timer_ponds)
    2756             : 
    2757   722474472 :       end subroutine icepack_step_therm1
    2758             : 
    2759             : !=======================================================================
    2760             : 
    2761             :       end module icepack_therm_vertical
    2762             : 
    2763             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd