Add GEOS coupling capability including heat and mass flux options in Icepack (CICE-Consortium#1012)

Add GEOS coupling updates. This allows coupling to the GEOS coupled system where a semi-implicit thermodynamic coupling scheme is introduced. Similar to the explicit case, the fields fsurfn are provided by the coupler, along with their derivatives with respect to surface temperature dfsurfn_dTs. In this case, calc_Tsfc is still set to true, allowing ice surface and internal temperature to be updated implicitly. The resultant surface temperature change is passed back to the atmosphere model via coupler to complete the full update of its temperature profiles. This middle-ground approach, enabled by semi_implicit_Tsfc=true, does not sacrifice accuracy because it does not need limiting effective conductivity as in the explicit case. In addition, in GEOS, the atmosphere model assumes vapor deposits or sublimates at 0 degC. In this case, mass conservation is enforced and the resulting discrepancy in energy is resolved by another term, de_vapor, and passed to ocean. This option is only on when vapor_flux_correction is true.

Add 4 new shortwave terms, uvrdr, uvrdf, pardr, pardf to the coupling. These terms represent a breakdown of the direct and diffuse visible shortwave terms into two components, par = photosynthetical active radiation (400-700nm) and uvr = rest of the visible shortwave term (>700nm). The current visible shortwave is exactly represented by these two components. Includes adding atm forcing and terms associated with radiation passthru to the ocean.

Add support for GEOS semi-implicit coupling of surface temperature. In GEOS, surface and latent heat flux is computed in the atmosphere at 0degC. The sea ice model has to respect that calculation, but then computes the d(dh)/dTs terms to correct the heatflux for the sea ice temperature which is then applied conservatively in the coupled system. Implementation includes turning off some of the heat flux calculations in Icepack. This is controlled by the semi_implicit_Tsfc namelist.

Add calculation of a vapor flux correction. A correction is needed for GEOS coupling to compute a mass and enthalpy correction for evaporation and sublimation. This is controlled by the vapor_flux_correction namelist.

Add mapl/geos coupling directory and coupling files

Add opmask (orphan mask) for points that are NOT ocean/ice in the ocean/ice model but are ocean/ice in the atmosphere model. This allows for thermodyanamic calculations on the orphan points while not being involved in any sea ice dynamics. opmask determined by ocn_gridcell_frac which should be set by coupling layer at initialization.

Add geosmom grid_type to read GEOS MOM grid files

Add discover port

Update documentation

Change use of grid_type='tripole' where ns_boundary_type='tripole' is more appropriate.

Add distribution_wght=blockfull option to move away from CPP in init_domain_blocks. This option turns off land block elimination and sets all blocks to maximum weight for distribution.

Clean up some declarations in ice_flux.F90 and ice_state.F90

Add some new initialization output for new GEOS options and to clean up grid_type output


---------

Co-authored-by: bzhao <[email protected]>

Author: apcraig

cicecore/cicedyn/dynamics/ice_dyn_eap.F90

@@ -2090,8 +2090,7 @@ subroutine read_restart_eap()
       use ice_boundary, only: ice_HaloUpdate_stress
       use ice_constants, only:  &
           field_loc_center, field_type_scalar
-      use ice_domain, only: nblocks, halo_info
-      use ice_grid, only: grid_type
+      use ice_domain, only: nblocks, halo_info, ns_boundary_type
       use ice_restart, only: read_restart_field
 
       ! local variables
@@ -2131,7 +2130,8 @@ subroutine read_restart_eap()
       call read_restart_field(nu_restart_eap,0,a12_4,'ruf8', &
            'a12_4',1,diag,field_loc_center,field_type_scalar) ! a12_4
 
-      if (trim(grid_type) == 'tripole') then
+      if (trim(ns_boundary_type) == 'tripole' .or. &
+          trim(ns_boundary_type) == 'tripoleT') then
 
          call ice_HaloUpdate_stress(a11_1, a11_3,      halo_info, &
                                     field_loc_center,  field_type_scalar)

cicecore/cicedyn/dynamics/ice_dyn_evp.F90

@@ -259,7 +259,8 @@ subroutine evp (dt)
       use ice_boundary, only: ice_halo, ice_HaloMask, ice_HaloUpdate, &
           ice_HaloDestroy, ice_HaloUpdate_stress
       use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
-      use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn
+      use ice_domain, only: nblocks, blocks_ice, halo_info, maskhalo_dyn, &
+          ns_boundary_type
       use ice_domain_size, only: max_blocks, ncat
       use ice_flux, only: rdg_conv, rdg_shear, strairxT, strairyT, &
           strairxU, strairyU, uocn, vocn, ss_tltx, ss_tlty, fmU, &
@@ -280,8 +281,7 @@ subroutine evp (dt)
       use ice_grid, only: tmask, umask, umaskCD, nmask, emask, uvm, epm, npm, &
           dxE, dxN, dxT, dxU, dyE, dyN, dyT, dyU, &
           tarear, uarear, earear, narear, grid_average_X2Y, uarea, &
-          grid_type, grid_ice, &
-          grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
+          grid_ice, grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
       use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, uvelN, vvelN, &
           uvelE, vvelE, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
@@ -1313,7 +1313,8 @@ subroutine evp (dt)
       endif
 
       ! Force symmetry across the tripole seam
-      if (trim(grid_type) == 'tripole') then
+      if (trim(ns_boundary_type) == 'tripole' .or. &
+          trim(ns_boundary_type) == 'tripoleT') then
          ! TODO: C/CD-grid
       if (maskhalo_dyn) then
          !-------------------------------------------------------

cicecore/cicedyn/dynamics/ice_dyn_shared.F90

@@ -240,7 +240,8 @@ subroutine init_dyn_shared (dt)
 
       use ice_blocks, only: block, get_block
       use ice_boundary, only: ice_halo, ice_haloUpdate
-      use ice_domain, only: nblocks, halo_dynbundle, blocks_ice, halo_info
+      use ice_domain, only: nblocks, halo_dynbundle, blocks_ice, halo_info, &
+          ns_boundary_type
       use ice_domain_size, only: max_blocks
       use ice_flux, only: &
           stressp_1, stressp_2, stressp_3, stressp_4, &
@@ -268,6 +269,13 @@ subroutine init_dyn_shared (dt)
 
       character(len=*), parameter :: subname = '(init_dyn_shared)'
 
+      ! checks
+      if (kdyn == 1 .and. evp_algorithm == 'shared_mem_1d' .and. &
+          (ns_boundary_type == 'tripole' .or. ns_boundary_type == 'tripoleT')) then
+         call abort_ice(subname//' ERROR: evp_alg shared mem 1d not supported with tripole', &
+            file=__FILE__, line=__LINE__)
+      endif
+
       call set_evp_parameters (dt)
       ! allocate dyn shared (init_uvel,init_vvel)
       call alloc_dyn_shared

cicecore/cicedyn/dynamics/ice_dyn_vp.F90

@@ -155,7 +155,8 @@ subroutine implicit_solver (dt)
       use ice_boundary, only: ice_HaloMask, ice_HaloUpdate, &
           ice_HaloDestroy, ice_HaloUpdate_stress
       use ice_blocks, only: block, get_block, nx_block, ny_block
-      use ice_domain, only: blocks_ice, halo_info, maskhalo_dyn
+      use ice_domain, only: blocks_ice, halo_info, maskhalo_dyn, &
+          ns_boundary_type
       use ice_domain_size, only: max_blocks, ncat
       use ice_dyn_shared, only: deformations, iceTmask, iceUmask, &
           cxp, cyp, cxm, cym
@@ -168,7 +169,7 @@ subroutine implicit_solver (dt)
           stressm_1, stressm_2, stressm_3, stressm_4, &
           stress12_1, stress12_2, stress12_3, stress12_4
       use ice_grid, only: tmask, umask, dxT, dyT, dxU, dyU, &
-          tarear, grid_type, grid_average_X2Y, &
+          tarear, grid_average_X2Y, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
       use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
@@ -546,7 +547,8 @@ subroutine implicit_solver (dt)
       endif
 
       ! Force symmetry across the tripole seam
-      if (trim(grid_type) == 'tripole') then
+      if (trim(ns_boundary_type) == 'tripole' .or. &
+          trim(ns_boundary_type) == 'tripoleT') then
       if (maskhalo_dyn) then
          !-------------------------------------------------------
          ! set halomask to zero because ice_HaloMask always keeps

cicecore/cicedyn/general/ice_flux.F90

@@ -114,8 +114,7 @@ module ice_flux
          dvirdgdt, & ! rate of ice volume ridged (m/s)
          opening     ! rate of opening due to divergence/shear (1/s)
 
-      real (kind=dbl_kind), &
-         dimension (:,:,:,:), allocatable, public :: &
+      real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
        ! ridging diagnostics in categories
          dardg1ndt, & ! rate of area loss by ridging ice (1/s)
          dardg2ndt, & ! rate of area gain by new ridges (1/s)
@@ -177,13 +176,26 @@ module ice_flux
        ! NOTE: when in CICE_IN_NEMO mode, these are gridbox mean fields,
        ! not per ice area. When in standalone mode, these are per ice area.
 
-      real (kind=dbl_kind), &
-         dimension (:,:,:,:), allocatable, public :: &
+      real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
          fsurfn_f   , & ! net flux to top surface, excluding fcondtop
          fcondtopn_f, & ! downward cond flux at top surface (W m-2)
          fsensn_f   , & ! sensible heat flux (W m-2)
          flatn_f        ! latent heat flux (W m-2)
 
+      ! in from atmosphere
+      ! required for coupling in GEOS
+
+      real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
+         evapn_f,      & ! evaporation/sublimation (kg m-2 s-1)
+         dflatndTsfc_f,  & ! derivative of latent flux w.r.t. Tsfc
+         dfsurfndTsfc_f    ! derivative of surface flux w.r.t. Tsfc
+
+      real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
+         swuvrdr     , & !  vis uvr flux, direct (W m-2)
+         swuvrdf     , & !  vis uvr flux, diffuse (W m-2)
+         swpardr     , & !  vis par flux, direct (W m-2)
+         swpardf         !  vis par flux, diffuse (W m-2)
+
        ! in from atmosphere
 
       real (kind=dbl_kind), dimension (:,:,:), allocatable, public :: &
@@ -258,7 +270,11 @@ module ice_flux
          fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
          fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2)
          fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2)
