Add GEOS coupling capability including heat and mass flux options in Icepack

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