AWI-ESM3 CO2 coupling

config/namelist.io.recom.2p3z2d

@@ -0,0 +1,95 @@
+&diag_list
+ldiag_solver     =.false.
+lcurt_stress_surf=.false.
+ldiag_curl_vel3  =.false.
+ldiag_Ri         =.false.
+ldiag_turbflux   =.false.
+ldiag_salt3D     =.false.
+ldiag_dMOC       =.false.
+ldiag_DVD        =.false.
+ldiag_forc       =.false.
+ldiag_extflds    =.false.
+ldiag_destine    =.false.
+ldiag_trflx      =.false.
+ldiag_uvw_sqr    =.false.
+ldiag_trgrd_xyz  =.false.
+/
+
+&nml_general
+io_listsize    =150 !number of streams to allocate. shallbe large or equal to the number of streams in &nml_list
+vec_autorotate =.false.
+compression_level = 1
+/
+
+! for sea ice related variables use_ice should be true, otherewise there will be no output
+! for 'curl_surf' to work lcurt_stress_surf must be .true. otherwise no output
+! for 'fer_C', 'bolus_u', 'bolus_v', 'bolus_w', 'fer_K' to work Fer_GM must be .true. otherwise no output
+! 'otracers' - all other tracers if applicable
+! for 'dMOC' to work ldiag_dMOC must be .true. otherwise no output
+&nml_list
+io_list =  'sst       ',1, 'm', 4,
+           'sss       ',1, 'm', 4,
+    	   'ssh       ',1, 'm', 4,
+           'uice      ',1, 'm', 4,
+           'vice      ',1, 'm', 4,
+           'a_ice     ',1, 'm', 4,
+           'm_ice     ',1, 'm', 4,
+           'm_snow    ',1, 'm', 4,
+           'MLD1      ',1, 'm', 4,
+           'MLD2      ',1, 'm', 4,
+           'MLD3      ',1, 'm', 4,
+           'tx_sur    ',1, 'm', 4,
+           'ty_sur    ',1, 'm', 4,
+           'prec      ',1, 'm', 4,
+           'snow      ',1, 'm', 4,
+           'evap      ',1, 'm', 4,
+           'runoff    ',1, 'm', 4,
+           'temp      ',1, 'm', 4,
+           'salt      ',1, 'm', 8,
+           'otracers  ',1, 'm', 4,
+           'N2        ',1, 'y', 4,
+           'Kv        ',1, 'y', 4,
+           'u         ',1, 'm', 4,
+           'v         ',1, 'm', 4,
+           'unod      ',1, 'm', 4,
+           'vnod      ',1, 'm', 4,
+           'w         ',1, 'm', 4,
+           'Av        ',1, 'y', 4,
+           'bolus_u   ',1, 'y', 4,
+           'bolus_v   ',1, 'y', 4,
+           'bolus_w   ',1, 'y', 4,
+           'fw        ',1, 'm', 4,
+           'fh        ',1, 'm', 4,
+           'dpCO2s    ',1, 'm', 4,
+           'pCO2s     ',1, 'm', 4,
+           'CO2f      ',1, 'm', 4,
+           'xCO2atm   ',4, 's', 4,
+           'Hp        ',1, 'm', 4,
+           'aFe       ',1, 'm', 4,
+           'aN        ',1, 'm', 4,
+           'denb      ',1, 'm', 4,
+           'benN      ',1, 'm', 4,
+           'benC      ',1, 'm', 4,
+           'benSi     ',1, 'm', 4,
+           'benCalc   ',1, 'm', 4,
+           'Chldegd   ',1, 'm', 4,
+           'Chldegn   ',1, 'm', 4,
+           'NNAd      ',1, 'm', 4,
+           'NNAn      ',1, 'm', 4,
+           'GPPd      ',1, 'm', 4,
+           'GPPn      ',1, 'm', 4,
+           'NPPd      ',1, 'm', 4,
+           'NPPn      ',1, 'm', 4,
+           'PAR       ',1, 'm', 4,
+           'respmeso',1, 'm', 4,
+           'respmacro',1, 'm', 4,
+           'respmicro',1, 'm', 4,
+           'calcdiss',1, 'm', 4,
+           'calcif',1, 'm', 4,
+           'aggn',1, 'm', 4,
+           'aggd',1, 'm', 4,
+           'docexn',1, 'm', 4,
+           'docexd',1, 'm', 4,
+           'respn',1, 'm', 4,
+           'respd',1, 'm', 4,
+/