-         fswthru_idf     ! nir dif shortwave penetrating to ocean (W/m^2)
+         fswthru_idf , & ! nir dif shortwave penetrating to ocean (W/m^2)
+         fswthru_uvrdr,& ! vis dir uvr SW penetrating to ocean (W/m^2)
+         fswthru_uvrdf,& ! vis dif uvr SW penetrating to ocean (W/m^2)
+         fswthru_pardr,& ! nir dir par SW penetrating to ocean (W/m^2)
+         fswthru_pardf   ! nir dif par SW penetrating to ocean (W/m^2)
 
        ! internal
 
@@ -326,16 +342,14 @@ module ice_flux
          frz_onset, &! day of year that freezing begins (congel or frazil)
          frazil_diag ! frazil ice growth diagnostic (m/step-->cm/day)
 
-      real (kind=dbl_kind), &
-         dimension (:,:,:,:), allocatable, public :: &
+      real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
          fsurfn,   & ! category fsurf
          fcondtopn,& ! category fcondtop
          fcondbotn,& ! category fcondbot
          fsensn,   & ! category sensible heat flux
          flatn       ! category latent heat flux
 
-      real (kind=dbl_kind), &
-         dimension (:,:,:,:), allocatable, public :: &
+      real (kind=dbl_kind), dimension (:,:,:,:), allocatable, public :: &
          snwcnt       ! counter for presence of snow
 
       ! As above but these remain grid box mean values i.e. they are not
