BGC update (CICE-Consortium#497)

This is a significant update in the BGC including refactoring Icepack interfaces.

Deprecate skl BGC but leave code alone for now hoping we get help from the community to validate the latest code.

Add BGC parameters to icepack_parameters.F90

Update BGC physics, consistent with E3SM-Project/E3SM#6457.  Remove redundant arguments in BGC interfaces.

    icepack_aerosol.F90
        revised subroutine update_snow_bgc
    icepack_algae.F90
        revised subroutine zbio
        revised subroutine z_biogeochemistry
        revised subroutine algal_dyn
        add subroutine bgc_carbon_sum
    icepack_brine.F90
        revise subroutine prepare_hbrine
        revise subroutine update_hbrine
    icepack_mechred.F90
        add mbio calculation to subroutine ridge_shift
        add flux_bio calculation to subroutine ridge_ice
    icepack_therm_itd.F90
        update calculation of dvssl and dvint in subroutine lateral_melt
    icepack_zbgc.F90
        lots of stuff
    icepack_zbgc_shared.F90
        lots of stuff

Remove redundant arguments in non-BGC interfaces.

    icepack_atmo.F90
    icepack_fsd.F90
    icepack_isotope.F90
    icepack_itd.F90
    icepack_meltpond_topo.F90
    icepack_mushy_physics.F90
    icepack_snow.F90
    icepack_therm_bl99.F90
    icepack_therm_mushy.F90
    icepack_therm_shared.F90
    icepack_therm_vertical.F90
    icepack_tracers.F90
    icepack_wavefracspec.F90

Generalize merge_fluxes to make all arguments optional

Fix bug in subroutine snow_redist computation of hsn_new when nslyr=1

Update the Icepack driver consistent with Icepack interface changes

Update zbgc initialization in the Icepack driver, move bgc parameter intiialization to icepack_init_parameters, update icepack_init_zbgc call.

Update some calls to icepack_warnings_setabort to add file and line number.

Update warning package to improve OpenMP robustness. Fixes potential race conditions that show up when lots of output is produced.

Modified congel implementation in subroutine thickness_changes in icepack_therm_vertical.F90 to recover bit-for-bit results. New congel and new bgc implementation were bit-for-bit independently, but when combining the changes, the intel compiler (with -O2) introduces non bit-for-bit changes (roundoff).

Update bgc namelist defaults and settings in icepack_in

Update bgc tracer sizing in set_env files

Update testing, remove skl tests, add zaero tests.

---------

Co-authored-by: Elizabeth Hunke <[email protected]>
Co-authored-by: David Bailey <[email protected]>
Co-authored-by: Nicole Jeffery <[email protected]>

columnphysics/icepack_atmo.F90

@@ -21,7 +21,7 @@ module icepack_atmo
       use icepack_parameters, only: pih, dragio, rhoi, rhos, rhow
       use icepack_parameters, only: atmbndy, calc_strair, formdrag
       use icepack_parameters, only: icepack_chkoptargflag
-      use icepack_tracers, only: n_iso
+      use icepack_tracers, only: ncat, n_iso
       use icepack_tracers, only: tr_iso
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
@@ -521,13 +521,10 @@ subroutine neutral_drag_coeffs (apnd,     hpnd,     &
                                       hdraft,   hridge,          &
                                       distrdg,  hkeel,           &
                                       dkeel,    lfloe,           &
-                                      dfloe,    ncat)
+                                      dfloe)
 
       use icepack_tracers, only: tr_pond
 
-      integer (kind=int_kind), intent(in) :: &
-         ncat
-
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          apnd     ,& ! melt pond fraction of sea ice
          hpnd     ,& ! mean melt pond depth over sea ice

columnphysics/icepack_flux.F90

