Update to Consortium #795280797a Sept 13, 2023

columnphysics/icepack_intfc.F90

@@ -1,8 +1,8 @@
 !=======================================================================
-! Copyright (c) 2022, Triad National Security, LLC
+! Copyright (c) 2023, Triad National Security, LLC
 ! All rights reserved.
 !
-! Copyright 2022. Triad National Security, LLC. This software was
+! Copyright 2023. Triad National Security, LLC. This software was
 ! produced under U.S. Government contract DE-AC52-06NA25396 for Los
 ! Alamos National Laboratory (LANL), which is operated by Triad
 ! National Security, LLC for the U.S. Department of Energy. The U.S.

columnphysics/icepack_parameters.F90

@@ -146,8 +146,9 @@ module icepack_parameters
                          ! 2 = mushy layer theory
 
       character (char_len), public :: &
-         conduct = 'bubbly', &      ! 'MU71' or 'bubbly'
-         fbot_xfer_type = 'constant' ! transfer coefficient type for ice-ocean heat flux
+         conduct = 'bubbly', &          ! 'MU71' or 'bubbly'
+         fbot_xfer_type = 'constant', & ! transfer coefficient type for ice-ocean heat flux
+         cpl_frazil = 'fresh_ice_correction' ! type of coupling for frazil ice
 
       logical (kind=log_kind), public :: &
          calc_Tsfc     = .true. ,&! if true, calculate surface temperature