@@ -453,6 +467,10 @@ subroutine alloc_flux
          swvdf      (nx_block,ny_block,max_blocks), & ! sw down, visible, diffuse (W/m^2)
          swidr      (nx_block,ny_block,max_blocks), & ! sw down, near IR, direct  (W/m^2)
          swidf      (nx_block,ny_block,max_blocks), & ! sw down, near IR, diffuse (W/m^2)
+         swuvrdr    (nx_block,ny_block,max_blocks), & ! vis uvr flux, direct (W m-2)
+         swuvrdf    (nx_block,ny_block,max_blocks), & ! vis uvr flux, diffuse (W m-2)
+         swpardr    (nx_block,ny_block,max_blocks), & ! vis par flux, direct (W m-2)
+         swpardf    (nx_block,ny_block,max_blocks), & ! vis par flux, diffuse (W m-2)
          flw        (nx_block,ny_block,max_blocks), & ! incoming longwave radiation (W/m^2)
          frain      (nx_block,ny_block,max_blocks), & ! rainfall rate (kg/m^2 s)
          fsnow      (nx_block,ny_block,max_blocks), & ! snowfall rate (kg/m^2 s)
@@ -499,11 +517,15 @@ subroutine alloc_flux
          fhocn      (nx_block,ny_block,max_blocks), & ! net heat flux to ocean (W/m^2)
          fsloss     (nx_block,ny_block,max_blocks), & ! rate of snow loss to leads (kg/m^2/s)
          fswthru    (nx_block,ny_block,max_blocks), & ! shortwave penetrating to ocean (W/m^2)
