Update Derecho Cray, rename subroutines, minor cleanup (CICE-Consorti…

…um#502)

* Update Derecho Cray, rename subroutines, minor cleanup

Update Derecho Cray, add -Rp to -O2 standard optimization to address SAME_TBS compiler bug that is triggered in limited cases. Reported issue to NCAR. -O2 -Rp changes answers relative to -O2.

Rename two public Icepack subroutines, but maintain backward compatibility. Update the calls in Icepack driver, but also confirmed it works fine with the old names as well as in CICE with the old names.
Closes CICE-Consortium#255 

    use icepack_therm_shared  , only: icepack_init_thermo => icepack_init_salinity 
    use icepack_therm_shared  , only: icepack_init_trcr => icepack_init_enthalpy 

Add a check and abort for negative values in the sqrt in computation of Tin in function calculate_Tin_from_qin.
Closes CICE-Consortium#482

Refactor calls to icepack_aggregate to make them consistent. This was part of the testing for the Derecho Cray bug, and decided to keep the implementation.

Update comments associated with floeshape constant attribution. Change from Steele to Rothrock 1984.
Closes CICE-Consortium#479

Clean up some of the variable declarations in subroutine set_state_var and module icedrv_state, merge multiple lines to one line and shift to assumed shape arrays where appropriate.

Clean up implementation error in icedrv_restart.F90, subroutine restartfile. This subroutine was using a parameter, ntrcr, directly from icepack_tracers. Switched that to a call to icepack_query_tracer_sizes.

Update documentation of kice noting it's use with BL99 and MU71. See CICE-Consortium#447.

Generate updated interface documentation (./icepack.setup --docintfc)

columnphysics/icepack_intfc.F90

@@ -112,9 +112,12 @@ module icepack_intfc
       use icepack_therm_shared  , only: icepack_snow_temperature
       use icepack_therm_shared  , only: icepack_liquidus_temperature
       use icepack_therm_shared  , only: icepack_sea_freezing_temperature
-      use icepack_therm_shared  , only: icepack_init_thermo
+      use icepack_therm_shared  , only: icepack_init_salinity
       use icepack_therm_shared  , only: icepack_salinity_profile
-      use icepack_therm_shared  , only: icepack_init_trcr
+      use icepack_therm_shared  , only: icepack_init_enthalpy
+      ! for backwards compatibilty, remove in the future
+      use icepack_therm_shared  , only: icepack_init_thermo => icepack_init_salinity
+      use icepack_therm_shared  , only: icepack_init_trcr => icepack_init_enthalpy
 
       use icepack_mushy_physics , only: icepack_enthalpy_snow
       use icepack_mushy_physics , only: icepack_enthalpy_mush

columnphysics/icepack_parameters.F90

@@ -117,8 +117,8 @@ module icepack_parameters
          Timelt    = 0.0_dbl_kind     ,&! melting temperature, ice top surface  (C)
          Tsmelt    = 0.0_dbl_kind     ,&! melting temperature, snow top surface (C)
          ice_ref_salinity =4._dbl_kind,&! (ppt)
-                                        ! kice is not used for mushy thermo
          kice      = 2.03_dbl_kind    ,&! thermal conductivity of fresh ice(W/m/deg)
+                                        ! kice is only used with ktherm=1 (BL99) and conduct='MU71'
          ksno      = 0.30_dbl_kind    ,&! thermal conductivity of snow  (W/m/deg)
          hs_min    = 1.e-4_dbl_kind   ,&! min snow thickness for computing zTsn (m)
          snowpatch = 0.02_dbl_kind    ,&! parameter for fractional snow area (m)
@@ -316,7 +316,7 @@ module icepack_parameters
          nfreq = 25                   ! number of frequencies
 
       real (kind=dbl_kind), public :: &
-         floeshape = 0.66_dbl_kind    ! constant from Steele (unitless)
+         floeshape = 0.66_dbl_kind    ! constant from Rothrock 1984 (unitless)
 
       real (kind=dbl_kind), public :: &
          floediam  = 300.0_dbl_kind   ! effective floe diameter for lateral melt (m)
@@ -856,7 +856,7 @@ subroutine icepack_init_parameters(   &
          nfreq_in           ! number of frequencies
 
       real (kind=dbl_kind), intent(in), optional :: &
-         floeshape_in       ! constant from Steele (unitless)
+         floeshape_in       ! constant from Rothrock 1984 (unitless)
 
       logical (kind=log_kind), intent(in), optional :: &
          wave_spec_in       ! if true, use wave forcing
@@ -1876,7 +1876,7 @@ subroutine icepack_query_parameters(   &
          nfreq_out          ! number of frequencies
 
       real (kind=dbl_kind), intent(out), optional :: &
-         floeshape_out      ! constant from Steele (unitless)
+         floeshape_out      ! constant from Rothrock 1984 (unitless)
 
       logical (kind=log_kind), intent(out), optional :: &
          wave_spec_out      ! if true, use wave forcing

columnphysics/icepack_therm_shared.F90

@@ -29,9 +29,9 @@ module icepack_therm_shared
       public :: calculate_Tin_from_qin, &
                 surface_heat_flux, &
                 dsurface_heat_flux_dTsf, &
-                icepack_init_thermo, &
+                icepack_init_salinity, &
                 icepack_salinity_profile, &
-                icepack_init_trcr, &
+                icepack_init_enthalpy, &
                 icepack_ice_temperature, &
                 icepack_snow_temperature, &
                 icepack_liquidus_temperature, &
@@ -65,21 +65,26 @@ function calculate_Tin_from_qin (qin, Tmltk) &
          Tmltk                  ! melting temperature at one level
 
       real (kind=dbl_kind) :: &
-         Tin                 ! internal temperature
+         Tin                    ! internal temperature
 
       ! local variables
 
       real (kind=dbl_kind) :: &
-         aa1,bb1,cc1         ! quadratic solvers
+         aa1,bb1,cc1,csqrt      ! quadratic solvers
 
       character(len=*),parameter :: subname='(calculate_Tin_from_qin)'
 
       if (l_brine) then
          aa1 = cp_ice
          bb1 = (cp_ocn-cp_ice)*Tmltk - qin/rhoi - Lfresh
          cc1 = Lfresh * Tmltk
-         Tin =  min((-bb1 - sqrt(bb1*bb1 - c4*aa1*cc1)) /  &
-                         (c2*aa1),Tmltk)
+         csqrt = bb1*bb1 - c4*aa1*cc1
+         if (csqrt < c0) then
+           call icepack_warnings_add(subname//' sqrt error: ')
+           call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+           return
+         endif
+         Tin =  min((-bb1 - sqrt(csqrt)) / (c2*aa1),Tmltk)
 
       else                ! fresh ice
          Tin = (Lfresh + qin/rhoi) / cp_ice
@@ -210,13 +215,14 @@ subroutine dsurface_heat_flux_dTsf(Tsf,  rhoa,      &
       end subroutine dsurface_heat_flux_dTsf
 
 !=======================================================================
-!autodocument_start icepack_init_thermo
-! Initialize the vertical profile of ice salinity and melting temperature.
+!autodocument_start icepack_init_salinity
+! Initialize the vertical profile of ice salinity.
+! This subroutine was renamed from icepack_init_thermo in Oct 2024
 !
 ! authors: C. M. Bitz, UW
 !          William H. Lipscomb, LANL
 
-      subroutine icepack_init_thermo(sprofile)
+      subroutine icepack_init_salinity(sprofile)
 
       real (kind=dbl_kind), dimension(:), intent(out) :: &
          sprofile                         ! vertical salinity profile
@@ -226,7 +232,7 @@ subroutine icepack_init_thermo(sprofile)
       integer (kind=int_kind) :: k        ! ice layer index
       real (kind=dbl_kind)    :: zn       ! normalized ice thickness
 
-      character(len=*),parameter :: subname='(icepack_init_thermo)'
+      character(len=*),parameter :: subname='(icepack_init_salinity)'
 
       !-----------------------------------------------------------------
       ! Determine l_brine based on saltmax.
@@ -239,7 +245,7 @@ subroutine icepack_init_thermo(sprofile)
       if (saltmax > min_salin) l_brine = .true.
 
       !-----------------------------------------------------------------
-      ! Prescibe vertical profile of salinity and melting temperature.
+      ! Prescibe vertical profile of salinity.
       ! Note this profile is only used for BL99 thermodynamics.
       !-----------------------------------------------------------------
 
@@ -259,7 +265,7 @@ subroutine icepack_init_thermo(sprofile)
          enddo
       endif ! l_brine
 
-      end subroutine icepack_init_thermo
+      end subroutine icepack_init_salinity
 
 !=======================================================================
 !autodocument_start icepack_salinity_profile
@@ -282,16 +288,17 @@ function icepack_salinity_profile(zn) result(salinity)
          nsal    = 0.407_dbl_kind, &
          msal    = 0.573_dbl_kind
 
-      character(len=*),parameter :: subname='(icepack_init_thermo)'
+      character(len=*),parameter :: subname='(icepack_salinity_profile)'
 
       salinity = (saltmax/c2)*(c1-cos(pi*zn**(nsal/(msal+zn))))
 
       end function icepack_salinity_profile
 
 !=======================================================================
-!autodocument_start icepack_init_trcr
+!autodocument_start icepack_init_enthalpy
+! This subroutine was renamed from icepack_init_trcr in Oct 2024
 !
-      subroutine icepack_init_trcr(Tair,     Tf,       &
+      subroutine icepack_init_enthalpy(Tair, Tf,      &
                                   Sprofile, Tprofile, &
                                   Tsfc,               &
                                   qin,      qsn)
@@ -320,7 +327,7 @@ subroutine icepack_init_trcr(Tair,     Tf,       &
       real (kind=dbl_kind) :: &
          slope, Ti
 
-      character(len=*),parameter :: subname='(icepack_init_trcr)'
+      character(len=*),parameter :: subname='(icepack_init_enthalpy)'
 
       ! surface temperature
       Tsfc = Tf ! default
@@ -346,7 +353,7 @@ subroutine icepack_init_trcr(Tair,     Tf,       &
           qsn(k) = -rhos*(Lfresh - cp_ice*Ti)
         enddo               ! nslyr
 
-    end subroutine icepack_init_trcr
+    end subroutine icepack_init_enthalpy
 
 !=======================================================================
 !autodocument_start icepack_liquidus_temperature
@@ -406,6 +413,7 @@ function icepack_sea_freezing_temperature(sss) result(Tf)
 
            call icepack_warnings_add(subname//' tfrz_option unsupported: '//trim(tfrz_option))
            call icepack_warnings_setabort(.true.,__FILE__,__LINE__)
+           return
 
         endif
 

configuration/driver/icedrv_InitMod.F90

@@ -11,6 +11,7 @@ module icedrv_InitMod
       use icedrv_constants, only: nu_diag
       use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
       use icepack_intfc, only: icepack_query_parameters, icepack_query_tracer_flags
+      use icepack_intfc, only: icepack_query_tracer_sizes
       use icepack_intfc, only: icepack_write_tracer_flags, icepack_write_tracer_indices
       use icepack_intfc, only: icepack_write_tracer_sizes, icepack_write_parameters
       use icedrv_system, only: icedrv_system_abort, icedrv_system_flush
@@ -187,12 +188,13 @@ subroutine init_restart
       use icedrv_state ! almost everything
 
       integer(kind=int_kind) :: &
-         i                            ! horizontal indices
+         i,          & ! horizontal indices
+         ntrcr         ! tracer count
 
       logical (kind=log_kind) :: &
-         skl_bgc, &    ! from icepack
-         z_tracers, &  ! from icepack
-         tr_brine, &   ! from icepack
+         skl_bgc,    & ! from icepack
+         z_tracers,  & ! from icepack
+         tr_brine,   & ! from icepack
          tr_fsd        ! from icepack
 
       character(len=*), parameter :: subname='(init_restart)'
@@ -204,6 +206,7 @@ subroutine init_restart
       call icepack_query_parameters(skl_bgc_out=skl_bgc)
       call icepack_query_parameters(z_tracers_out=z_tracers)
       call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_fsd_out=tr_fsd)
+      call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__,line= __LINE__)
@@ -235,20 +238,20 @@ subroutine init_restart
       !-----------------------------------------------------------------
       do i = 1, nx
          if (tmask(i)) &
-         call icepack_aggregate(aicen=aicen(i,:),   &
-                                vicen=vicen(i,:),   &
-                                vsnon=vsnon(i,:),   &
-                                trcrn=trcrn(i,:,:), &
-                                aice=aice (i),      &
-                                vice=vice (i),      &
-                                vsno=vsno (i),      &
-                                trcr=trcr (i,:),    &
-                                aice0=aice0(i),     &
-                                trcr_depend=trcr_depend, &
-                                trcr_base=trcr_base,     &
-                                n_trcr_strata=n_trcr_strata, &
-                                nt_strata=nt_strata, &
-                                Tf=Tf(i))
+         call icepack_aggregate(trcrn=trcrn(i,1:ntrcr,:),     &
+                                aicen=aicen(i,:),             &
+                                vicen=vicen(i,:),             &
+                                vsnon=vsnon(i,:),             &
+                                trcr=trcr (i,1:ntrcr),        &
+                                aice=aice (i),                &
+                                vice=vice (i),                &
+                                vsno=vsno (i),                &
+                                aice0=aice0(i),               &
+                                trcr_depend=trcr_depend(1:ntrcr),     &
+                                trcr_base=trcr_base    (1:ntrcr,:),   &
+                                n_trcr_strata=n_trcr_strata(1:ntrcr), &
+                                nt_strata=nt_strata    (1:ntrcr,:), &
+                                Tf = Tf(i))
       enddo
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &

configuration/driver/icedrv_init.F90

@@ -15,7 +15,7 @@ module icedrv_init
       use icepack_intfc, only: icepack_init_tracer_flags
       use icepack_intfc, only: icepack_init_tracer_sizes
       use icepack_intfc, only: icepack_init_tracer_indices
-      use icepack_intfc, only: icepack_init_trcr
+      use icepack_intfc, only: icepack_init_enthalpy
       use icepack_intfc, only: icepack_query_parameters
       use icepack_intfc, only: icepack_query_tracer_flags
       use icepack_intfc, only: icepack_query_tracer_sizes
@@ -1301,27 +1301,24 @@ subroutine set_state_var (nx, &
       integer (kind=int_kind), intent(in) :: &
          nx          ! number of grid cells
 
-      real (kind=dbl_kind), dimension (nx), intent(in) :: &
+      real (kind=dbl_kind), dimension (:), intent(in) :: &
          Tair       ! air temperature  (K)
 
       ! ocean values may be redefined here, unlike in CICE
-      real (kind=dbl_kind), dimension (nx), intent(inout) :: &
+      real (kind=dbl_kind), dimension (:), intent(inout) :: &
          Tf     , & ! freezing temperature (C)
          sst        ! sea surface temperature (C)
 
-      real (kind=dbl_kind), dimension (nx,nilyr), &
-         intent(in) :: &
+      real (kind=dbl_kind), dimension (:,:), intent(in) :: &
          salinz , & ! initial salinity profile
          Tmltz      ! initial melting temperature profile
 
-      real (kind=dbl_kind), dimension (nx,ncat), &
-         intent(out) :: &
+      real (kind=dbl_kind), dimension (:,:), intent(out) :: &
          aicen , & ! concentration of ice
          vicen , & ! volume per unit area of ice          (m)
          vsnon     ! volume per unit area of snow         (m)
 
-      real (kind=dbl_kind), dimension (nx,max_ntrcr,ncat), &
-         intent(out) :: &
+      real (kind=dbl_kind), dimension (:,:,:), intent(out) :: &
          trcrn     ! ice tracers
                    ! 1: surface temperature of ice/snow (C)
 
@@ -1431,7 +1428,7 @@ subroutine set_state_var (nx, &
          vicen(i,n) = hinit(n) * ainit(n) ! m
          vsnon(i,n) = c0
          ! tracers
-         call icepack_init_trcr(Tair     = Tair(i),     &
+         call icepack_init_enthalpy(Tair = Tair(i),     &
                                 Tf       = Tf(i),       &
                                 Sprofile = salinz(i,:), &
                                 Tprofile = Tmltz(i,:),  &
@@ -1502,7 +1499,7 @@ subroutine set_state_var (nx, &
          vicen(i,n) = hinit(n) * ainit(n) ! m
          vsnon(i,n) = min(aicen(i,n)*hsno_init,p2*vicen(i,n))
          ! tracers
-         call icepack_init_trcr(Tair     = Tair(i),     &
+         call icepack_init_enthalpy(Tair = Tair(i),     &
                                 Tf       = Tf(i),       &
                                 Sprofile = salinz(i,:), &
                                 Tprofile = Tmltz(i,:),  &

configuration/driver/icedrv_init_column.F90

@@ -24,7 +24,7 @@ module icedrv_init_column
       use icepack_intfc, only: icepack_query_tracer_indices
       use icepack_intfc, only: icepack_query_parameters
       use icepack_intfc, only: icepack_init_zbgc
-      use icepack_intfc, only: icepack_init_thermo, icepack_init_radiation
+      use icepack_intfc, only: icepack_init_salinity, icepack_init_radiation
       use icepack_intfc, only: icepack_step_radiation, icepack_init_orbit
       use icepack_intfc, only: icepack_init_bgc
       use icepack_intfc, only: icepack_init_ocean_bio, icepack_load_ocean_bio_array
@@ -69,7 +69,7 @@ subroutine init_thermo_vertical
       !-----------------------------------------------------------------
 
       call icepack_query_parameters(depressT_out=depressT)
-      call icepack_init_thermo(sprofile=sprofile)
+      call icepack_init_salinity(sprofile=sprofile)
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__, line=__LINE__)