@@ -73,7 +73,9 @@ subroutine merge_fluxes (aicen,                &
 
       ! single category fluxes
       real (kind=dbl_kind), intent(in) :: &
-          aicen   , & ! concentration of ice
+          aicen       ! concentration of ice
+
+      real (kind=dbl_kind), optional, intent(in) :: &
           flw     , & ! downward longwave flux          (W/m**2)
           strairxn, & ! air/ice zonal  strss,           (N/m**2)
           strairyn, & ! air/ice merdnl strss,           (N/m**2)
@@ -115,7 +117,7 @@ subroutine merge_fluxes (aicen,                &
           Urefn       ! air speed reference level       (m/s)
 
       ! cumulative fluxes
-      real (kind=dbl_kind), intent(inout) :: &
+      real (kind=dbl_kind), optional, intent(inout) :: &
           strairxT, & ! air/ice zonal  strss,           (N/m**2)
           strairyT, & ! air/ice merdnl strss,           (N/m**2)
           Cdn_atm_ratio, & ! ratio of total drag over neutral drag
@@ -148,25 +150,23 @@ subroutine merge_fluxes (aicen,                &
           ilpnd     ! pond loss/gain (+/-) to ice lid freezing/melting (m/step)
 
       real (kind=dbl_kind), intent(inout), optional :: &
-          fswthru_vdr , & ! vis dir sw radiation through ice bot    (W/m**2)
-          fswthru_vdf , & ! vis dif sw radiation through ice bot    (W/m**2)
-          fswthru_idr , & ! nir dir sw radiation through ice bot    (W/m**2)
-          fswthru_idf     ! nir dif sw radiation through ice bot    (W/m**2)
-
-      real (kind=dbl_kind), intent(inout), optional :: &
+          fswthru_vdr, & ! vis dir sw radiation through ice bot    (W/m**2)
+          fswthru_vdf, & ! vis dif sw radiation through ice bot    (W/m**2)
+          fswthru_idr, & ! nir dir sw radiation through ice bot    (W/m**2)
+          fswthru_idf, & ! nir dif sw radiation through ice bot    (W/m**2)
           dsnow,    & ! change in snow depth            (m)
           Uref        ! air speed reference level       (m/s)
 
-      real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
-          Qref_iso, & ! isotope air sp hum ref level    (kg/kg)
-          fiso_ocn, & ! isotope fluxes to ocean         (kg/m2/s)
-          fiso_evap   ! isotope evaporation             (kg/m2/s)
-
       real (kind=dbl_kind), dimension(:), intent(in), optional :: &
           Qrefn_iso, & ! isotope air sp hum ref level   (kg/kg)
           fiso_ocnn, & ! isotope fluxes to ocean        (kg/m2/s)
           fiso_evapn   ! isotope evaporation            (kg/m2/s)
 
+      real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
+          Qref_iso, & ! isotope air sp hum ref level    (kg/kg)
+          fiso_ocn, & ! isotope fluxes to ocean         (kg/m2/s)
+          fiso_evap   ! isotope evaporation             (kg/m2/s)
+
       character(len=*),parameter :: subname='(merge_fluxes)'
 
       !-----------------------------------------------------------------
@@ -178,23 +178,38 @@ subroutine merge_fluxes (aicen,                &
 
       ! atmo fluxes
 
-      strairxT   = strairxT + strairxn  * aicen
-      strairyT   = strairyT + strairyn  * aicen
-      Cdn_atm_ratio = Cdn_atm_ratio + &
-                      Cdn_atm_ratio_n   * aicen
-      fsurf      = fsurf    + fsurfn    * aicen
-      fcondtop   = fcondtop + fcondtopn * aicen
-      fcondbot   = fcondbot + fcondbotn * aicen
-      fsens      = fsens    + fsensn    * aicen
-      flat       = flat     + flatn     * aicen
-      fswabs     = fswabs   + fswabsn   * aicen
-      flwout     = flwout   &
-           + (flwoutn - (c1-emissivity)*flw) * aicen
-      evap       = evap     + evapn     * aicen
-      evaps      = evaps    + evapsn    * aicen
-      evapi      = evapi    + evapin    * aicen
-      Tref       = Tref     + Trefn     * aicen
-      Qref       = Qref     + Qrefn     * aicen
+      if (present(strairxn) .and. present(strairxT)) &
+         strairxT   = strairxT + strairxn  * aicen
+      if (present(strairyn) .and. present(strairyT)) &
+         strairyT   = strairyT + strairyn  * aicen
+      if (present(Cdn_atm_ratio_n) .and. present(Cdn_atm_ratio)) &
+         Cdn_atm_ratio = Cdn_atm_ratio + &
+                         Cdn_atm_ratio_n   * aicen
+      if (present(fsurfn) .and. present(fsurf)) &
+         fsurf      = fsurf    + fsurfn    * aicen
+      if (present(fcondtopn) .and. present(fcondtop)) &
+         fcondtop   = fcondtop + fcondtopn * aicen
+      if (present(fcondbotn) .and. present(fcondbot)) &
+         fcondbot   = fcondbot + fcondbotn * aicen
+      if (present(fsensn) .and. present(fsens)) &
+         fsens      = fsens    + fsensn    * aicen
+      if (present(flatn) .and. present(flat)) &
+         flat       = flat     + flatn     * aicen
+      if (present(fswabsn) .and. present(fswabs)) &
+         fswabs     = fswabs   + fswabsn   * aicen
+      if (present(flwoutn) .and. present(flwout) .and. present(flw)) &
+         flwout     = flwout   &
+              + (flwoutn - (c1-emissivity)*flw) * aicen
+      if (present(evapn) .and. present(evap)) &
+         evap       = evap     + evapn     * aicen
+      if (present(evapsn) .and. present(evaps)) &
+         evaps      = evaps    + evapsn    * aicen
+      if (present(evapin) .and. present(evapi)) &
+         evapi      = evapi    + evapin    * aicen
+      if (present(Trefn) .and. present(Tref)) &
+         Tref       = Tref     + Trefn     * aicen
+      if (present(Qrefn) .and. present(Qref)) &
+         Qref       = Qref     + Qrefn     * aicen
 
       ! Isotopes
       if (tr_iso) then
@@ -211,35 +226,46 @@ subroutine merge_fluxes (aicen,                &
 
       ! ocean fluxes
       if (present(Urefn) .and. present(Uref)) then
-         Uref = Uref     + Urefn     * aicen
+         Uref      = Uref      + Urefn     * aicen
       endif
 
-      fresh     = fresh     + freshn    * aicen
-      fsalt     = fsalt     + fsaltn    * aicen
-      fhocn     = fhocn     + fhocnn    * aicen
-      fswthru   = fswthru   + fswthrun  * aicen
-      if (present(fswthru_vdr)) &
-         fswthru_vdr   = fswthru_vdr   + fswthrun_vdr  * aicen
-      if (present(fswthru_vdf)) &
-         fswthru_vdf   = fswthru_vdf   + fswthrun_vdf  * aicen
-      if (present(fswthru_idr)) &
-         fswthru_idr   = fswthru_idr   + fswthrun_idr  * aicen
-      if (present(fswthru_idf)) &
-         fswthru_idf   = fswthru_idf   + fswthrun_idf  * aicen
+      if (present(freshn) .and. present(fresh)) &
+         fresh     = fresh     + freshn    * aicen
+      if (present(fsaltn) .and. present(fsalt)) &
+         fsalt     = fsalt     + fsaltn    * aicen
+      if (present(fhocnn) .and. present(fhocn)) &
+         fhocn     = fhocn     + fhocnn    * aicen
+      if (present(fswthrun) .and. present(fswthru)) &
+         fswthru   = fswthru   + fswthrun  * aicen
+
+      if (present(fswthrun_vdr) .and. present(fswthru_vdr)) &
+         fswthru_vdr = fswthru_vdr + fswthrun_vdr  * aicen
+      if (present(fswthrun_vdf) .and. present(fswthru_vdf)) &
+         fswthru_vdf = fswthru_vdf + fswthrun_vdf  * aicen
+      if (present(fswthrun_idr) .and. present(fswthru_idr)) &
+         fswthru_idr = fswthru_idr + fswthrun_idr  * aicen
+      if (present(fswthrun_idf) .and. present(fswthru_idf)) &
+         fswthru_idf = fswthru_idf + fswthrun_idf  * aicen
 
       ! ice/snow thickness
 
-      meltt     = meltt     + melttn    * aicen
-      meltb     = meltb     + meltbn    * aicen
-      melts     = melts     + meltsn    * aicen
+      if (present(melttn) .and. present(meltt)) &
+         meltt     = meltt     + melttn    * aicen
+      if (present(meltbn) .and. present(meltb)) &
+         meltb     = meltb     + meltbn    * aicen
+      if (present(meltsn) .and. present(melts)) &
+         melts     = melts     + meltsn    * aicen
       if (snwgrain) then
-         meltsliq  = meltsliq  + meltsliqn * aicen
+         if (present(meltsliqn) .and. present(meltsliq)) &
+            meltsliq  = meltsliq  + meltsliqn * aicen
       endif
-      if (present(dsnow)) then
+      if (present(dsnown) .and. present(dsnow)) then
          dsnow     = dsnow     + dsnown    * aicen
       endif
-      congel    = congel    + congeln   * aicen
-      snoice    = snoice    + snoicen   * aicen
+      if (present(congeln) .and. present(congel)) &
+         congel    = congel    + congeln   * aicen
+      if (present(snoicen) .and. present(snoice)) &
+         snoice    = snoice    + snoicen   * aicen
 
       ! Meltwater fluxes
       if (tr_pond) then