-         fswthru_vdr (nx_block,ny_block,max_blocks), & ! vis dir shortwave penetrating to ocean (W/m^2)
-         fswthru_vdf (nx_block,ny_block,max_blocks), & ! vis dif shortwave penetrating to ocean (W/m^2)
-         fswthru_idr (nx_block,ny_block,max_blocks), & ! nir dir shortwave penetrating to ocean (W/m^2)
-         fswthru_idf (nx_block,ny_block,max_blocks), & ! nir dif shortwave penetrating to ocean (W/m^2)
-         scale_factor (nx_block,ny_block,max_blocks), & ! scaling factor for shortwave components
+         fswthru_vdr(nx_block,ny_block,max_blocks), & ! vis dir shortwave penetrating to ocean (W/m^2)
+         fswthru_vdf(nx_block,ny_block,max_blocks), & ! vis dif shortwave penetrating to ocean (W/m^2)
+         fswthru_idr(nx_block,ny_block,max_blocks), & ! nir dir shortwave penetrating to ocean (W/m^2)
+         fswthru_idf(nx_block,ny_block,max_blocks), & ! nir dif shortwave penetrating to ocean (W/m^2)
+         fswthru_uvrdr (nx_block,ny_block,max_blocks), & ! vis dir uvr SW penetrating to ocean (W/m^2)
+         fswthru_uvrdf (nx_block,ny_block,max_blocks), & ! vis dir uvr SW penetrating to ocean (W/m^2)
+         fswthru_pardr (nx_block,ny_block,max_blocks), & ! vis dir par SW penetrating to ocean (W/m^2)
+         fswthru_pardf (nx_block,ny_block,max_blocks), & ! vis dir par SW penetrating to ocean (W/m^2)
+         scale_factor  (nx_block,ny_block,max_blocks), & ! scaling factor for shortwave components
          strairx_ocn(nx_block,ny_block,max_blocks), & ! stress on ocean by air, x-direction
          strairy_ocn(nx_block,ny_block,max_blocks), & ! stress on ocean by air, y-direction
          fsens_ocn  (nx_block,ny_block,max_blocks), & ! sensible heat flux (W/m^2)
@@ -566,6 +588,9 @@ subroutine alloc_flux
          fcondtopn_f(nx_block,ny_block,ncat,max_blocks), & ! downward cond flux at top surface (W m-2)
          fsensn_f   (nx_block,ny_block,ncat,max_blocks), & ! sensible heat flux (W m-2)
          flatn_f    (nx_block,ny_block,ncat,max_blocks), & ! latent heat flux (W m-2)
+         evapn_f    (nx_block,ny_block,ncat,max_blocks), & ! evaporative water flux (kg/m^2/s) by atmosphere model
+         dflatndTsfc_f (nx_block,ny_block,ncat,max_blocks), & ! derivative of flatn with respect to Tsfc
+         dfsurfndTsfc_f(nx_block,ny_block,ncat,max_blocks), & ! derivative of fsurfn with respect to Tsfc
          meltsn     (nx_block,ny_block,ncat,max_blocks), & ! snow melt in category n (m)
          melttn     (nx_block,ny_block,ncat,max_blocks), & ! top melt in category n (m)
          meltbn     (nx_block,ny_block,ncat,max_blocks), & ! bottom melt in category n (m)
@@ -584,6 +609,11 @@ subroutine alloc_flux
          stat=ierr)
       if (ierr/=0) call abort_ice('(alloc_flux): Out of memory')
 
+      swuvrdr(:,:,:) = c0
+      swuvrdf(:,:,:) = c0
+      swpardr(:,:,:) = c0
+      swpardf(:,:,:) = c0
+
       if (grid_ice == "CD" .or. grid_ice == "C") &
          allocate( &
          taubxN     (nx_block,ny_block,max_blocks), & ! seabed stress (x) at N points (N/m^2)
@@ -792,6 +822,10 @@ subroutine init_coupler_flux
       fswthru_vdf (:,:,:) = c0
       fswthru_idr (:,:,:) = c0
       fswthru_idf (:,:,:) = c0
+      fswthru_uvrdr (:,:,:) = c0
+      fswthru_uvrdf (:,:,:) = c0
+      fswthru_pardr (:,:,:) = c0
+      fswthru_pardf (:,:,:) = c0
       fresh_da(:,:,:) = c0    ! data assimilation
       fsalt_da(:,:,:) = c0
       flux_bio (:,:,:,:) = c0 ! bgc
@@ -853,6 +887,8 @@ subroutine init_flux_atm
       ! strairxT(:,:,:) = 0.15_dbl_kind
       ! strairyT(:,:,:) = 0.15_dbl_kind
 
+      fsurf   (:,:,:) = c0
+      fcondtop(:,:,:) = c0
       fsens   (:,:,:) = c0
       flat    (:,:,:) = c0
       fswabs  (:,:,:) = c0
@@ -898,6 +934,10 @@ subroutine init_flux_ocn
       fswthru_vdf  (:,:,:)   = c0
       fswthru_idr  (:,:,:)   = c0
       fswthru_idf  (:,:,:)   = c0
+      fswthru_uvrdr(:,:,:)   = c0
+      fswthru_uvrdf(:,:,:)   = c0
+      fswthru_pardr(:,:,:)   = c0
+      fswthru_pardf(:,:,:)   = c0
 
       faero_ocn (:,:,:,:) = c0
       fiso_ocn  (:,:,:,:) = c0