config/namelist.tra.recom.2p3z2d

@@ -0,0 +1,90 @@
+&tracer_listsize
+num_tracers=100 !number of tracers to allocate. shallbe large or equal to the number of streams in &nml_list
+/
+
+&tracer_list
+nml_tracer_list =  
+1  ,   'MFCT',       'QR4C',       'FCT ',        1.,          1., 
+2  ,   'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1001,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1002,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1003,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1004,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1005,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1006,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1007,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1008,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1009,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1010,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1011,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1012,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1013,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1014,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1015,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1016,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1017,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1018,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1019,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1020,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1021,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1022,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1023,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1024,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1025,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1026,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1027,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1028,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1029,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+1030,  'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+!1031, 'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+!1032, 'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+!1033, 'MFCT',       'QR4C',       'FCT ',        1.,          1.,
+!101, 'UPW1',       'UPW1',       'NON ',        0.,          0.
+/
+
+&tracer_init3d                            ! initial conditions for tracers
+n_ic3d   = 8                              ! number of tracers to initialize
+idlist   = 1019, 1022, 1018, 1003, 1002, 1001, 2, 1 ! their IDs (0 is temperature, 1 is salinity, etc.). The reading order is defined here!
+filelist = 'fe_pisces_opa_eq_init_3D_changed_name.nc', 'woa18_all_o00_01_mmol_fesom2.nc', 'woa13_all_i00_01_fesom2.nc', 'GLODAPv2.2016b.TAlk_fesom2_mmol_fix_z_Fillvalue.nc', 'GLODAPv2.2016b.TCO2_fesom2_mmol_fix_z_Fillvalue.nc', 'woa13_all_n00_01_fesom2.nc', 'phc3.0_winter.nc', 'phc3.0_winter.nc'  ! list of files in ClimateDataPath to read (one file per tracer), same order as idlist
+varlist  = 'Fe', 'oxygen_mmol', 'i_an', 'TAlk_mmol', 'TCO2_mmol', 'n_an', 'salt', 'temp'                 ! variables to read from specified files
+t_insitu = .true.                         ! if T is insitu it will be converted to potential after reading it
+/
+
+&tracer_init2d                                      ! initial conditions for 2D tracers (sea ice)
+n_ic2d   = 3                                        ! number of tracers to initialize
+idlist   = 1, 2, 3                                  ! their IDs (0 is a_ice, 1 is m_ice, 3 m_snow). The reading order is defined here!
+filelist = 'a_ice.nc', 'm_ice.nc', 'm_snow.nc' ! list of files in ClimateDataPath to read (one file per tracer), same order as idlist
+varlist  = 'a_ice', 'm_ice', 'm_snow'                 ! variables to read from specified files
+ini_ice_from_file=.false.
+/
+
+&tracer_general
+! bharmonic diffusion for tracers. We recommend to use this option in very high resolution runs (Redi is generally off there).
+smooth_bh_tra =.false.     ! use biharmonic diffusion (filter implementation) for tracers
+gamma0_tra    = 0.0005     ! gammaX_tra are analogous to those in the dynamical part
+gamma1_tra    = 0.0125
+gamma2_tra    = 0.
+i_vert_diff   =.true.
+/
+
+&tracer_phys
+use_momix     = .true.     ! switch on/off  !Monin-Obukhov -> TB04 mixing
+momix_lat     = -50.0      ! latitidinal treshhold for TB04, =90 --> global
+momix_kv      = 0.01       ! PP/KPP, mixing coefficient within MO length
+use_instabmix = .true.     ! enhance convection in case of instable stratification
+instabmix_kv  = 0.1
+use_windmix   = .false.    ! enhance mixing trough wind only for PP mixing (for stability)
+windmix_kv    = 1.e-3
+windmix_nl    = 2
+diff_sh_limit=5.0e-3       ! for KPP, max diff due to shear instability
+Kv0_const=.true.
+double_diffusion=.false.   ! for KPP,dd switch
+K_ver=1.0e-5
+K_hor=3000.
+surf_relax_T=0.0
+surf_relax_S=1.929e-06   ! 50m/300days 6.43e-07! m/s 10./(180.*86400.)
+balance_salt_water =.true. ! balance virtual-salt or freshwater flux or not
+clim_relax=0.0	           ! 1/s, geometrical information has to be supplied
+ref_sss_local=.true.
+ref_sss=34.
+/

