From cd863c45878f95d79ec678a7b00b574de3f3a07b Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Tue, 14 Sep 2021 18:59:39 +0000
Subject: [PATCH 01/20] 'heatflux_linear' flag: sensible+latent heatfluxes
 using traditional linear bulk formula

---
 columnphysics/icepack_atmo.F90       | 14 +++++++++++---
 columnphysics/icepack_parameters.F90 | 22 ++++++++++++++--------
 configuration/driver/icedrv_init.F90 | 10 ++++++----
 3 files changed, 31 insertions(+), 15 deletions(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index 94d30c518..176ad5faa 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -18,7 +18,7 @@ module icepack_atmo
       use icepack_parameters,  only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
       use icepack_parameters,  only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
       use icepack_parameters,  only: pih, dragio, rhoi, rhos, rhow
-      use icepack_parameters, only: atmbndy, calc_strair, formdrag
+      use icepack_parameters, only: atmbndy, heatflux_linear, calc_strair, formdrag
       use icepack_tracers, only: n_iso
       use icepack_tracers, only: tr_iso
       use icepack_warnings, only: warnstr, icepack_warnings_add
@@ -359,8 +359,16 @@ subroutine atmo_boundary_layer (sfctype,            &
       ! as in Jordan et al (JGR, 1999)
       !------------------------------------------------------------
 
-      shcoef = rhoa * ustar * cp * rh + c1
-      lhcoef = rhoa * ustar * Lheat  * re
+      if (heatflux_linear) then
+         !- Use constant coefficients for sensible and latent heat fluxes
+         !    similar to atmo_boundary_const but using (wind-Vice) instead of wind-only
+         shcoef = (1.20e-3_dbl_kind)*cp_air*rhoa*vmag
+         lhcoef = (1.50e-3_dbl_kind)*Lheat *rhoa*vmag
+      else
+         !- Use Jordan et al (JGR, 1999) formulation
+         shcoef = rhoa * ustar * cp * rh + c1
+         lhcoef = rhoa * ustar * Lheat  * re
+      endif
 
       !------------------------------------------------------------
       ! Compute diagnostics: 2m ref T, Q, U
diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index 6191de996..9996716a5 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -250,11 +250,12 @@ module icepack_parameters
          TTTocn = 5107.4_dbl_kind      ! for qsat over ocn
 
       character (len=char_len), public :: &
-         atmbndy = 'default' ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy = 'default'           ! atmo boundary method, 'default' ('ccsm3') or 'constant'
 
       logical (kind=log_kind), public :: &
          calc_strair     = .true.  , & ! if true, calculate wind stress
          formdrag        = .false. , & ! if true, calculate form drag
+         heatflux_linear = .false. , & ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
          highfreq        = .false.     ! if true, calculate high frequency coupling
 
       integer (kind=int_kind), public :: &
@@ -446,7 +447,7 @@ subroutine icepack_init_parameters(   &
          albicev_in, albicei_in, albsnowv_in, &
          ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, &
          kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, &
-         atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
+         atmbndy_in, heatflux_linear_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
          atmiter_conv_in, calc_dragio_in, &
          tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
          floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
@@ -646,12 +647,13 @@ subroutine icepack_init_parameters(   &
          TTTocn_in        ! for qsat over ocn
 
       character (len=*), intent(in), optional :: &
-         atmbndy_in ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy_in   ! atmo boundary method, 'default' ('ccsm3') or 'constant'
         
       logical (kind=log_kind), intent(in), optional :: &
-         calc_strair_in, &  ! if true, calculate wind stress components
-         formdrag_in,    &  ! if true, calculate form drag
-         highfreq_in        ! if true, use high frequency coupling
+         calc_strair_in,     & ! if true, calculate wind stress components
+         formdrag_in,        & ! if true, calculate form drag
+         heatflux_linear_in, & ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
+         highfreq_in           ! if true, use high frequency coupling
         
       integer (kind=int_kind), intent(in), optional :: &
          natmiter_in        ! number of iterations for boundary layer calculations
@@ -902,6 +904,7 @@ subroutine icepack_init_parameters(   &
       if (present(krdg_redist_in)       ) krdg_redist      = krdg_redist_in
       if (present(mu_rdg_in)            ) mu_rdg           = mu_rdg_in
       if (present(atmbndy_in)           ) atmbndy          = atmbndy_in
+      if (present(heatflux_linear_in)   ) heatflux_linear  = heatflux_linear_in
       if (present(calc_strair_in)       ) calc_strair      = calc_strair_in
       if (present(formdrag_in)          ) formdrag         = formdrag_in
       if (present(highfreq_in)          ) highfreq         = highfreq_in
@@ -1117,7 +1120,7 @@ subroutine icepack_query_parameters(   &
          albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, &
          rsnw_mlt_out, dEdd_algae_out, &
          kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, &
-         atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
+         atmbndy_out, heatflux_linear_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
          atmiter_conv_out, calc_dragio_out, &
          tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
          floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
@@ -1326,7 +1329,8 @@ subroutine icepack_query_parameters(   &
          TTTocn_out        ! for qsat over ocn
 
       character (len=*), intent(out), optional :: &
-         atmbndy_out ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy_out,        & ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         heatflux_linear_out   ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
         
       logical (kind=log_kind), intent(out), optional :: &
          calc_strair_out, &  ! if true, calculate wind stress components
@@ -1622,6 +1626,7 @@ subroutine icepack_query_parameters(   &
       if (present(krdg_redist_out)       ) krdg_redist_out  = krdg_redist
       if (present(mu_rdg_out)            ) mu_rdg_out       = mu_rdg
       if (present(atmbndy_out)           ) atmbndy_out      = atmbndy
+      if (present(heatflux_linear_out)   ) heatflux_linear  = heatflux_linear
       if (present(calc_strair_out)       ) calc_strair_out  = calc_strair
       if (present(formdrag_out)          ) formdrag_out     = formdrag
       if (present(highfreq_out)          ) highfreq_out     = highfreq
@@ -1826,6 +1831,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  krdg_redist   = ", krdg_redist
         write(iounit,*) "  mu_rdg        = ", mu_rdg
         write(iounit,*) "  atmbndy       = ", atmbndy
+        write(iounit,*) "  heatflux_linear = ", heatflux_linear
         write(iounit,*) "  calc_strair   = ", calc_strair
         write(iounit,*) "  formdrag      = ", formdrag
         write(iounit,*) "  highfreq      = ", highfreq
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index c4a2c7a17..321bf0f9f 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -105,7 +105,8 @@ subroutine input_data
       ! Flux convergence tolerance
       real (kind=dbl_kind) :: atmiter_conv
 
-      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio
+      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio, &
+                                 heatflux_linear
       logical (kind=log_kind) :: conserv_check
 
       integer (kind=int_kind) :: ntrcr
@@ -166,7 +167,7 @@ subroutine input_data
         update_ocn_f,    l_mpond_fresh,   ustar_min,       &
         fbot_xfer_type,  oceanmixed_ice,  emissivity,      &
         formdrag,        highfreq,        natmiter,        &
-        atmiter_conv,    calc_dragio,                      &
+        atmiter_conv,    calc_dragio,     heatflux_linear, &
         tfrz_option,     default_season,  wave_spec_type,  &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
@@ -200,7 +201,7 @@ subroutine input_data
            R_snw_out=R_snw, dT_mlt_out=dT_mlt, rsnw_mlt_out=rsnw_mlt, &
            kstrength_out=kstrength, krdg_partic_out=krdg_partic, &
            krdg_redist_out=krdg_redist, mu_rdg_out=mu_rdg, &
-           atmbndy_out=atmbndy, calc_strair_out=calc_strair, &
+           atmbndy_out=atmbndy, heatflux_linear_out=heatflux_linear, calc_strair_out=calc_strair, &
            formdrag_out=formdrag, highfreq_out=highfreq, &
            emissivity_out=emissivity, &
            kitd_out=kitd, kcatbound_out=kcatbound, hs0_out=hs0, & 
@@ -704,6 +705,7 @@ subroutine input_data
 
          write(nu_diag,1030) ' atmbndy                   = ', &
                                trim(atmbndy)
+         write(nu_diag,1010) ' heatflux_linear           = ', heatflux_linear
          write(nu_diag,1010) ' formdrag                  = ', formdrag
          write(nu_diag,1010) ' highfreq                  = ', highfreq
          write(nu_diag,1020) ' natmiter                  = ', natmiter
@@ -911,7 +913,7 @@ subroutine input_data
            R_snw_in=R_snw, dT_mlt_in=dT_mlt, rsnw_mlt_in=rsnw_mlt, &
            kstrength_in=kstrength, krdg_partic_in=krdg_partic, &
            krdg_redist_in=krdg_redist, mu_rdg_in=mu_rdg, &
-           atmbndy_in=atmbndy, calc_strair_in=calc_strair, &
+           atmbndy_in=atmbndy, heatflux_linear_in=heatflux_linear, calc_strair_in=calc_strair, &
            formdrag_in=formdrag, highfreq_in=highfreq, &
            emissivity_in=emissivity, &
            kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, &

From aa871efc8c2a35664d8ba3f4dbb3d51760b16d73 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Tue, 14 Sep 2021 20:37:32 +0000
Subject: [PATCH 02/20] ERROR correction: heatflux_linear is logical

---
 columnphysics/icepack_parameters.F90 | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index 9996716a5..86b043f9c 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -1329,13 +1329,13 @@ subroutine icepack_query_parameters(   &
          TTTocn_out        ! for qsat over ocn
 
       character (len=*), intent(out), optional :: &
-         atmbndy_out,        & ! atmo boundary method, 'default' ('ccsm3') or 'constant'
-         heatflux_linear_out   ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
+         atmbndy_out     & ! atmo boundary method, 'default' ('ccsm3') or 'constant'
         
       logical (kind=log_kind), intent(out), optional :: &
-         calc_strair_out, &  ! if true, calculate wind stress components
-         formdrag_out,    &  ! if true, calculate form drag
-         highfreq_out        ! if true, use high frequency coupling
+         calc_strair_out,     & ! if true, calculate wind stress components
+         formdrag_out,        & ! if true, calculate form drag
+         heatflux_linear_out, & ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
+         highfreq_out           ! if true, use high frequency coupling
         
       integer (kind=int_kind), intent(out), optional :: &
          natmiter_out        ! number of iterations for boundary layer calculations

From d602c2c2e70bbb365c3e0c60a21b5e9fa68f7284 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Tue, 14 Sep 2021 20:51:46 +0000
Subject: [PATCH 03/20] ERROR/syntax correction.

---
 columnphysics/icepack_parameters.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index 86b043f9c..4ddf226e9 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -1329,7 +1329,7 @@ subroutine icepack_query_parameters(   &
          TTTocn_out        ! for qsat over ocn
 
       character (len=*), intent(out), optional :: &
-         atmbndy_out     & ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy_out   ! atmo boundary method, 'default' ('ccsm3') or 'constant'
         
       logical (kind=log_kind), intent(out), optional :: &
          calc_strair_out,     & ! if true, calculate wind stress components

From 2da0b98935e17cfaa51c66dd884ad00a26195f74 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Wed, 15 Sep 2021 13:13:53 +0000
Subject: [PATCH 04/20] Add option atmbndy='mixed' boundary layer condition

---
 columnphysics/icepack_atmo.F90       | 12 ++++++------
 columnphysics/icepack_parameters.F90 | 12 +++---------
 configuration/driver/icedrv_init.F90 | 14 ++++++--------
 3 files changed, 15 insertions(+), 23 deletions(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index 176ad5faa..20559f0d9 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -18,7 +18,7 @@ module icepack_atmo
       use icepack_parameters,  only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
       use icepack_parameters,  only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
       use icepack_parameters,  only: pih, dragio, rhoi, rhos, rhow
-      use icepack_parameters, only: atmbndy, heatflux_linear, calc_strair, formdrag
+      use icepack_parameters, only: atmbndy, calc_strair, formdrag
       use icepack_tracers, only: n_iso
       use icepack_tracers, only: tr_iso
       use icepack_warnings, only: warnstr, icepack_warnings_add
@@ -359,13 +359,13 @@ subroutine atmo_boundary_layer (sfctype,            &
       ! as in Jordan et al (JGR, 1999)
       !------------------------------------------------------------
 
-      if (heatflux_linear) then
+      if (trim(atmbndy) == 'mixed') then
          !- Use constant coefficients for sensible and latent heat fluxes
-         !    similar to atmo_boundary_const but using (wind-Vice) instead of wind-only
+         !    similar to atmo_boundary_const but using vmag instead of wind
          shcoef = (1.20e-3_dbl_kind)*cp_air*rhoa*vmag
          lhcoef = (1.50e-3_dbl_kind)*Lheat *rhoa*vmag
-      else
-         !- Use Jordan et al (JGR, 1999) formulation
+      else ! 'similarity'
+         !- Monin-Obukhov similarity theory for boundary layer
          shcoef = rhoa * ustar * cp * rh + c1
          lhcoef = rhoa * ustar * Lheat  * re
       endif
@@ -960,7 +960,7 @@ subroutine icepack_atm_boundary(sfctype,                   &
                                    delt,     delq,     &
                                    lhcoef,   shcoef    )
          if (icepack_warnings_aborted(subname)) return
-      else ! default
+      else ! 'similarity' or 'mixed'
          call atmo_boundary_layer (sfctype,                 &
                                    calc_strair, formdrag,   &
                                    Tsf,      potT,          &
diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index 4ddf226e9..d3de23dd6 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -250,12 +250,11 @@ module icepack_parameters
          TTTocn = 5107.4_dbl_kind      ! for qsat over ocn
 
       character (len=char_len), public :: &
-         atmbndy = 'default'           ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy = 'similarity'        ! atmo boundary method, 'similarity' ('ccsm3'), 'constant' or 'mixed'
 
       logical (kind=log_kind), public :: &
          calc_strair     = .true.  , & ! if true, calculate wind stress
          formdrag        = .false. , & ! if true, calculate form drag
-         heatflux_linear = .false. , & ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
          highfreq        = .false.     ! if true, calculate high frequency coupling
 
       integer (kind=int_kind), public :: &
@@ -447,7 +446,7 @@ subroutine icepack_init_parameters(   &
          albicev_in, albicei_in, albsnowv_in, &
          ahmax_in, R_ice_in, R_pnd_in, R_snw_in, dT_mlt_in, rsnw_mlt_in, &
          kalg_in, kstrength_in, krdg_partic_in, krdg_redist_in, mu_rdg_in, &
-         atmbndy_in, heatflux_linear_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
+         atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
          atmiter_conv_in, calc_dragio_in, &
          tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
          floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
@@ -652,7 +651,6 @@ subroutine icepack_init_parameters(   &
       logical (kind=log_kind), intent(in), optional :: &
          calc_strair_in,     & ! if true, calculate wind stress components
          formdrag_in,        & ! if true, calculate form drag
-         heatflux_linear_in, & ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
          highfreq_in           ! if true, use high frequency coupling
         
       integer (kind=int_kind), intent(in), optional :: &
@@ -904,7 +902,6 @@ subroutine icepack_init_parameters(   &
       if (present(krdg_redist_in)       ) krdg_redist      = krdg_redist_in
       if (present(mu_rdg_in)            ) mu_rdg           = mu_rdg_in
       if (present(atmbndy_in)           ) atmbndy          = atmbndy_in
-      if (present(heatflux_linear_in)   ) heatflux_linear  = heatflux_linear_in
       if (present(calc_strair_in)       ) calc_strair      = calc_strair_in
       if (present(formdrag_in)          ) formdrag         = formdrag_in
       if (present(highfreq_in)          ) highfreq         = highfreq_in
@@ -1120,7 +1117,7 @@ subroutine icepack_query_parameters(   &
          albsnowi_out, ahmax_out, R_ice_out, R_pnd_out, R_snw_out, dT_mlt_out, &
          rsnw_mlt_out, dEdd_algae_out, &
          kalg_out, kstrength_out, krdg_partic_out, krdg_redist_out, mu_rdg_out, &
-         atmbndy_out, heatflux_linear_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
+         atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
          atmiter_conv_out, calc_dragio_out, &
          tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
          floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
@@ -1334,7 +1331,6 @@ subroutine icepack_query_parameters(   &
       logical (kind=log_kind), intent(out), optional :: &
          calc_strair_out,     & ! if true, calculate wind stress components
          formdrag_out,        & ! if true, calculate form drag
-         heatflux_linear_out, & ! if true, calculate sensible+latent heatfluxes using traditional linear bulk formula
          highfreq_out           ! if true, use high frequency coupling
         
       integer (kind=int_kind), intent(out), optional :: &
@@ -1626,7 +1622,6 @@ subroutine icepack_query_parameters(   &
       if (present(krdg_redist_out)       ) krdg_redist_out  = krdg_redist
       if (present(mu_rdg_out)            ) mu_rdg_out       = mu_rdg
       if (present(atmbndy_out)           ) atmbndy_out      = atmbndy
-      if (present(heatflux_linear_out)   ) heatflux_linear  = heatflux_linear
       if (present(calc_strair_out)       ) calc_strair_out  = calc_strair
       if (present(formdrag_out)          ) formdrag_out     = formdrag
       if (present(highfreq_out)          ) highfreq_out     = highfreq
@@ -1831,7 +1826,6 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  krdg_redist   = ", krdg_redist
         write(iounit,*) "  mu_rdg        = ", mu_rdg
         write(iounit,*) "  atmbndy       = ", atmbndy
-        write(iounit,*) "  heatflux_linear = ", heatflux_linear
         write(iounit,*) "  calc_strair   = ", calc_strair
         write(iounit,*) "  formdrag      = ", formdrag
         write(iounit,*) "  highfreq      = ", highfreq
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 321bf0f9f..2a57d7d68 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -105,8 +105,7 @@ subroutine input_data
       ! Flux convergence tolerance
       real (kind=dbl_kind) :: atmiter_conv
 
-      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio, &
-                                 heatflux_linear
+      logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio
       logical (kind=log_kind) :: conserv_check
 
       integer (kind=int_kind) :: ntrcr
@@ -167,7 +166,7 @@ subroutine input_data
         update_ocn_f,    l_mpond_fresh,   ustar_min,       &
         fbot_xfer_type,  oceanmixed_ice,  emissivity,      &
         formdrag,        highfreq,        natmiter,        &
-        atmiter_conv,    calc_dragio,     heatflux_linear, &
+        atmiter_conv,    calc_dragio,                      &
         tfrz_option,     default_season,  wave_spec_type,  &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
@@ -201,7 +200,7 @@ subroutine input_data
            R_snw_out=R_snw, dT_mlt_out=dT_mlt, rsnw_mlt_out=rsnw_mlt, &
            kstrength_out=kstrength, krdg_partic_out=krdg_partic, &
            krdg_redist_out=krdg_redist, mu_rdg_out=mu_rdg, &
-           atmbndy_out=atmbndy, heatflux_linear_out=heatflux_linear, calc_strair_out=calc_strair, &
+           atmbndy_out=atmbndy, calc_strair_out=calc_strair, &
            formdrag_out=formdrag, highfreq_out=highfreq, &
            emissivity_out=emissivity, &
            kitd_out=kitd, kcatbound_out=kcatbound, hs0_out=hs0, & 
@@ -564,8 +563,8 @@ subroutine input_data
       if (formdrag) then
       if (trim(atmbndy) == 'constant') then
          write (nu_diag,*) 'WARNING: atmbndy = constant not allowed with formdrag'
-         write (nu_diag,*) 'WARNING: Setting atmbndy = default'
-         atmbndy = 'default'
+         write (nu_diag,*) 'WARNING: Setting atmbndy = similarity'
+         atmbndy = 'similarity'
       endif
 
       if (.not. calc_strair) then
@@ -705,7 +704,6 @@ subroutine input_data
 
          write(nu_diag,1030) ' atmbndy                   = ', &
                                trim(atmbndy)
-         write(nu_diag,1010) ' heatflux_linear           = ', heatflux_linear
          write(nu_diag,1010) ' formdrag                  = ', formdrag
          write(nu_diag,1010) ' highfreq                  = ', highfreq
          write(nu_diag,1020) ' natmiter                  = ', natmiter
@@ -913,7 +911,7 @@ subroutine input_data
            R_snw_in=R_snw, dT_mlt_in=dT_mlt, rsnw_mlt_in=rsnw_mlt, &
            kstrength_in=kstrength, krdg_partic_in=krdg_partic, &
            krdg_redist_in=krdg_redist, mu_rdg_in=mu_rdg, &
-           atmbndy_in=atmbndy, heatflux_linear_in=heatflux_linear, calc_strair_in=calc_strair, &
+           atmbndy_in=atmbndy, calc_strair_in=calc_strair, &
            formdrag_in=formdrag, highfreq_in=highfreq, &
            emissivity_in=emissivity, &
            kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, &

From dd63009531c8a0443c6efd7e173ae807a70788ba Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Wed, 15 Sep 2021 21:24:11 +0000
Subject: [PATCH 05/20] For 'atmbndy': change 'default' (obsolete) to
 'similarity'. Same physics

---
 configuration/scripts/icepack_in | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in
index ff5f4b166..dad73618e 100644
--- a/configuration/scripts/icepack_in
+++ b/configuration/scripts/icepack_in
@@ -97,7 +97,7 @@
 
 &forcing_nml
     formdrag        = .false.
-    atmbndy         = 'default'
+    atmbndy         = 'similarity'
     calc_strair     = .true.
     calc_Tsfc       = .true.
     highfreq        = .false.

From bd91a7d159a3ea6b5818073f0ce16eec05b77e10 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Wed, 15 Sep 2021 21:34:36 +0000
Subject: [PATCH 06/20] For 'atmbndy': change 'default' (obsolete) to
 'similarity'. Same physics

---
 configuration/scripts/options/set_nml.dyn | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configuration/scripts/options/set_nml.dyn b/configuration/scripts/options/set_nml.dyn
index f95f0c60e..623438a10 100644
--- a/configuration/scripts/options/set_nml.dyn
+++ b/configuration/scripts/options/set_nml.dyn
@@ -3,5 +3,5 @@ kstrength = 0
 krdg_partic = 0
 krdg_redist = 0
 formdrag = .false.
-atmbndy = 'constant'
+atmbndy = 'similarity'
 highfreq = .false.

From b785a57ec39dfedea4ffed8f745bcd56d994fe39 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Wed, 15 Sep 2021 21:35:13 +0000
Subject: [PATCH 07/20] New options for 'atmbndy'

---
 doc/source/icepack_index.rst               | 2 +-
 doc/source/user_guide/interfaces.include   | 4 ++--
 doc/source/user_guide/ug_case_settings.rst | 6 ++++--
 3 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst
index 293ddd797..544d8fe54 100755
--- a/doc/source/icepack_index.rst
+++ b/doc/source/icepack_index.rst
@@ -52,7 +52,7 @@ either Celsius or Kelvin units).
    "atmiter_conv", ":math:`\bullet` convergence criteria for ustar", "0.0"
    "atm_data_format", ":math:`\bullet` format of atmospheric forcing files", ""
    "atm_data_type", ":math:`\bullet` type of atmospheric forcing", ""
-   "atmbndy", ":math:`\bullet` atmo boundary layer parameterization (‘default’ or ‘constant’)", ""
+   "atmbndy", ":math:`\bullet` atmo boundary layer parameterization (‘similarity’,‘constant’ or 'mixed')", ""
    "awtidf", "weighting factor for near-ir, diffuse albedo", "0.36218"
    "awtidr", "weighting factor for near-ir, direct albedo", "0.00182"
    "awtvdf", "weighting factor for visible, diffuse albedo", "0.63282"
diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include
index 907959bf4..6cbe1ce34 100644
--- a/doc/source/user_guide/interfaces.include
+++ b/doc/source/user_guide/interfaces.include
@@ -944,7 +944,7 @@ icepack_init_parameters
            TTTocn_in        ! for qsat over ocn
   
         character (len=char_len), intent(in), optional :: &
-           atmbndy_in ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+           atmbndy_in ! atmo boundary method, 'similarity' ('ccsm3') 'constant' or 'mixed'
           
         logical (kind=log_kind), intent(in), optional :: &
            calc_strair_in, &  ! if true, calculate wind stress components
@@ -1306,7 +1306,7 @@ icepack_query_parameters
            TTTocn_out        ! for qsat over ocn
   
         character (len=char_len), intent(out), optional :: &
-           atmbndy_out ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+           atmbndy_out ! atmo boundary method, 'similarity' ('ccsm3'), 'constant' or 'mixed'
           
         logical (kind=log_kind), intent(out), optional :: &
            calc_strair_out, &  ! if true, calculate wind stress components
diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst
index 7ee316339..768f64f1f 100755
--- a/doc/source/user_guide/ug_case_settings.rst
+++ b/doc/source/user_guide/ug_case_settings.rst
@@ -314,8 +314,10 @@ forcing_nml
    :widths: 15, 15, 30, 15 
 
    "", "", "", ""
-   "``atmbndy``", "``constant``", "bulk transfer coefficients", "``default``"
-   "", "``default``", "stability-based boundary layer", ""
+   "``atmbndy``", "string", "bulk transfer coefficients", "``similarity``"
+   "", "``similarity``", "stability-based boundary layer", ""
+   "", "``constant``", "constant-based boundary layer", ""
+   "", "``mixed``", "stability-based, but constant for sensible+latent heatfluxes", ""
    "``atmiter_conv``", "real", "convergence criteria for ustar", "0.0"
    "``atm_data_file``", "string", "file containing atmospheric data", "' '"
    "``atm_data_format``", "``bin``", "read direct access binary forcing files", "``bin``"

From ffb374d8b5dc759a9a08b8dec24b433761614578 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 16 Sep 2021 12:55:32 +0000
Subject: [PATCH 08/20] "atmbndy" options: 'similarity', 'constant', 'mixed'

---
 doc/source/science_guide/sg_boundary_forcing.rst | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst
index e3006ed5f..5d8a454fa 100755
--- a/doc/source/science_guide/sg_boundary_forcing.rst
+++ b/doc/source/science_guide/sg_boundary_forcing.rst
@@ -120,7 +120,12 @@ parameterizations. Rain and all melted snow end up in the ocean.
 Wind stress and transfer coefficients for the
 turbulent heat fluxes are computed in subroutine
 *atmo\_boundary\_layer* following :cite:`Kauffman02`, with additions and changes as detailed in Appendix A of :cite:`Roberts15` for high frequency coupling (namelist variable ``highfreq``).
-The resulting equations are provided here.
+The resulting equations are provided here for the default boundary layer 
+scheme, which is based on Monin-Obukhov theory (``atmbndy`` = ‘stability’). 
+Alternatively, ``atmbndy`` = ‘constant’ provides constant coefficients for
+wind stress, sensible heat and latent heat calculations (computed in subroutine
+*atmo\_boundary\_const*); ``atmbndy`` = ‘mixed’ uses the stability based 
+calculation for wind stress and constant coefficients for sensible and latent heat fluxes.
 
 The wind stress and turbulent heat flux calculation accounts for both
 stable and unstable atmosphere–ice boundary layers. We first define the

From 4cfa3abdf7e0a35ddd49d138519e2611c93ae453 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 16 Sep 2021 12:58:10 +0000
Subject: [PATCH 09/20] Abort if "atmbndy" option is unknown + set
 default=similarity for backward compability

---
 columnphysics/icepack_atmo.F90 | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index 20559f0d9..49e15f8c5 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -950,6 +950,10 @@ subroutine icepack_atm_boundary(sfctype,                   &
 
       Cdn_atm_ratio_n = c1
 
+      if (trim(atmbndy) == 'default') then
+         atmbndy = 'similarity' ! For backward compability only ...
+      endif
+
       if (trim(atmbndy) == 'constant') then
          call atmo_boundary_const (sfctype,  calc_strair, &
                                    uatm,     vatm,     &
@@ -960,7 +964,8 @@ subroutine icepack_atm_boundary(sfctype,                   &
                                    delt,     delq,     &
                                    lhcoef,   shcoef    )
          if (icepack_warnings_aborted(subname)) return
-      else ! 'similarity' or 'mixed'
+      elseif (trim(atmbndy) == 'similarity' .or. &
+              trim(atmbndy) == 'mixed') then
          call atmo_boundary_layer (sfctype,                 &
                                    calc_strair, formdrag,   &
                                    Tsf,      potT,          &
@@ -979,6 +984,10 @@ subroutine icepack_atm_boundary(sfctype,                   &
                                    uvel=l_uvel, vvel=l_vvel,  &
                                    Uref=l_Uref, zlvs=zlvs)
          if (icepack_warnings_aborted(subname)) return
+      else
+         call icepack_warnings_add( &
+            subname//' atmbndy = '//trim(atmbndy)//' : unknown value')
+         call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
       endif ! atmbndy
 
       if (present(Uref)) then

From 644d32407219eb137e0879cc245069fb93b5d8e5 Mon Sep 17 00:00:00 2001
From: Mads Hvid Ribergaard <38077893+mhrib@users.noreply.github.com>
Date: Thu, 16 Sep 2021 16:22:30 +0200
Subject: [PATCH 10/20] Update doc/source/science_guide/sg_boundary_forcing.rst

Co-authored-by: Philippe Blain <levraiphilippeblain@gmail.com>
---
 doc/source/science_guide/sg_boundary_forcing.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst
index 5d8a454fa..e006c162b 100755
--- a/doc/source/science_guide/sg_boundary_forcing.rst
+++ b/doc/source/science_guide/sg_boundary_forcing.rst
@@ -121,10 +121,10 @@ Wind stress and transfer coefficients for the
 turbulent heat fluxes are computed in subroutine
 *atmo\_boundary\_layer* following :cite:`Kauffman02`, with additions and changes as detailed in Appendix A of :cite:`Roberts15` for high frequency coupling (namelist variable ``highfreq``).
 The resulting equations are provided here for the default boundary layer 
-scheme, which is based on Monin-Obukhov theory (``atmbndy`` = ‘stability’). 
-Alternatively, ``atmbndy`` = ‘constant’ provides constant coefficients for
+scheme, which is based on Monin-Obukhov theory (``atmbndy = 'stability'``). 
+Alternatively, ``atmbndy = 'constant'`` provides constant coefficients for
 wind stress, sensible heat and latent heat calculations (computed in subroutine
-*atmo\_boundary\_const*); ``atmbndy`` = ‘mixed’ uses the stability based 
+*atmo\_boundary\_const*); ``atmbndy = 'mixed'`` uses the stability based 
 calculation for wind stress and constant coefficients for sensible and latent heat fluxes.
 
 The wind stress and turbulent heat flux calculation accounts for both

From ef391d5ed236a2ee9a831edd2e168db13ddf50b2 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 16 Sep 2021 14:31:03 +0000
Subject: [PATCH 11/20] atmbndy option 'ccsm' not an option. Previous given as
 'default'

---
 columnphysics/icepack_parameters.F90 | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index d3de23dd6..f3f81c253 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -250,7 +250,7 @@ module icepack_parameters
          TTTocn = 5107.4_dbl_kind      ! for qsat over ocn
 
       character (len=char_len), public :: &
-         atmbndy = 'similarity'        ! atmo boundary method, 'similarity' ('ccsm3'), 'constant' or 'mixed'
+         atmbndy = 'similarity'        ! atmo boundary method, 'similarity', 'constant' or 'mixed'
 
       logical (kind=log_kind), public :: &
          calc_strair     = .true.  , & ! if true, calculate wind stress
@@ -646,7 +646,7 @@ subroutine icepack_init_parameters(   &
          TTTocn_in        ! for qsat over ocn
 
       character (len=*), intent(in), optional :: &
-         atmbndy_in   ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy_in   ! atmo boundary method, 'similarity', 'constant' or 'mixed'
         
       logical (kind=log_kind), intent(in), optional :: &
          calc_strair_in,     & ! if true, calculate wind stress components
@@ -1326,7 +1326,7 @@ subroutine icepack_query_parameters(   &
          TTTocn_out        ! for qsat over ocn
 
       character (len=*), intent(out), optional :: &
-         atmbndy_out   ! atmo boundary method, 'default' ('ccsm3') or 'constant'
+         atmbndy_out   ! atmo boundary method, 'similarity', 'constant' or 'mixed'
         
       logical (kind=log_kind), intent(out), optional :: &
          calc_strair_out,     & ! if true, calculate wind stress components

From 2c10349d01df06e6505bd97efd63426e34730f63 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 16 Sep 2021 14:32:12 +0000
Subject: [PATCH 12/20] namelist parameter plus string inside double quotes

---
 doc/source/science_guide/sg_boundary_forcing.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst
index 5d8a454fa..5d99962d3 100755
--- a/doc/source/science_guide/sg_boundary_forcing.rst
+++ b/doc/source/science_guide/sg_boundary_forcing.rst
@@ -121,10 +121,10 @@ Wind stress and transfer coefficients for the
 turbulent heat fluxes are computed in subroutine
 *atmo\_boundary\_layer* following :cite:`Kauffman02`, with additions and changes as detailed in Appendix A of :cite:`Roberts15` for high frequency coupling (namelist variable ``highfreq``).
 The resulting equations are provided here for the default boundary layer 
-scheme, which is based on Monin-Obukhov theory (``atmbndy`` = ‘stability’). 
-Alternatively, ``atmbndy`` = ‘constant’ provides constant coefficients for
+scheme, which is based on Monin-Obukhov theory (``atmbndy = ‘stability’``). 
+Alternatively, ``atmbndy = ‘constant’`` provides constant coefficients for
 wind stress, sensible heat and latent heat calculations (computed in subroutine
-*atmo\_boundary\_const*); ``atmbndy`` = ‘mixed’ uses the stability based 
+*atmo\_boundary\_const*); ``atmbndy = ‘mixed’`` uses the stability based 
 calculation for wind stress and constant coefficients for sensible and latent heat fluxes.
 
 The wind stress and turbulent heat flux calculation accounts for both

From a44c806086957574784107c4e308205a05ccc2c7 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 16 Sep 2021 14:46:34 +0000
Subject: [PATCH 13/20] Re-interduce file, after accidentally git rm command

---
 .../science_guide/sg_boundary_forcing.rst     | 461 ++++++++++++++++++
 1 file changed, 461 insertions(+)
 create mode 100755 doc/source/science_guide/sg_boundary_forcing.rst

diff --git a/doc/source/science_guide/sg_boundary_forcing.rst b/doc/source/science_guide/sg_boundary_forcing.rst
new file mode 100755
index 000000000..5d99962d3
--- /dev/null
+++ b/doc/source/science_guide/sg_boundary_forcing.rst
@@ -0,0 +1,461 @@
+:tocdepth: 3
+
+.. _boundary_forcing:
+
+Atmosphere and ocean boundary forcing
+=====================================
+
+.. _tab-flux-cpl:
+
+.. csv-table:: *External forcing data that are relevant to Icepack*
+   :header: "Variable", "Description", "External Interactions"
+   :widths: 10, 25, 25
+     
+   ":math:`z_o`", "Atmosphere level height", "From *atmosphere model* to *sea ice model*"
+   ":math:`z_{o,s}`", "Atmosphere level height (scalar quantities) (optional)", "From *atmosphere model* to *sea ice model*"
+   ":math:`\vec{U}_a`", "Wind velocity", "From *atmosphere model* to *sea ice model*"
+   ":math:`Q_a`", "Specific humidity", "From *atmosphere model* to *sea ice model*"
+   ":math:`\rho_a`", "Air density", "From *atmosphere model* to *sea ice model*"
+   ":math:`\Theta_a`", "Air potential temperature", "From *atmosphere model* to *sea ice model*"
+   ":math:`T_a`", "Air temperature", "From *atmosphere model* to *sea ice model*"
+   ":math:`F_{sw\downarrow}`", "Incoming shortwave radiation (4 bands)", "From *atmosphere model* to *sea ice model*"
+   ":math:`F_{L\downarrow}`", "Incoming longwave radiation", "From *atmosphere model* to *sea ice model*"
+   ":math:`F_{rain}`", "Rainfall rate", "From *atmosphere model* to *sea ice model*"
+   ":math:`F_{snow}`", "Snowfall rate", "From *atmosphere model* to *sea ice model*"
+   ":math:`F_{frzmlt}`", "Freezing/melting potential", "From *ocean model* to *sea ice model*"
+   ":math:`T_w`", "Sea surface temperature", "From *ocean model* to *sea ice model*"
+   ":math:`S`", "Sea surface salinity", "From *ocean model* to *sea ice model*"
+   ":math:`\nabla H_o`", "Sea surface slope", "From *ocean model* via flux coupler to *sea ice model*"
+   ":math:`h_1`", "Thickness of first ocean level (optional)", "From *ocean model* to *sea ice model*"
+   ":math:`\vec{U}_w`", "Surface ocean currents", "From *ocean model* to *sea ice model* (available in Icepack driver, not used directly in column physics)"
+   ":math:`\vec{\tau}_a`", "Wind stress", "From *sea ice model* to *atmosphere model*"
+   ":math:`F_s`", "Sensible heat flux", "From *sea ice model* to *atmosphere model*"
+   ":math:`F_l`", "Latent heat flux", "From *sea ice model* to *atmosphere model*"
+   ":math:`F_{L\uparrow}`", "Outgoing longwave radiation", "From *sea ice model* to *atmosphere model*"
+   ":math:`F_{evap}`", "Evaporated water", "From *sea ice model* to *atmosphere model*"
+   ":math:`\alpha`", "Surface albedo (4 bands)", "From *sea ice model* to *atmosphere model*"
+   ":math:`T_{sfc}`", "Surface temperature", "From *sea ice model* to *atmosphere model*"
+   ":math:`F_{sw\Downarrow}`", "Penetrating shortwave radiation", "From *sea ice model* to *ocean model*"
+   ":math:`F_{water}`", "Fresh water flux", "From *sea ice model* to *ocean model*"
+   ":math:`F_{hocn}`", "Net heat flux to ocean", "From *sea ice model* to *ocean model*"
+   ":math:`F_{salt}`", "Salt flux", "From *sea ice model* to *ocean model*"
+   ":math:`\vec{\tau}_w`", "Ice-ocean stress", "From *sea ice model* to *ocean model*"
+   ":math:`F_{bio}`", "Biogeochemical fluxes", "From *sea ice model* to *ocean model*"
+   ":math:`a_{i}`", "Ice fraction", "From *sea ice model* to both *ocean and atmosphere models*"
+   ":math:`T^{ref}_{a}`", "2m reference temperature (diagnostic)", "From *sea ice model* to both *ocean and atmosphere models*"
+   ":math:`Q^{ref}_{a}`", "2m reference humidity (diagnostic)", "From *sea ice model* to both *ocean and atmosphere models*"
+   ":math:`F_{swabs}`", "Absorbed shortwave (diagnostic)", "From *sea ice model* to both *ocean and atmosphere models*"
+   ":math:`E(f)`", "Wave spectrum as a function of frequency", "From *ocean surface wave model* to *sea ice model*"
+
+The ice fraction :math:`a_i` (aice) is the total fractional ice
+coverage of a grid cell. That is, in each cell,
+
+.. math::
+   \begin{array}{cl}
+                  a_{i}=0 & \mbox{if there is no ice} \\ 
+                  a_{i}=1 & \mbox{if there is no open water} \\ 
+                  0<a_{i}<1 & \mbox{if there is both ice and open water,}
+   \end{array}
+
+where :math:`a_{i}` is the sum of fractional ice areas for each category
+of ice. The ice fraction is used by the flux coupler to merge fluxes
+from the sea ice model with fluxes from the other earth system components. For example,
+the penetrating shortwave radiation flux, weighted by :math:`a_i`, is
+combined with the net shortwave radiation flux through ice-free leads,
+weighted by (:math:`1-a_i`), to obtain the net shortwave flux into the
+ocean over the entire grid cell. The CESM flux coupler requires the fluxes to
+be divided by the total ice area so that the ice and land models are
+treated identically (land also may occupy less than 100% of an
+atmospheric grid cell). These fluxes are "per unit ice area" rather than
+"per unit grid cell area."
+
+In some coupled climate models (for example, recent versions of the U.K.	
+Hadley Centre model) the surface air temperature and fluxes are computed	
+within the atmosphere model and are passed to CICE for use in the column physics. In this case the	
+logical parameter ``calc_Tsfc`` in *ice_therm_vertical* is set to false.	
+The fields ``fsurfn`` (the net surface heat flux from the atmosphere), ``flatn``	
+(the surface latent heat flux), and ``fcondtopn`` (the conductive flux at	
+the top surface) for each ice thickness category are copied or derived	
+from the input coupler fluxes and are passed to the thermodynamic driver	
+subroutine, *thermo_vertical*. At the end of the time step, the surface	
+temperature and effective conductivity (i.e., thermal conductivity	
+divided by thickness) of the top ice/snow layer in each category are	
+returned to the atmosphere model via the coupler. Since the ice surface	
+temperature is treated explicitly, the effective conductivity may need	
+to be limited to ensure stability. As a result, accuracy may be	
+significantly reduced, especially for thin ice or snow layers. A more	
+stable and accurate procedure would be to compute the temperature	
+profiles for both the atmosphere and ice, together with the surface	
+fluxes, in a single implicit calculation. This was judged impractical,	
+however, given that the atmosphere and sea ice models generally exist on	
+different grids and/or processor sets.
+
+.. _atmo:
+
+Atmosphere
+----------
+
+The wind velocity, specific humidity, air density and potential
+temperature at the given level height :math:`z_\circ` (optionally :math:`z_{\circ,s}`, see below) are used to
+compute transfer coefficients used in formulas for the surface wind
+stress and turbulent heat fluxes :math:`\vec\tau_a`, :math:`F_s`, and
+:math:`F_l`, as described below. The sensible and latent heat fluxes,
+:math:`F_s` and :math:`F_l`, along with shortwave and longwave
+radiation, :math:`F_{sw\downarrow}`, :math:`F_{L\downarrow}`
+and :math:`F_{L\uparrow}`, are included in the flux balance that
+determines the ice or snow surface temperature when ``calc_Tsfc`` is true.
+As described in the :ref:`thermo` section, these fluxes depend nonlinearly
+on the ice surface temperature :math:`T_{sfc}`. The balance
+equation is iterated until convergence, and the resulting fluxes and
+:math:`T_{sfc}` are then passed to the flux coupler.
+
+The snowfall precipitation rate (provided as liquid water equivalent and
+converted by the ice model to snow depth) also contributes to the heat
+and water mass budgets of the ice layer. Melt ponds generally form on
+the ice surface in the Arctic and refreeze later in the fall, reducing
+the total amount of fresh water that reaches the ocean and altering the
+heat budget of the ice; this version includes two new melt pond
+parameterizations. Rain and all melted snow end up in the ocean.
+
+Wind stress and transfer coefficients for the
+turbulent heat fluxes are computed in subroutine
+*atmo\_boundary\_layer* following :cite:`Kauffman02`, with additions and changes as detailed in Appendix A of :cite:`Roberts15` for high frequency coupling (namelist variable ``highfreq``).
+The resulting equations are provided here for the default boundary layer 
+scheme, which is based on Monin-Obukhov theory (``atmbndy = ‘stability’``). 
+Alternatively, ``atmbndy = ‘constant’`` provides constant coefficients for
+wind stress, sensible heat and latent heat calculations (computed in subroutine
+*atmo\_boundary\_const*); ``atmbndy = ‘mixed’`` uses the stability based 
+calculation for wind stress and constant coefficients for sensible and latent heat fluxes.
+
+The wind stress and turbulent heat flux calculation accounts for both
+stable and unstable atmosphere–ice boundary layers. We first define the
+"stability"
+
+.. math::
+   \Upsilon = {\kappa g z_\circ\over u^{*2}}
+   \left({\Theta^*\over\Theta_a\left(1+0.606Q_a\right)}  +
+   {Q^*\over 1/0.606 + Q_a}\right),
+   :label: upsilon
+
+where :math:`\kappa` is the von Karman constant, :math:`g` is
+gravitational acceleration, and :math:`u^*`, :math:`\Theta^*` and
+:math:`Q^*` are turbulent scales for velocity difference, temperature, and humidity,
+respectively, computed as (given the ice velocity :math:`\vec{U}_i`):
+
+.. math::
+   \begin{aligned}
+   u^*&=&c_u\;\textrm{max}\left(U_{\Delta\textrm{min}}, \left|\vec{U}_a - \vec{U}_i \right|\right), \\
+   \Theta^*&=& c_\theta\left(\Theta_a-T_{sfc}\right), \\
+   Q^*&=&c_q\left(Q_a-Q_{sfc}\right).
+   \end{aligned}
+   :label: stars
+
+
+Note that *atmo_boundary_layer* also accepts an optional argument, ``zlvs``, to support staggered atmospheric levels, i.e. receiving scalar quantities from the atmospheric model (humidity and temperature)
+at a different vertical level than the winds. In that case a separate stability :math:`\Upsilon_s` is computed using the same formula as above but substituting :math:`z_o` by :math:`z_{o,s}`.
+
+Within the :math:`u^*` expression, :math:`U_{\Delta\textrm{min}}` is the minimum allowable value of :math:`|\vec{U}_{a} - \vec{U}_{i}|` , which is set to of 0.5 m/s for high frequency coupling (``highfreq`` =.true.). 
+When high frequency coupling is turned off (``highfreq`` =.false.), it is assumed in equation :eq:`stars` that:
+
+.. math::
+ \vec{U}_{a} - \vec{U}_{i} \approx  \vec{U}_{a} 
+ :label: lowfreq 
+
+and a higher threshold is taken for :math:`U_{\Delta\textrm{min}}` of 1m/s. Equation :eq:`lowfreq` is a poor assumption when resolving inertial oscillations in ice-ocean configurations where the ice velocity vector may make a complete rotation over a period of :math:`\ge` 11.96 hours, as discussed in :cite:`Roberts15`.
+However, :eq:`lowfreq`  is acceptable for low frequency ice-ocean coupling on the order of a day or more, when transient ice-ocean Ekman transport is effectively filtered from the model solution.
+For the :math:`\Theta^*` and :math:`Q^*` terms in :eq:`stars`, :math:`T_{sfc}` and :math:`Q_{sfc}` are the surface temperature and specific
+humidity, respectively.  The latter is calculated by assuming a saturated
+surface, as described in the :ref:`sfc-forcing` section.
+
+Neglecting form drag, the exchange coefficients :math:`c_u`,
+:math:`c_\theta` and :math:`c_q` are initialized as
+
+.. math:: 
+   \kappa\over \ln(z_{ref}/z_{ice})
+   :label: kappa
+
+and updated during a short iteration, as they depend upon the turbulent
+scales. The number of iterations is set by the namelist variable
+``natmiter``, nominally set to five but sometimes increased by users employing the ``highfreq`` option.
+A convergence tolerance ``atmiter_conv`` on ``ustar`` can be set to exit the ``natmiter``
+loop early if desired.  Here, :math:`z_{ref}` is a reference height of 10m and
+:math:`z_{ice}` is the roughness length scale for the given
+sea ice category. :math:`\Upsilon` is constrained to have magnitude less
+than 10. Further, defining
+:math:`\chi = \left(1-16\Upsilon\right)^{0.25}` and :math:`\chi \geq 1`,
+the "integrated flux profiles" for momentum and stability in the
+unstable (:math:`\Upsilon <0`) case are given by
+
+.. math::
+   \begin{aligned}
+   \psi_m = &\mbox{}&2\ln\left[0.5(1+\chi)\right] +
+            \ln\left[0.5(1+\chi^2)\right] -2\tan^{-1}\chi +
+            {\pi\over 2}, \\
+   \psi_s = &\mbox{}&2\ln\left[0.5(1+\chi^2)\right].\end{aligned}
+   :label: psi1
+
+In a departure from the parameterization used in
+:cite:`Kauffman02`, we use profiles for the stable case
+following :cite:`Jordan99`,
+
+.. math::
+   \psi_m = \psi_s = -\left[0.7\Upsilon + 0.75\left(\Upsilon-14.3\right)
+            \exp\left(-0.35\Upsilon\right) + 10.7\right].
+   :label: psi2
+
+The coefficients are then updated as
+
+.. math::
+   \begin{aligned}
+   c_u^\prime&=&{c_u\over 1+c_u\left(\lambda-\psi_m\right)/\kappa} \\
+   c_\theta^\prime&=& {c_\theta\over 1+c_\theta\left(\lambda_s-\psi_s\right)/\kappa}\\
+   c_q^\prime&=&c_\theta^\prime\end{aligned}
+   :label: coeff1
+
+where :math:`\lambda = \ln\left(z_\circ/z_{ref}\right)` and :math:`\lambda_s = \ln\left(z_{\circ,s}/z_{ref}\right)` if staggered atmospheric levels are used, else :math:`\lambda_s=\lambda`. The
+first iteration ends with new turbulent scales from
+equations :eq:`stars`. After ``natmiter`` iterations the latent and sensible
+heat flux coefficients are computed, along with the wind stress:
+
+.. math::
+   \begin{aligned}
+   C_l&=&\rho_a \left(L_{vap}+L_{ice}\right) u^* c_q \\
+   C_s&=&\rho_a c_p u^* c_\theta^* + 1 \\
+   \vec{\tau}_a&=&{\rho_a (u^{*})^2 \left( \vec{U}_{a} - \vec{U}_{i} \right) \over  \left| \vec{U}_{a} - \vec{U}_{i} \right|}
+   \end{aligned}
+   :label: coeff2
+
+where :math:`L_{vap}` and :math:`L_{ice}` are
+latent heats of vaporization and fusion, :math:`\rho_a` is the density
+of air and :math:`c_p` is its specific heat. Again following
+:cite:`Jordan99`, we have added a constant to the sensible
+heat flux coefficient in order to allow some heat to pass between the
+atmosphere and the ice surface in stable, calm conditions. 
+For the atmospheric stress term in :eq:`coeff2`, we make the assumption in :eq:`lowfreq` when ``highfreq`` =.false..
+
+The atmospheric reference temperature :math:`T_a^{ref}` is computed from
+:math:`T_a` and :math:`T_{sfc}` using the coefficients
+:math:`c_u`, :math:`c_\theta` and :math:`c_q`. Although the sea ice
+model does not use this quantity, it is convenient for the ice model to
+perform this calculation. The atmospheric reference temperature is
+returned to the flux coupler as a climate diagnostic. The same is true
+for the reference humidity, :math:`Q_a^{ref}`.
+
+Additional details about the latent and sensible heat fluxes and other
+quantities referred to here can be found in
+the :ref:`sfc-forcing` section.
+
+.. _ocean:
+
+Ocean
+-----
+
+New sea ice forms when the ocean temperature drops below its freezing
+temperature. In the Bitz and Lipscomb thermodynamics,
+:cite:`Bitz99` :math:`T_f=-\mu S`, where :math:`S` is the
+seawater salinity and :math:`\mu=0.054^\circ`/ppt is the ratio of the
+freezing temperature of brine to its salinity (linear liquidus
+approximation). For the mushy thermodynamics, :math:`T_f` is given by a
+piecewise linear liquidus relation. The ocean model calculates the new
+ice formation; if the freezing/melting potential
+:math:`F_{frzmlt}` is positive, its value represents a certain
+amount of frazil ice that has formed in one or more layers of the ocean
+and floated to the surface. (The ocean model assumes that the amount of
+new ice implied by the freezing potential actually forms.)
+
+If :math:`F_{frzmlt}` is negative, it is used to heat already
+existing ice from below. In particular, the sea surface temperature and
+salinity are used to compute an oceanic heat flux :math:`F_w`
+(:math:`\left|F_w\right| \leq \left|F_{frzmlt}\right|`) which
+is applied at the bottom of the ice. The portion of the melting
+potential actually used to melt ice is returned to the coupler in
+:math:`F_{hocn}`. The ocean model adjusts its own heat budget
+with this quantity, assuming that the rest of the flux remained in the
+ocean.
+
+In addition to runoff from rain and melted snow, the fresh water flux
+:math:`F_{water}` includes ice melt water from the top surface
+and water frozen (a negative flux) or melted at the bottom surface of
+the ice. This flux is computed as the net change of fresh water in the
+ice and snow volume over the coupling time step, excluding frazil ice
+formation and newly accumulated snow. Setting the namelist option
+``update_ocn_f`` to true causes frazil ice to be included in the fresh
+water and salt fluxes.
+
+There is a flux of salt into the ocean under melting conditions, and a
+(negative) flux when sea water is freezing. However, melting sea ice
+ultimately freshens the top ocean layer, since the ocean is much more
+saline than the ice. The ice model passes the net flux of salt
+:math:`F_{salt}` to the flux coupler, based on the net change
+in salt for ice in all categories. In the present configuration,
+``ice_ref_salinity`` is used for computing the salt flux, although the ice
+salinity used in the thermodynamic calculation has differing values in
+the ice layers.
+
+A fraction of the incoming shortwave :math:`F_{sw\Downarrow}`
+penetrates the snow and ice layers and passes into the ocean, as
+described in the :ref:`sfc-forcing` section.
+
+A thermodynamic slab ocean mixed-layer parameterization is available 
+in **icepack\_ocean.F90** and can be run in the full CICE configuration.
+The turbulent fluxes are computed above the water surface using the same
+parameterizations as for sea ice, but with parameters appropriate for
+the ocean. The surface flux balance takes into account the turbulent
+fluxes, oceanic heat fluxes from below the mixed layer, and shortwave
+and longwave radiation, including that passing through the sea ice into
+the ocean. If the resulting sea surface temperature falls below the
+salinity-dependent freezing point, then new ice (frazil) forms.
+Otherwise, heat is made available for melting the ice.
+
+The ice-ocean drag coefficient, :math:`c_w`, can optionally be computed from the thickness of the first ocean level, :math:`h_1`, and an under-ice roughness length, :math:`z_{io}`.
+The computation follows :cite:`Roy15` :
+
+.. math::
+   c_w = c_w^* \lambda^2
+   :label: dragio
+
+where
+
+.. math::
+   \begin{aligned}
+   c_w^* &= \frac{\kappa^2} {\ln^2\left( h_1 /  z_{io} \right)}, \\
+   \lambda &= \frac{h_1 - z_{io}} {h_1 \left[ \sqrt{c_w^*} \kappa^{-1} \left( \ln(2) - 1 +  z_{io} / h_1 \right) + 1 \right] }
+   \end{aligned}
+   :label: dragio-defs
+
+.. _formdrag:
+
+Variable exchange coefficients
+------------------------------
+
+In the default configuration, atmospheric and oceanic neutral drag
+coefficients (:math:`c_u` and :math:`c_w`) are assumed constant in time
+and space. These constants are chosen to reflect friction associated
+with an effective sea ice surface roughness at the ice–atmosphere and
+ice–ocean interfaces. Sea ice (in both Arctic and Antarctic) contains
+pressure ridges as well as floe and melt pond edges that act as discrete
+obstructions to the flow of air or water past the ice, and are a source
+of form drag. Following :cite:`Tsamados14` and based on
+recent theoretical developments :cite:`Lupkes12,Lu11`, the
+neutral drag coefficients can now be estimated from properties of the
+ice cover such as ice concentration, vertical extent and area of the
+ridges, freeboard and floe draft, and size of floes and melt ponds. The
+new parameterization allows the drag coefficients to be coupled to the
+sea ice state and therefore to evolve spatially and temporally. This
+parameterization is contained in the subroutine *neutral\_drag\_coeffs*
+and is accessed by setting ``formdrag`` = true in the namelist.
+(Note:  see also :ref:`bugs`.)
+
+Following :cite:`Tsamados14`, consider the general case of
+fluid flow obstructed by N randomly oriented obstacles of height
+:math:`H` and transverse length :math:`L_y`, distributed on a domain
+surface area :math:`S_T`. Under the assumption of a logarithmic fluid
+velocity profile, the general formulation of the form drag coefficient
+can be expressed as
+
+.. math:: 
+   C_d=\frac{N c S_c^2 \gamma L_y  H}{2 S_T}\left[\frac{\ln(H/z_0)}{\ln(z_{ref}/z_0)}\right]^2,
+   :label: formdrag
+
+where :math:`z_0` is a roughness length parameter at the top or bottom
+surface of the ice, :math:`\gamma` is a geometric factor, :math:`c` is
+the resistance coefficient of a single obstacle, and :math:`S_c` is a
+sheltering function that takes into account the shielding effect of the
+obstacle,
+
+.. math:: 
+   S_{c}=\left(1-\exp(-s_l D/H)\right)^{1/2},
+   :label: shelter
+
+with :math:`D` the distance between two obstacles and :math:`s_l` an
+attenuation parameter.
+
+As in the original drag formulation in CICE (:ref:`atmo` and
+:ref:`ocean` sections), :math:`c_u` and :math:`c_w` along with the transfer
+coefficients for sensible heat, :math:`c_{\theta}`, and latent heat,
+:math:`c_{q}`, are initialized to a situation corresponding to neutral
+atmosphere–ice and ocean–ice boundary layers. The corresponding neutral
+exchange coefficients are then replaced by coefficients that explicitly
+account for form drag, expressed in terms of various contributions as
+
+.. math::
+   \tt{Cdn\_atm}  = \tt{Cdn\_atm\_rdg} + \tt{Cdn\_atm\_floe} + \tt{Cdn\_atm\_skin} + \tt{Cdn\_atm\_pond} ,
+   :label: Cda
+
+.. math::
+   \tt{Cdn\_ocn}  =  \tt{Cdn\_ocn\_rdg} + \tt{Cdn\_ocn\_floe} + \tt{Cdn\_ocn\_skin}. 
+   :label: Cdw
+
+The contributions to form drag from ridges (and keels underneath the
+ice), floe edges and melt pond edges can be expressed using the general
+formulation of equation :eq:`formdrag` (see :cite:`Tsamados14` for
+details). Individual terms in equation :eq:`Cdw` are fully described in
+:cite:`Tsamados14`. Following :cite:`Arya75`
+the skin drag coefficient is parametrized as
+
+.. math:: 
+   { \tt{Cdn\_(atm/ocn)\_skin}}=a_{i} \left(1-m_{(s/k)} \frac{H_{(s/k)}}{D_{(s/k)}}\right)c_{s(s/k)}, \mbox{       if  $\displaystyle\frac{H_{(s/k)}}{D_{(s/k)}}\ge\frac{1}{m_{(s/k)}}$,}
+   :label: skindrag
+
+where :math:`m_s` (:math:`m_k`) is a sheltering parameter that depends
+on the average sail (keel) height, :math:`H_s` (:math:`H_k`), but is
+often assumed constant, :math:`D_s` (:math:`D_k`) is the average
+distance between sails (keels), and :math:`c_{ss}` (:math:`c_{sk}`) is
+the unobstructed atmospheric (oceanic) skin drag that would be attained
+in the absence of sails (keels) and with complete ice coverage,
+:math:`a_{ice}=1`.
+
+Calculation of equations :eq:`formdrag` – :eq:`skindrag` requires that small-scale geometrical
+properties of the ice cover be related to average grid cell quantities
+already computed in the sea ice model. These intermediate quantities are
+briefly presented here and described in more detail in
+:cite:`Tsamados14`. The sail height is given by
+
+.. math:: 
+   H_{s} = \displaystyle 2\frac{v_{rdg}}{a_{rdg}}\left(\frac{\alpha\tan \alpha_{k} R_d+\beta \tan \alpha_{s} R_h}{\phi_r\tan \alpha_{k} R_d+\phi_k \tan \alpha_{s} R_h^2}\right),
+   :label: Hs
+
+and the distance between sails\ 
+
+.. math:: 
+   D_{s} = \displaystyle 2 H_s\frac{a_{i}}{a_{rdg}} \left(\frac{\alpha}{\tan \alpha_s}+\frac{\beta}{\tan \alpha_k}\frac{R_h}{R_d}\right),
+   :label: Ds
+
+where :math:`0<\alpha<1` and :math:`0<\beta<1` are weight functions,
+:math:`\alpha_{s}` and :math:`\alpha_{k}` are the sail and keel slope,
+:math:`\phi_s` and :math:`\phi_k` are constant porosities for the sails
+and keels, and we assume constant ratios for the average keel depth and
+sail height (:math:`H_k/H_s=R_h`) and for the average distances between
+keels and between sails (:math:`D_k/D_s=R_d`). With the assumption of
+hydrostatic equilibrium, the effective ice plus snow freeboard is
+:math:`H_{f}=\bar{h_i}(1-\rho_i/\rho_w)+\bar{h_s}(1-\rho_s/\rho_w)`,
+where :math:`\rho_i`, :math:`\rho_w` and :math:`\rho_s` are
+respectively the densities of sea ice, water and snow, :math:`\bar{h_i}`
+is the mean ice thickness and :math:`\bar{h_s}` is the mean snow
+thickness (means taken over the ice covered regions). For the melt pond
+edge elevation we assume that the melt pond surface is at the same level
+as the ocean surface surrounding the floes
+:cite:`Flocco07,Flocco10,Flocco12` and use the simplification
+:math:`H_p = H_f`. Finally to estimate the typical floe size
+:math:`L_A`, distance between floes, :math:`D_F`, and melt pond size,
+:math:`L_P` we use the parameterizations of :cite:`Lupkes12`
+to relate these quantities to the ice and pond concentrations. All of
+these intermediate quantities are available for output, along
+with ``Cdn_atm``, ``Cdn_ocn`` and the ratio ``Cdn_atm_ratio_n`` between the
+total atmospheric drag and the atmospheric neutral drag coefficient.
+
+We assume that the total neutral drag coefficients are thickness
+category independent, but through their dependance on the diagnostic
+variables described above, they vary both spatially and temporally. The
+total drag coefficients and heat transfer coefficients will also depend
+on the type of stratification of the atmosphere and the ocean, and we
+use the parameterization described in the :ref:`atmo` section that accounts
+for both stable and unstable atmosphere–ice boundary layers. In contrast
+to the neutral drag coefficients the stability effect of the atmospheric
+boundary layer is calculated separately for each ice thickness category.
+
+The transfer coefficient for oceanic heat flux to the bottom of the ice
+may be varied based on form drag considerations by setting the namelist
+variable ``fbot_xfer_type`` to ``Cdn_ocn``; this is recommended when using
+the form drag parameterization. The default value of the transfer
+coefficient is 0.006 (``fbot_xfer_type = ’constant’``).

From 56ccc779b516422bab56a7a15784e94f368c31a3 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 16 Sep 2021 14:58:11 +0000
Subject: [PATCH 14/20] Remove "ccsm" as an option to atmbndy

---
 doc/source/user_guide/interfaces.include | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include
index 6cbe1ce34..e0a39dc4d 100644
--- a/doc/source/user_guide/interfaces.include
+++ b/doc/source/user_guide/interfaces.include
@@ -944,7 +944,7 @@ icepack_init_parameters
            TTTocn_in        ! for qsat over ocn
   
         character (len=char_len), intent(in), optional :: &
-           atmbndy_in ! atmo boundary method, 'similarity' ('ccsm3') 'constant' or 'mixed'
+           atmbndy_in ! atmo boundary method, 'similarity', 'constant' or 'mixed'
           
         logical (kind=log_kind), intent(in), optional :: &
            calc_strair_in, &  ! if true, calculate wind stress components
@@ -1306,7 +1306,7 @@ icepack_query_parameters
            TTTocn_out        ! for qsat over ocn
   
         character (len=char_len), intent(out), optional :: &
-           atmbndy_out ! atmo boundary method, 'similarity' ('ccsm3'), 'constant' or 'mixed'
+           atmbndy_out ! atmo boundary method, 'similarity', 'constant' or 'mixed'
           
         logical (kind=log_kind), intent(out), optional :: &
            calc_strair_out, &  ! if true, calculate wind stress components

From 4807cc49982c401f82ddfca660e08910c5943bb4 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Fri, 17 Sep 2021 11:13:01 +0000
Subject: [PATCH 15/20] atmbndy: 'constant' or 'mixed'. Anything else end in
 'similarity' behaviour

---
 columnphysics/icepack_atmo.F90 | 11 +----------
 1 file changed, 1 insertion(+), 10 deletions(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index 49e15f8c5..ec67d58a2 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -950,10 +950,6 @@ subroutine icepack_atm_boundary(sfctype,                   &
 
       Cdn_atm_ratio_n = c1
 
-      if (trim(atmbndy) == 'default') then
-         atmbndy = 'similarity' ! For backward compability only ...
-      endif
-
       if (trim(atmbndy) == 'constant') then
          call atmo_boundary_const (sfctype,  calc_strair, &
                                    uatm,     vatm,     &
@@ -964,8 +960,7 @@ subroutine icepack_atm_boundary(sfctype,                   &
                                    delt,     delq,     &
                                    lhcoef,   shcoef    )
          if (icepack_warnings_aborted(subname)) return
-      elseif (trim(atmbndy) == 'similarity' .or. &
-              trim(atmbndy) == 'mixed') then
+      else
          call atmo_boundary_layer (sfctype,                 &
                                    calc_strair, formdrag,   &
                                    Tsf,      potT,          &
@@ -984,10 +979,6 @@ subroutine icepack_atm_boundary(sfctype,                   &
                                    uvel=l_uvel, vvel=l_vvel,  &
                                    Uref=l_Uref, zlvs=zlvs)
          if (icepack_warnings_aborted(subname)) return
-      else
-         call icepack_warnings_add( &
-            subname//' atmbndy = '//trim(atmbndy)//' : unknown value')
-         call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
       endif ! atmbndy
 
       if (present(Uref)) then

From 78660ed2c58ee021eabae528839f57d98b7538c9 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Mon, 20 Sep 2021 09:09:44 +0000
Subject: [PATCH 16/20] Error correction. Revert atmbndy to constant as before

---
 configuration/scripts/options/set_nml.dyn | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/configuration/scripts/options/set_nml.dyn b/configuration/scripts/options/set_nml.dyn
index 623438a10..f95f0c60e 100644
--- a/configuration/scripts/options/set_nml.dyn
+++ b/configuration/scripts/options/set_nml.dyn
@@ -3,5 +3,5 @@ kstrength = 0
 krdg_partic = 0
 krdg_redist = 0
 formdrag = .false.
-atmbndy = 'similarity'
+atmbndy = 'constant'
 highfreq = .false.

From 1ec0b3f23306b25a799f0f6e2d1676f64f76645f Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Mon, 20 Sep 2021 09:12:43 +0000
Subject: [PATCH 17/20] Put constants intp icepack_parameters

---
 columnphysics/icepack_atmo.F90       | 9 +++++----
 columnphysics/icepack_parameters.F90 | 2 ++
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index ec67d58a2..4c410a3e9 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -15,6 +15,7 @@ module icepack_atmo
       use icepack_kinds
       use icepack_parameters,  only: c0, c1, c2, c4, c5, c8, c10
       use icepack_parameters,  only: c16, c20, p001, p01, p2, p4, p5, p75, puny
+      use icepack_parameters,  only: p001, p0012
       use icepack_parameters,  only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
       use icepack_parameters,  only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
       use icepack_parameters,  only: pih, dragio, rhoi, rhos, rhow
@@ -362,8 +363,8 @@ subroutine atmo_boundary_layer (sfctype,            &
       if (trim(atmbndy) == 'mixed') then
          !- Use constant coefficients for sensible and latent heat fluxes
          !    similar to atmo_boundary_const but using vmag instead of wind
-         shcoef = (1.20e-3_dbl_kind)*cp_air*rhoa*vmag
-         lhcoef = (1.50e-3_dbl_kind)*Lheat *rhoa*vmag
+         shcoef = p0012*cp_air*rhoa*vmag
+         lhcoef = p0015*Lheat *rhoa*vmag
       else ! 'similarity'
          !- Monin-Obukhov similarity theory for boundary layer
          shcoef = rhoa * ustar * cp * rh + c1
@@ -511,8 +512,8 @@ subroutine atmo_boundary_const (sfctype,  calc_strair, &
       ! coefficients for turbulent flux calculation
       !------------------------------------------------------------
 
-      shcoef = (1.20e-3_dbl_kind)*cp_air*rhoa*wind
-      lhcoef = (1.50e-3_dbl_kind)*Lheat *rhoa*wind
+      shcoef = p0012*cp_air*rhoa*wind
+      lhcoef = p0015*Lheat *rhoa*wind
 
       end subroutine atmo_boundary_const
 
diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index f3f81c253..36546c775 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -54,6 +54,8 @@ module icepack_parameters
          p15  = 0.15_dbl_kind, &
          p25  = 0.25_dbl_kind, &
          p75  = 0.75_dbl_kind, &
+         p0012= 0.0012_dbl_kind, &
+         p0015= 0.0015_dbl_kind, &
          p333 = c1/c3, &
          p666 = c2/c3, &
          spval_const= -1.0e36_dbl_kind

From bee1627eb1c5f231594b3cd7db0dca8b8789902f Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Mon, 20 Sep 2021 09:36:31 +0000
Subject: [PATCH 18/20] Introduce p0012, p0015 from icepack_parameters

---
 columnphysics/icepack_atmo.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index 4c410a3e9..41e1068e4 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -15,7 +15,7 @@ module icepack_atmo
       use icepack_kinds
       use icepack_parameters,  only: c0, c1, c2, c4, c5, c8, c10
       use icepack_parameters,  only: c16, c20, p001, p01, p2, p4, p5, p75, puny
-      use icepack_parameters,  only: p001, p0012
+      use icepack_parameters,  only: p0012, p0015
       use icepack_parameters,  only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
       use icepack_parameters,  only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
       use icepack_parameters,  only: pih, dragio, rhoi, rhos, rhow

From 2640405757969e784bf65b4b318fe4a32cc0c498 Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 23 Sep 2021 07:28:51 +0000
Subject: [PATCH 19/20] Changed p0012/p0015 to specific constants
 senscoef/latncoef

---
 columnphysics/icepack_atmo.F90       | 10 +++++-----
 columnphysics/icepack_parameters.F90 |  6 +++---
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index 41e1068e4..b18726f63 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/columnphysics/icepack_atmo.F90
@@ -15,7 +15,7 @@ module icepack_atmo
       use icepack_kinds
       use icepack_parameters,  only: c0, c1, c2, c4, c5, c8, c10
       use icepack_parameters,  only: c16, c20, p001, p01, p2, p4, p5, p75, puny
-      use icepack_parameters,  only: p0012, p0015
+      use icepack_parameters,  only: senscoef, latncoef
       use icepack_parameters,  only: cp_wv, cp_air, iceruf, zref, qqqice, TTTice, qqqocn, TTTocn
       use icepack_parameters,  only: Lsub, Lvap, vonkar, Tffresh, zvir, gravit
       use icepack_parameters,  only: pih, dragio, rhoi, rhos, rhow
@@ -363,8 +363,8 @@ subroutine atmo_boundary_layer (sfctype,            &
       if (trim(atmbndy) == 'mixed') then
          !- Use constant coefficients for sensible and latent heat fluxes
          !    similar to atmo_boundary_const but using vmag instead of wind
-         shcoef = p0012*cp_air*rhoa*vmag
-         lhcoef = p0015*Lheat *rhoa*vmag
+         shcoef = senscoef*cp_air*rhoa*vmag
+         lhcoef = latncoef*Lheat *rhoa*vmag
       else ! 'similarity'
          !- Monin-Obukhov similarity theory for boundary layer
          shcoef = rhoa * ustar * cp * rh + c1
@@ -512,8 +512,8 @@ subroutine atmo_boundary_const (sfctype,  calc_strair, &
       ! coefficients for turbulent flux calculation
       !------------------------------------------------------------
 
-      shcoef = p0012*cp_air*rhoa*wind
-      lhcoef = p0015*Lheat *rhoa*wind
+      shcoef = senscoef*cp_air*rhoa*wind
+      lhcoef = latncoef*Lheat *rhoa*wind
 
       end subroutine atmo_boundary_const
 
diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index 36546c775..46759d83e 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -54,8 +54,6 @@ module icepack_parameters
          p15  = 0.15_dbl_kind, &
          p25  = 0.25_dbl_kind, &
          p75  = 0.75_dbl_kind, &
-         p0012= 0.0012_dbl_kind, &
-         p0015= 0.0015_dbl_kind, &
          p333 = c1/c3, &
          p666 = c2/c3, &
          spval_const= -1.0e36_dbl_kind
@@ -249,7 +247,9 @@ module icepack_parameters
          qqqice = 11637800._dbl_kind ,&! for qsat over ice
          TTTice = 5897.8_dbl_kind    ,&! for qsat over ice
          qqqocn = 627572.4_dbl_kind  ,&! for qsat over ocn
-         TTTocn = 5107.4_dbl_kind      ! for qsat over ocn
+         TTTocn = 5107.4_dbl_kind    ,&! for qsat over ocn
+         senscoef= 0.0012_dbl_kind   ,&! Sensible heat flux coefficient for constant-based boundary layer
+         latncoef= 0.0015_dbl_kind   , ! Latent heat flux coefficient for constant-based boundary layer
 
       character (len=char_len), public :: &
          atmbndy = 'similarity'        ! atmo boundary method, 'similarity', 'constant' or 'mixed'

From 81b80baa94f3268fc066573dc19222f5762b2f7b Mon Sep 17 00:00:00 2001
From: mhrib <mhri@dmi.dk>
Date: Thu, 23 Sep 2021 07:43:43 +0000
Subject: [PATCH 20/20] senscoef/latncoef introduced

---
 columnphysics/icepack_parameters.F90 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/columnphysics/icepack_parameters.F90 b/columnphysics/icepack_parameters.F90
index 46759d83e..50aa3a47f 100644
--- a/columnphysics/icepack_parameters.F90
+++ b/columnphysics/icepack_parameters.F90
@@ -249,7 +249,7 @@ module icepack_parameters
          qqqocn = 627572.4_dbl_kind  ,&! for qsat over ocn
          TTTocn = 5107.4_dbl_kind    ,&! for qsat over ocn
          senscoef= 0.0012_dbl_kind   ,&! Sensible heat flux coefficient for constant-based boundary layer
-         latncoef= 0.0015_dbl_kind   , ! Latent heat flux coefficient for constant-based boundary layer
+         latncoef= 0.0015_dbl_kind     ! Latent heat flux coefficient for constant-based boundary layer
 
       character (len=char_len), public :: &
          atmbndy = 'similarity'        ! atmo boundary method, 'similarity', 'constant' or 'mixed'