@@ -438,7 +439,8 @@ subroutine icepack_init_parameters(   &
          argcheck_in, puny_in, bignum_in, pi_in, secday_in, &
          rhos_in, rhoi_in, rhow_in, cp_air_in, emissivity_in, &
          cp_ice_in, cp_ocn_in, hfrazilmin_in, floediam_in, &
-         depressT_in, dragio_in, thickness_ocn_layer1_in, iceruf_ocn_in, albocn_in, gravit_in, viscosity_dyn_in, &
+         depressT_in, dragio_in, thickness_ocn_layer1_in, iceruf_ocn_in, &
+         albocn_in, gravit_in, viscosity_dyn_in, &
          Tocnfrz_in, rhofresh_in, zvir_in, vonkar_in, cp_wv_in, &
          stefan_boltzmann_in, ice_ref_salinity_in, &
          Tffresh_in, Lsub_in, Lvap_in, Timelt_in, Tsmelt_in, &
@@ -452,6 +454,7 @@ subroutine icepack_init_parameters(   &
          qqqice_in, TTTice_in, qqqocn_in, TTTocn_in, &
          ktherm_in, conduct_in, fbot_xfer_type_in, calc_Tsfc_in, dts_b_in, &
          update_ocn_f_in, ustar_min_in, hi_min_in, a_rapid_mode_in, &
+         cpl_frazil_in, &
          Rac_rapid_mode_in, aspect_rapid_mode_in, &
          dSdt_slow_mode_in, phi_c_slow_mode_in, &
          phi_i_mushy_in, shortwave_in, albedo_type_in, albsnowi_in, &
@@ -545,8 +548,9 @@ subroutine icepack_init_parameters(   &
                             ! 2 = mushy layer theory
 
       character (len=*), intent(in), optional :: &
-         conduct_in, &      ! 'MU71' or 'bubbly'
-         fbot_xfer_type_in  ! transfer coefficient type for ice-ocean heat flux
+         conduct_in, &        ! 'MU71' or 'bubbly'
+         fbot_xfer_type_in, & ! transfer coefficient type for ice-ocean heat flux
+         cpl_frazil_in        ! type of coupling for frazil ice
 
       logical (kind=log_kind), intent(in), optional :: &
          calc_Tsfc_in    , &! if true, calculate surface temperature
@@ -908,6 +912,7 @@ subroutine icepack_init_parameters(   &
       if (present(conduct_in)           ) conduct          = conduct_in
       if (present(fbot_xfer_type_in)    ) fbot_xfer_type   = fbot_xfer_type_in
       if (present(calc_Tsfc_in)         ) calc_Tsfc        = calc_Tsfc_in
+      if (present(cpl_frazil_in)        ) cpl_frazil       = cpl_frazil_in
       if (present(update_ocn_f_in)      ) update_ocn_f     = update_ocn_f_in
       if (present(dts_b_in)             ) dts_b            = dts_b_in
       if (present(ustar_min_in)         ) ustar_min        = ustar_min_in
@@ -1170,7 +1175,7 @@ subroutine icepack_query_parameters(   &
          saltmax_out, phi_init_out, min_salin_out, salt_loss_out, &
          Tliquidus_max_out, &
          min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, &
-         awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, &
+         awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, cpl_frazil_out, &
          qqqice_out, TTTice_out, qqqocn_out, TTTocn_out, update_ocn_f_out, &
          Lfresh_out, cprho_out, Cp_out, ustar_min_out, hi_min_out, a_rapid_mode_out, &
          ktherm_out, conduct_out, fbot_xfer_type_out, calc_Tsfc_out, dts_b_out, &
@@ -1276,8 +1281,9 @@ subroutine icepack_query_parameters(   &
                             ! 2 = mushy layer theory
 
       character (len=*), intent(out), optional :: &
-         conduct_out, &     ! 'MU71' or 'bubbly'
-         fbot_xfer_type_out ! transfer coefficient type for ice-ocean heat flux
+         conduct_out, &        ! 'MU71' or 'bubbly'
+         fbot_xfer_type_out, & ! transfer coefficient type for ice-ocean heat flux
+         cpl_frazil_out        ! type of coupling for frazil ice
 
       logical (kind=log_kind), intent(out), optional :: &
          calc_Tsfc_out    ,&! if true, calculate surface temperature
@@ -1673,6 +1679,7 @@ subroutine icepack_query_parameters(   &
       if (present(conduct_out)           ) conduct_out      = conduct
       if (present(fbot_xfer_type_out)    ) fbot_xfer_type_out = fbot_xfer_type
       if (present(calc_Tsfc_out)         ) calc_Tsfc_out    = calc_Tsfc
+      if (present(cpl_frazil_out)        ) cpl_frazil_out   = cpl_frazil
       if (present(update_ocn_f_out)      ) update_ocn_f_out = update_ocn_f
       if (present(dts_b_out)             ) dts_b_out        = dts_b
       if (present(ustar_min_out)         ) ustar_min_out    = ustar_min
@@ -1881,6 +1888,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  conduct    = ", trim(conduct)
         write(iounit,*) "  fbot_xfer_type = ", trim(fbot_xfer_type)
         write(iounit,*) "  calc_Tsfc  = ", calc_Tsfc
+        write(iounit,*) "  cpl_frazil = ", cpl_frazil
         write(iounit,*) "  update_ocn_f = ", update_ocn_f
         write(iounit,*) "  dts_b      = ", dts_b
         write(iounit,*) "  ustar_min  = ", ustar_min

columnphysics/icepack_therm_itd.F90

@@ -28,7 +28,7 @@ module icepack_therm_itd
       use icepack_parameters, only: rhosi, conserv_check, rhosmin, snwredist
       use icepack_parameters, only: kitd, ktherm
       use icepack_parameters, only: z_tracers, hfrazilmin, hi_min
-      use icepack_parameters, only: saltflux_option
+      use icepack_parameters, only: cpl_frazil, update_ocn_f, saltflux_option
       use icepack_parameters, only: icepack_chkoptargflag
 
       use icepack_tracers, only: ntrcr, nbtrcr
@@ -1308,7 +1308,6 @@ subroutine add_new_ice (ncat,      nilyr,      &
                               aice0,     aice,       &
                               frzmlt,    frazil,     &
                               frz_onset, yday,       &
-                              update_ocn_f,          &
                               fresh,     fsalt,      &
                               Tf,        sss,        &
                               salinz,    phi_init,   &
@@ -1381,9 +1380,6 @@ subroutine add_new_ice (ncat,      nilyr,      &
          phi_init     , & ! initial frazil liquid fraction
          dSin0_frazil     ! initial frazil bulk salinity reduction from sss
 
-      logical (kind=log_kind), intent(in) :: &
-         update_ocn_f ! if true, update fresh water and salt fluxes
-
       ! BGC
       real (kind=dbl_kind), dimension (nblyr+2), intent(in) :: &
          bgrid              ! biology nondimensional vertical grid points
@@ -1619,29 +1615,29 @@ subroutine add_new_ice (ncat,      nilyr,      &
       !       is NOT included in fluxes fresh and fsalt.
       !-----------------------------------------------------------------
 
-      if (update_ocn_f) then
-         dfresh = -rhoi*vi0new/dt
+      dfresh = c0
+      dfsalt = c0
+      if (cpl_frazil == 'external') then
+         ! do nothing here, calculations are in the coupler or elsewhere
+      else
+         if (update_ocn_f) then
+            dfresh = -rhoi*vi0new/dt
+         elseif (cpl_frazil == 'fresh_ice_correction' .and. ktherm == 2) then
+            ! correct frazil fluxes for mushy
+            vi0tmp = fnew*dt / (rhoi*Lfresh) ! ocn/cpl assumes frazil volume is pure, fresh ice
+            dfresh = -rhoi*(vi0new - vi0tmp)/dt
+            frazil_diag = frazil - vi0tmp
+!        else 
+!           do nothing - other correction options could be implemented in the future             
+         endif
+
          if (saltflux_option == 'prognostic') then
             dfsalt = Si0new*p001*dfresh
          else
             dfsalt = ice_ref_salinity*p001*dfresh
          endif
          fresh  = fresh + dfresh
-         fsalt  = fsalt + dfsalt
-      else ! update_ocn_f = false
-         if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
-            vi0tmp = fnew*dt / (rhoi*Lfresh)
-            dfresh = -rhoi*(vi0new - vi0tmp)/dt
-            if (saltflux_option == 'prognostic') then
-               dfsalt = Si0new*p001*dfresh
-            else
-               dfsalt = ice_ref_salinity*p001*dfresh
-            endif
-            fresh  = fresh + dfresh
-            fsalt  = fsalt + dfsalt
-            frazil_diag = frazil - vi0tmp
-         ! elseif ktherm==1 do nothing
-         endif
+         fsalt  = fsalt + dfsalt      
       endif
 
       !-----------------------------------------------------------------
@@ -2011,6 +2007,8 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
                                      d_afsd_latm,  d_afsd_weld,   &
                                      floe_rad_c,   floe_binwidth)
 
+      use icepack_parameters, only: icepack_init_parameters
+
       integer (kind=int_kind), intent(in) :: &
          ncat     , & ! number of thickness categories
          nltrcr   , & ! number of zbgc tracers
@@ -2021,7 +2019,7 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
       integer (kind=int_kind), intent(in), optional :: &
          nfsd         ! number of floe size categories
 
-      logical (kind=log_kind), intent(in) :: &
+      logical (kind=log_kind), intent(in), optional :: &
          update_ocn_f ! if true, update fresh water and salt fluxes
 
       real (kind=dbl_kind), dimension(0:ncat), intent(in) :: &
@@ -2141,6 +2139,9 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
       ! Check optional arguments and set local values
       !-----------------------------------------------------------------
 
+       if (present(update_ocn_f)) then
+          call icepack_init_parameters(update_ocn_f_in=update_ocn_f)
+       endif
        if (icepack_chkoptargflag(first_call)) then
           if (tr_iso) then
              if (.not.(present(fiso_ocn)   .and. &
@@ -2239,7 +2240,6 @@ subroutine icepack_step_therm2 (dt, ncat, nltrcr,           &
                            aice0,         aice,         &
                            frzmlt,        frazil,       &
                            frz_onset,     yday,         &
-                           update_ocn_f,                &
                            fresh,         fsalt,        &
                            Tf,            sss,          &
                            salinz,        phi_init,     &

columnphysics/version.txt

@@ -1 +1 @@
-ICEPACK 1.3.3
+ICEPACK 1.3.4

configuration/driver/icedrv_MAIN.F90

@@ -1,8 +1,8 @@
 !=======================================================================
-! Copyright (c) 2022, Triad National Security, LLC
+! Copyright (c) 2023, Triad National Security, LLC
 ! All rights reserved.
 !
-! Copyright 2022. Triad National Security, LLC. This software was
+! Copyright 2023. Triad National Security, LLC. This software was
 ! produced under U.S. Government contract DE-AC52-06NA25396 for Los
 ! Alamos National Laboratory (LANL), which is operated by Triad
 ! National Security, LLC for the U.S. Department of Energy. The U.S.

configuration/driver/icedrv_flux.F90

@@ -19,7 +19,6 @@ module icedrv_flux
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_query_parameters
       use icepack_intfc, only: icepack_query_tracer_flags, icepack_query_tracer_indices
-      use icepack_intfc, only: icepack_query_parameters
       use icedrv_system, only: icedrv_system_abort
 
       implicit none
@@ -193,7 +192,6 @@ module icedrv_flux
          scale_factor! scaling factor for shortwave components
 
       logical (kind=log_kind), public :: &
-         update_ocn_f, & ! if true, update fresh water and salt fluxes
          l_mpond_fresh   ! if true, include freshwater feedback from meltponds
                          ! when running in ice-ocean or coupled configuration
 

configuration/driver/icedrv_init.F90

@@ -64,7 +64,7 @@ subroutine input_data
       use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years
       use icedrv_history, only: history_format
       use icedrv_restart_shared, only: restart, restart_dir, restart_file, restart_format
-      use icedrv_flux, only: update_ocn_f, l_mpond_fresh, cpl_bgc
+      use icedrv_flux, only: l_mpond_fresh, cpl_bgc
       use icedrv_flux, only: default_season
       use icedrv_forcing, only: precip_units,    fyear_init,      ycycle
       use icedrv_forcing, only: atm_data_type,   ocn_data_type,   bgc_data_type
@@ -99,9 +99,10 @@ subroutine input_data
          natmiter, kitd, kcatbound
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
-         tfrz_option, saltflux_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
+         cpl_frazil, tfrz_option, saltflux_option, &
+         frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
 
-      logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain
+      logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain, update_ocn_f
       real (kind=dbl_kind)    :: sw_frac, sw_dtemp
 
       ! Flux convergence tolerance
@@ -175,7 +176,7 @@ subroutine input_data
         formdrag,        highfreq,        natmiter,        &
         atmiter_conv,    calc_dragio,                      &
         tfrz_option,     saltflux_option, ice_ref_salinity, &
-        default_season,  wave_spec_type,  &
+        default_season,  wave_spec_type,  cpl_frazil,      &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
         lateral_flux_type,                                &
@@ -217,7 +218,7 @@ subroutine input_data
            pndaspect_out=pndaspect, hs1_out=hs1, hp1_out=hp1, &
            ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc, &
            floediam_out=floediam, hfrazilmin_out=hfrazilmin, &
-           update_ocn_f_out = update_ocn_f, &
+           update_ocn_f_out = update_ocn_f, cpl_frazil_out = cpl_frazil, &
            conduct_out=conduct, a_rapid_mode_out=a_rapid_mode, &
            Rac_rapid_mode_out=Rac_rapid_mode, &
            aspect_rapid_mode_out=aspect_rapid_mode, &
@@ -780,6 +781,7 @@ subroutine input_data
          if (trim(atm_data_type)=='default') &
          write(nu_diag,1030) ' default_season            = ', trim(default_season)
 
+         write(nu_diag,1030) ' cpl_frazil                = ', trim(cpl_frazil)
          write(nu_diag,1010) ' update_ocn_f              = ', update_ocn_f
          write(nu_diag,1010) ' wave_spec                 = ', wave_spec
          if (wave_spec) &
@@ -971,6 +973,7 @@ subroutine input_data
            floediam_in=floediam, hfrazilmin_in=hfrazilmin, &
            ktherm_in=ktherm, calc_Tsfc_in=calc_Tsfc, &
            conduct_in=conduct, a_rapid_mode_in=a_rapid_mode, &
+           update_ocn_f_in=update_ocn_f, cpl_frazil_in=cpl_frazil, &
            Rac_rapid_mode_in=Rac_rapid_mode, &
            aspect_rapid_mode_in=aspect_rapid_mode, &
            dSdt_slow_mode_in=dSdt_slow_mode, &

configuration/driver/icedrv_step.F90

@@ -437,7 +437,7 @@ subroutine step_therm2 (dt)
       use icedrv_domain_size, only: ncat, nilyr, nslyr, n_aero, nblyr, &
                                     nltrcr, nx, nfsd
       use icedrv_flux, only: fresh, frain, fpond, frzmlt, frazil, frz_onset
-      use icedrv_flux, only: update_ocn_f, fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat
+      use icedrv_flux, only: fsalt, Tf, sss, salinz, fhocn, rside, fside, wlat
       use icedrv_flux, only: meltl, frazil_diag, flux_bio, faero_ocn, fiso_ocn
       use icedrv_flux, only: HDO_ocn, H2_16O_ocn, H2_18O_ocn
       use icedrv_init, only: tmask
@@ -505,7 +505,7 @@ subroutine step_therm2 (dt)
                          frzmlt=frzmlt(i), frazil=frazil(i),          &
                          frain=frain(i),   fpond=fpond(i),            &
                          fresh=fresh(i),   fsalt=fsalt(i),            &
-                         fhocn=fhocn(i),   update_ocn_f=update_ocn_f, &
+                         fhocn=fhocn(i),                              &
                          bgrid=bgrid,      cgrid=cgrid,               &
                          igrid=igrid,      faero_ocn=faero_ocn(i,:),  &
                          first_ice=first_ice(i,:),                    &

configuration/scripts/icepack_in

@@ -110,6 +110,7 @@
     calc_dragio     = .false.
     emissivity      = 0.985
     fbot_xfer_type  = 'constant'
+    cpl_frazil      = 'fresh_ice_correction'
     update_ocn_f    = .false.
     l_mpond_fresh   = .false.
     tfrz_option     = 'mushy'

configuration/scripts/machines/Macros.derecho_cray

@@ -0,0 +1,47 @@
+#==============================================================================
+# Makefile macros for NCAR derecho, cray compiler
+#==============================================================================
+
+CPP        := ftn -e P
+CPPDEFS    := -DFORTRANUNDERSCORE -DNO_R16 ${ICE_CPPDEFS}
+CFLAGS     := -c -O2
+
+FIXEDFLAGS := -132
+FREEFLAGS  := 
+FFLAGS     := -hbyteswapio
+FFLAGS_NOOPT:= -O0
+
+ifeq ($(ICE_BLDDEBUG), true)
+  FFLAGS     += -O0 -hfp0 -g -Rbcdps -Ktrap=fp
+else
+  FFLAGS     += -O2 -hfp0   # -eo
+endif
+
+SCC   := cc
+SFC   := ftn
+CC := $(SCC)
+FC := $(SFC)
+LD := $(FC)
+
+NETCDF_PATH := $(NETCDF)
+
+#PIO_CONFIG_OPTS:= --enable-filesystem-hints=gpfs
+
+#PNETCDF_PATH := $(PNETCDF)
+#PNETCDF_PATH := /glade/u/apps/ch/opt/pio/2.2/mpt/2.15f/intel/17.0.1/lib
+
+INCLDIR := $(INCLDIR)
+
+LIB_NETCDF := $(NETCDF)/lib
+#LIB_PNETCDF := $(PNETCDF_PATH)/lib
+#LIB_MPI := $(IMPILIBDIR)
+
+#SLIBS   := -L$(LIB_NETCDF) -lnetcdf -lnetcdff -L$(LIB_PNETCDF) -lpnetcdf -lgptl
+SLIBS   := -L$(LIB_NETCDF) -lnetcdf -lnetcdff
+
+ifeq ($(ICE_THREADED), true)
+   LDFLAGS += -qopenmp
+   CFLAGS += -qopenmp
+   FFLAGS += -qopenmp
+endif
+