src/cpl_driver.F90

@@ -17,6 +17,8 @@ module cpl_driver
   use o_param,  only : rad
   USE MOD_PARTIT
   use mpi
+
+
   implicit none
   save   
   !
@@ -28,8 +30,15 @@ module cpl_driver
   ! (final number of fields depends now on lwiso switch and is set in subroutine cpl_oasis3mct_define_unstr)
 
 #if defined (__oifs)
+#if defined(__recom)
+  integer                    :: nsend = 8
+  integer                    :: nrecv = 16
+!With oifs, without recom:
+#else
   integer                    :: nsend = 7
   integer                    :: nrecv = 15
+#endif
+!Without oifs
 #else
   integer                    :: nsend = 4
   integer                    :: nrecv = 12
@@ -661,6 +670,9 @@ subroutine cpl_oasis3mct_define_unstr(partit, mesh)
     cpl_send( 5)='sia_feom' ! 5. sea ice albedo [%-100]            ->
     cpl_send( 6)='u_feom'   ! 6. eastward  surface velocity [m/s]  ->
     cpl_send( 7)='v_feom'   ! 7. northward surface velocity [m/s]  ->
+#if defined (__recom)
+    cpl_send( 8)='FCO2_feom'! 8. CO2 flux [kgCO2 m-2 s-1]             ->
+#endif
 #else
     cpl_send( 1)='sst_feom' ! 1. sea surface temperature [°C]      ->
     cpl_send( 2)='sit_feom' ! 2. sea ice thickness [m]             ->
@@ -701,6 +713,10 @@ subroutine cpl_oasis3mct_define_unstr(partit, mesh)
     cpl_recv(13) = 'calv_oce'
     cpl_recv(14) = 'u10w_oce'
     cpl_recv(15) = 'v10w_oce'
+#if defined (__recom)
+    cpl_recv(16) = 'XCO2_oce'
+#endif
+!Not oifs
 #else
     cpl_recv(1)  = 'taux_oce'
     cpl_recv(2)  = 'tauy_oce'

src/fesom_module.F90

@@ -48,6 +48,7 @@ module fesom_main_storage_module
   use iceberg_step
   use iceberg_ocean_coupling
   use Toy_Channel_Soufflet, only: compute_zonal_mean
+  use ieee_arithmetic
   ! Define icepack module
 
 #if defined (__icepack)

src/gen_forcing_couple.F90

@@ -273,6 +273,9 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
 #endif
   use gen_bulk
   use force_flux_consv_interface
+#if defined (__recom)
+  use REcoM_GloVar, only: x_co2atm, GloCO2flux_seaicemask
+#endif
 
   implicit none
   integer,        intent(in)            :: istep
@@ -386,6 +389,13 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
               do n=1,myDim_nod2D+eDim_nod2D
                 exchange(n) = UVnode(2,1,n)
               end do
+#if defined (__recom)
+            elseif (i.eq.8) then
+              ! GloCO2flux_seaicemask is in [mmol/m2/s], need [kg/m2/s]
+              ! Conversion: 1.0e-3_WP -> mol/s -> kg/s
+              ! 1 mol CO2 = 44.0095 g/mol = 0.0440095 kg/mol (NIST 2018)
+              exchange(:) = GloCO2flux_seaicemask(:) * 1.0e-3_WP * 0.0440095_WP  ! [kg m⁻² s⁻¹]
+#endif
             else    
             print *, 'not installed yet or error in cpl_oasis3mct_send', mype
 #else
@@ -565,6 +575,15 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh)
              if (action) then
                 v_wind(:)                     = exchange(:)        ! meridional wind
              end if
+#if defined (__recom)
+         elseif (i.eq.16) then
+             if (action) then
+                ! Convert mass mixing ratio (kg/kg) to mole fraction (dimensionless)
+                ! MW_CO2 = 44.0095 g/mol (NIST 2018)
+                ! MW_dry_air = 28.9647 g/mol (standard atmosphere composition)
+                x_co2atm(:) = exchange(:) * ((28.9647_WP/44.0095_WP)*1e6_WP)  ! [mole fraction]
+             end if
+#endif
 #else
          elseif (i.eq.13) then
             if (action) then

src/gen_ic3d.F90

@@ -436,7 +436,9 @@ SUBROUTINE getcoeffld(tracers, partit, mesh)
          if (x<0.)   x=x+360.
          if (x>360.) x=x-360.
          if ( min(i,j)>0 ) then
-         if (any(ncdata(i:ip1,j:jp1,1) > dummy*0.99_WP)) cycle
+         ! CAVITY FIX: Check for NaN in climatology data which can occur near cavity regions
+         if (any(ncdata(i:ip1,j:jp1,1) > dummy*0.99_WP) .or. &
+             any(.not. ieee_is_finite(ncdata(i:ip1,j:jp1,1)))) cycle
             x1 = nc_lon(i)
             x2 = nc_lon(ip1)
             y1 = nc_lat(j)
@@ -447,7 +449,9 @@ SUBROUTINE getcoeffld(tracers, partit, mesh)
             data1d(:) = ( ncdata(i,j,:)   * (x2-x)*(y2-y)   + ncdata(ip1,j,:)     * (x-x1)*(y2-y) + &
                         ncdata(i,jp1,:) * (x2-x)*(y-y1)   + ncdata(ip1, jp1, :) * (x-x1)*(y-y1)     ) / denom
             where (ncdata(i,j,:)   > 0.99_WP*dummy .OR. ncdata(ip1,j,:)   > 0.99_WP*dummy .OR. &
-                    ncdata(i,jp1,:) > 0.99_WP*dummy .OR. ncdata(ip1,jp1,:) > 0.99_WP*dummy)
+                    ncdata(i,jp1,:) > 0.99_WP*dummy .OR. ncdata(ip1,jp1,:) > 0.99_WP*dummy .OR. &
+                    .not. ieee_is_finite(ncdata(i,j,:)) .OR. .not. ieee_is_finite(ncdata(ip1,j,:)) .OR. &
+                    .not. ieee_is_finite(ncdata(i,jp1,:)) .OR. .not. ieee_is_finite(ncdata(ip1,jp1,:)))
                 data1d(:)=dummy
             end where   
 
@@ -505,6 +509,14 @@ SUBROUTINE getcoeffld(tracers, partit, mesh)
                     end if
                 enddo
             end if ! --> if (use_cavity) then
+         else
+            ! CAVITY FIX: If bilinear interpolation fails (missing data), set fallback values
+            ! This prevents NaN tracers when climatology has dummy values near cavity regions
+            do k= ul1, nl1
+               if (.not. ieee_is_finite(tracers%data(current_tracer)%values(k,ii))) then
+                  tracers%data(current_tracer)%values(k,ii) = 0.0_WP
+               endif
+            enddo
          end if ! --> if ( min(i,j)>0 ) then
       end do !ii
       if (mype==0) then
@@ -546,7 +558,11 @@ SUBROUTINE do_ic3d(tracers, partit, mesh)
             ! get first coeficients for time inerpolation on model grid for all datas
             call getcoeffld(tracers, partit, mesh)
             call nc_end ! deallocate arrqays associated with netcdf file
-            call extrap_nod(tracers%data(current_tracer)%values(:,:), partit, mesh)
+            ! CAVITY FIX: Initialize to 0.0 before extrapolation to prevent NaN
+      where (.not. ieee_is_finite(tracers%data(current_tracer)%values(:,:)))
+         tracers%data(current_tracer)%values(:,:) = 0.0_WP
+      end where
+      call extrap_nod(tracers%data(current_tracer)%values(:,:), partit, mesh)
             exit
          elseif (current_tracer==tracers%num_tracers) then
             if (partit%mype==0) write(*,*) "idlist contains tracer which is not listed in tracer_id!"
@@ -560,6 +576,19 @@ SUBROUTINE do_ic3d(tracers, partit, mesh)
 
       do current_tracer=1, tracers%num_tracers
          !_________________________________________________________________________
+         ! CAVITY FIX: Clean up any remaining NaN values before dummy check
+         where (.not. ieee_is_finite(tracers%data(current_tracer)%values(:,:)))
+               tracers%data(current_tracer)%values(:,:) = 0.0_WP
+         end where
+         !_________________________________________________________________________
+         ! set remaining dummy values from bottom topography to 0.0_WP
+         where (tracers%data(current_tracer)%values > 0.9_WP*dummy)
+               tracers%data(current_tracer)%values=0.0_WP
+         end where
+
+         !_________________________________________________________________________
+         ! eliminate values within cavity that result from the extrapolation of 
+         ! initialisation
          ! set remaining dummy values and NaN from interpolation to 0.0_WP
          do n=1,partit%myDim_nod2d + partit%eDim_nod2D
             do i=1, mesh%nl-1