From 0a7b6cd5bee37f8b443fd9e22e1fcbecbfc336d2 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Mon, 21 Jul 2025 17:06:40 +0200 Subject: [PATCH 01/17] save wip towards co2 coupling with openifs48r1 --- src/cpl_driver.F90 | 16 ++++++++++++++++ src/gen_forcing_couple.F90 | 13 +++++++++++++ 2 files changed, 29 insertions(+) diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index a26b842f2..edc8bf824 100644 --- a/src/cpl_driver.F90 +++ b/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 @@ -658,6 +667,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] -> @@ -698,6 +710,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' diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 67ba51c43..0889b3de5 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -98,6 +98,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_locVar, only: x_co2atm, co2flux +#endif implicit none integer, intent(in) :: istep @@ -209,6 +212,10 @@ 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 + exchange(n) = co2flux(:) ! CO2 concentration +#endif else print *, 'not installed yet or error in cpl_oasis3mct_send', mype #else ! oifs @@ -384,6 +391,12 @@ 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 + x_co2atm(:) = exchange(:) ! Mass Mixing ratio CO2 + end if +#endif #else ! oifs elseif (i.eq.13) then if (action) then From b2f3de1f8547aabb49ad908ff474e70d82b61b17 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Tue, 22 Jul 2025 10:35:06 +0200 Subject: [PATCH 02/17] minor fixes to co2 coupling --- src/gen_forcing_couple.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 0889b3de5..f0804b9d4 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -99,7 +99,7 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) use gen_bulk use force_flux_consv_interface #if defined (__recom) - use REcoM_locVar, only: x_co2atm, co2flux + use REcoM_locVar, only: LocAtmCO2, co2flux #endif implicit none @@ -214,7 +214,7 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) end do #if defined (__recom) elseif (i.eq.8) then - exchange(n) = co2flux(:) ! CO2 concentration + exchange(:) = co2flux(:) ! CO2 concentration #endif else print *, 'not installed yet or error in cpl_oasis3mct_send', mype @@ -394,7 +394,7 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) #if defined (__recom) elseif (i.eq.16) then if (action) then - x_co2atm(:) = exchange(:) ! Mass Mixing ratio CO2 + LocAtmCO2(1) = exchange(1) ! Atmospheric CO2 [uatm] end if #endif #else ! oifs From 5e8a5255aedb897cf2bcfb4ab0278036fc3288fc Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Tue, 22 Jul 2025 11:51:28 +0200 Subject: [PATCH 03/17] maybe the right fields now? --- src/gen_forcing_couple.F90 | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index f0804b9d4..84b94ae60 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -99,7 +99,7 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) use gen_bulk use force_flux_consv_interface #if defined (__recom) - use REcoM_locVar, only: LocAtmCO2, co2flux + use recom_modules, only: x_co2atm, GloCO2flux #endif implicit none @@ -214,7 +214,12 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) end do #if defined (__recom) elseif (i.eq.8) then - exchange(:) = co2flux(:) ! CO2 concentration + ! Convert CO2 flux from mmol/(m² day) to kg/(m² s) + ! GloCO2flux is in [mmol/m2/day], need [kg/m2/s] + ! Conversion: mmol/day -> mol/s -> kg/s + ! 1 mmol/day = 1e-3 mol / (24*3600 s) = 1.157e-8 mol/s + ! 1 mol CO2 = 0.04401 kg + exchange(:) = GloCO2flux(:) * 1.157e-8_WP * 0.04401_WP ! [kg m⁻² s⁻¹] #endif else print *, 'not installed yet or error in cpl_oasis3mct_send', mype @@ -394,7 +399,9 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) #if defined (__recom) elseif (i.eq.16) then if (action) then - LocAtmCO2(1) = exchange(1) ! Atmospheric CO2 [uatm] + ! Convert mass mixing ratio (kg/kg) to mole fraction (dimensionless) + ! MW_CO2/MW_air = 44.01/28.97 + x_co2atm(:) = exchange(:) * (28.97_WP/44.01_WP) ! [mole fraction] end if #endif #else ! oifs From ed7431d513561ce6328d5438571923169cb5d72b Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Tue, 22 Jul 2025 11:54:50 +0200 Subject: [PATCH 04/17] higher precision conversion --- src/gen_forcing_couple.F90 | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 84b94ae60..13b6c4e47 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -217,9 +217,9 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) ! Convert CO2 flux from mmol/(m² day) to kg/(m² s) ! GloCO2flux is in [mmol/m2/day], need [kg/m2/s] ! Conversion: mmol/day -> mol/s -> kg/s - ! 1 mmol/day = 1e-3 mol / (24*3600 s) = 1.157e-8 mol/s - ! 1 mol CO2 = 0.04401 kg - exchange(:) = GloCO2flux(:) * 1.157e-8_WP * 0.04401_WP ! [kg m⁻² s⁻¹] + ! 1 mmol/day = 1e-3 mol / (86400 s) = 1.157407407407407e-8 mol/s + ! 1 mol CO2 = 44.0095 g/mol = 0.0440095 kg/mol (NIST 2018) + exchange(:) = GloCO2flux(:) * 1.157407407407407e-8_WP * 0.0440095_WP ! [kg m⁻² s⁻¹] #endif else print *, 'not installed yet or error in cpl_oasis3mct_send', mype @@ -400,8 +400,9 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) elseif (i.eq.16) then if (action) then ! Convert mass mixing ratio (kg/kg) to mole fraction (dimensionless) - ! MW_CO2/MW_air = 44.01/28.97 - x_co2atm(:) = exchange(:) * (28.97_WP/44.01_WP) ! [mole fraction] + ! 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) ! [mole fraction] end if #endif #else ! oifs From 89bb9b2bd1543e167a9fbb0ba6be3d08abd63069 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 23 Jul 2025 10:38:34 +0200 Subject: [PATCH 05/17] fix git merge mistake --- config/bin_2p3z2d/namelist.recom | 4 ---- 1 file changed, 4 deletions(-) diff --git a/config/bin_2p3z2d/namelist.recom b/config/bin_2p3z2d/namelist.recom index 07f2d3f46..5d18e9295 100644 --- a/config/bin_2p3z2d/namelist.recom +++ b/config/bin_2p3z2d/namelist.recom @@ -49,11 +49,7 @@ useErosion = .false. NitrogenSS = .false. ! When set to true, external sources and sinks of nitrogen are activated (Riverine, aeolian and denitrification) useAeolianN = .false. ! When set to true, aeolian nitrogen deposition is activated firstyearoffesomcycle = 1958 ! The first year of the actual physical forcing (e.g. JRA-55) used -<<<<<<< HEAD -lastyearoffesomcycle = 2023 ! Last year of the actual physical forcing used -======= lastyearoffesomcycle = 2024 ! Last year of the actual physical forcing used ->>>>>>> origin/fesom2.6_recom_tra_diags numofCO2cycles = 1 ! Number of cycles of the forcing planned currentCO2cycle = 1 ! Which CO2 cycle we are currently running DIC_PI = .true. From d06a657039c34c3a45950b1b5b36bf4e92abcb67 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 23 Jul 2025 14:15:57 +0200 Subject: [PATCH 06/17] get recom to compile on gnu by removing unused and unallocated vars --- src/int_recom/recom_init.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/int_recom/recom_init.F90 b/src/int_recom/recom_init.F90 index 3348b512b..361980ef1 100644 --- a/src/int_recom/recom_init.F90 +++ b/src/int_recom/recom_init.F90 @@ -97,7 +97,6 @@ subroutine recom_init(tracers, partit, mesh) allocate(decayBenthos ( benthos_num )) ! [1/day] Decay rate of detritus in the benthic layer allocate(PAR3D ( nl-1, node_size )) - GloFeDust = 0.d0 AtmFeInput = 0.d0 GloNDust = 0.d0 @@ -131,8 +130,6 @@ subroutine recom_init(tracers, partit, mesh) LocBenthos = 0.d0 decayBenthos = 0.d0 - wFluxPhy = 0.d0 - wFluxDia = 0.d0 PAR3D = 0.d0 ! pco2surf = 0.d0 From 663a67cfb271c495bc01559e9c40cb1f21ae9bce Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 23 Jul 2025 14:16:15 +0200 Subject: [PATCH 07/17] default recom namelists for AWI-ESM3 --- config/namelist.recom.2p3z2d | 355 +++++++++++++++++++++++++++++++ config/namelist.tra.recom.2p3z2d | 90 ++++++++ 2 files changed, 445 insertions(+) create mode 100644 config/namelist.recom.2p3z2d create mode 100644 config/namelist.tra.recom.2p3z2d diff --git a/config/namelist.recom.2p3z2d b/config/namelist.recom.2p3z2d new file mode 100644 index 000000000..5d18e9295 --- /dev/null +++ b/config/namelist.recom.2p3z2d @@ -0,0 +1,355 @@ +! This is the namelist file for recom + +&nam_rsbc +fe_data_source ='Albani' +nm_fe_data_file ='/albedo/work/projects/p_pool_recom/input/mesh_CORE2_finaltopo_mean/DustClimMonthlyAlbani.nc' +nm_aen_data_file ='/albedo/work/projects/p_pool_recom/input/mesh_CORE2_finaltopo_mean/AeolianNitrogenDep.nc' +nm_river_data_file ='/albedo/work/projects/p_pool_recom/input/mesh_CORE2_finaltopo_mean/RiverineInput.nc' +nm_erosion_data_file ='/albedo/work/projects/p_pool_recom/input/mesh_CORE2_finaltopo_mean/ErosionInput.nc' +nm_co2_data_file ='/albedo/work/projects/p_pool_recom/input/mesh_CORE2_finaltopo_mean/MonthlyAtmCO2_gcb2023.nc' +/ + +&pavariables +use_REcoM =.true. +REcoM_restart =.true. + +bgc_num = 30 !22 !33 !24 !38 +diags3d_num = 28 ! Number of diagnostic 3d tracers to be saved +bgc_base_num = 22 ! standard tracers +VDet = 20.d0 ! Sinking velocity, constant through the water column and positive downwards +VDet_zoo2 = 200.d0 ! Sinking velocity, constant through the water column +VPhy = 0.d0 !!! If the number of sinking velocities are different from 3, code needs to be changed !!! +VDia = 0.d0 +VCocco = 0.d0 +allow_var_sinking = .true. +biostep = 1 ! Number of times biology should be stepped forward for each time step +REcoM_Geider_limiter = .false. ! Decides what routine should be used to calculate limiters in sms +REcoM_Grazing_Variable_Preference = .true. ! Decides if grazing should have preference for phyN or DiaN +Grazing_detritus = .true. +het_resp_noredfield = .true. ! Decides respiratation of copepod group +diatom_mucus = .true. ! Decides nutrient limitation effect on aggregation +O2dep_remin = .false. ! O2remin Add option for O2 dependency of organic matter remineralization +use_ballasting = .false. ! BALL +use_density_scaling = .false. ! BALL +use_viscosity_scaling = .false. ! BALL +OmegaC_diss = .false. ! DISS Use OmegaC from Mocsy to compute calcite dissolution (after Aumont et al. 2015 and Gehlen et al. 2007) +CO2lim = .false. ! CO2 dependence of growth and calcification +Diags = .true. +constant_CO2 = .true. +UseFeDust = .true. ! Turns dust input of iron off when set to.false. +UseDustClim = .true. +UseDustClimAlbani = .true. ! Use Albani dustclim field (If it is false Mahowald will be used) +use_photodamage = .true. ! use Alvarez et al (2018) for chlorophyll degradation +HetRespFlux_plus = .true. !MB More stable computation of zooplankton respiration fluxes adding a small number to HetN +REcoMDataPath = '/albedo/work/projects/MarESys/ogurses/input/mesh_CORE2_finaltopo_mean/' +restore_alkalinity = .true. +useRivers = .false. +useRivFe = .false. ! When set to true, riverine Fe source is activated +useErosion = .false. +NitrogenSS = .false. ! When set to true, external sources and sinks of nitrogen are activated (Riverine, aeolian and denitrification) +useAeolianN = .false. ! When set to true, aeolian nitrogen deposition is activated +firstyearoffesomcycle = 1958 ! The first year of the actual physical forcing (e.g. JRA-55) used +lastyearoffesomcycle = 2024 ! Last year of the actual physical forcing used +numofCO2cycles = 1 ! Number of cycles of the forcing planned +currentCO2cycle = 1 ! Which CO2 cycle we are currently running +DIC_PI = .true. +Nmocsy = 1 ! Length of the vector that is passed to mocsy (always one for recom) +recom_debug =.false. +ciso =.false. ! Main switch to enable/disable carbon isotopes (13|14C) +benthos_num = 4 ! number of benthic BGC tracers -> 8 if (ciso == .true.) otherwise -> 4 +use_MEDUSA = .false. ! Main switch for the sediment model MEDUSA +sedflx_num = 0 ! if 0: no file from MEDUSA is read but default sediment +bottflx_num = 4 ! if ciso&ciso_14: =8; if .not.ciso_14: =6; no ciso: =4 +use_atbox = .false. +add_loopback = .false. ! add loopback fluxes through rivers to the surface +lb_tscale = 1.0 ! /year: fraction of loopback fluxes yearly added to the surface +/ + +&pasinking +Vdet_a = 0.0288 ! [1/day] +Vcalc = 0.0144 ! [1/day] +/ + +&painitialization_N +cPhyN = 0.2d0 +cHetN = 0.2d0 +cZoo2N = 0.2d0 +/ + +&paArrhenius +recom_Tref = 288.15d0 ! [K] +C2K = 273.15d0 ! Conversion from degrees C to K +Ae = 4500.d0 ! [K] Slope of the linear part of the Arrhenius function +reminSi = 0.02d0 +k_o2_remin = 15.d0 ! NEW O2remin mmol m-3; Table 1 in Cram 2018 cites DeVries & Weber 2017 for a range of 0-30 mmol m-3 +/ + +&palimiter_function +NMinSlope = 50.d0 +SiMinSlope = 1000.d0 +NCmin = 0.04d0 !0.05d0 +NCmin_d = 0.04d0 !0.05d0 +NCmin_c = 0.04d0 ! NEW +SiCmin = 0.04d0 +k_Fe = 0.04d0 +k_Fe_d = 0.12d0 +k_Fe_c = 0.09d0 ! NEW +k_si = 4.d0 +P_cm = 3.0d0 ! [1/day] Rate of C-specific photosynthesis +P_cm_d = 3.5d0 +P_cm_c = 2.8d0 ! NEW +/ + +&palight_calculations +k_w = 0.04d0 ! [1/m] Light attenuation coefficient +a_chl = 0.03d0 ! [1/m * 1/(mg Chl)] Chlorophyll specific attenuation coefficients +/ + +&paphotosynthesis +alfa = 0.14d0 ! [(mmol C*m2)/(mg Chl*W*day)] +alfa_d = 0.19d0 ! An initial slope of the P-I curve +alfa_c = 0.10d0 ! NEW +parFrac = 0.43d0 +/ + +&paassimilation +V_cm_fact = 0.7d0 ! scaling factor for temperature dependent maximum of C-specific N-uptake +V_cm_fact_d = 0.7d0 +V_cm_fact_c = 0.7d0 ! NEW +NMaxSlope = 1000.d0 ! Max slope for limiting function +SiMaxSlope = 1000.d0 +NCmax = 0.2d0 ! [mmol N/mmol C] Maximum cell quota of nitrogen (N:C) +NCmax_d = 0.2d0 +NCmax_c = 0.15d0 ! NEW +SiCmax = 0.8d0 +NCuptakeRatio = 0.2d0 ! [mmol N/mmol C] Maximum uptake ratio of N:C +NCUptakeRatio_d = 0.2d0 +NCUptakeRatio_c = 0.2d0 ! NEW +SiCUptakeRatio = 0.2d0 +k_din = 0.55d0 ! [mmol N/m3] Half-saturation constant for nitrate uptake +k_din_d = 1.0d0 +k_din_c = 0.9d0 ! NEW +Chl2N_max = 3.15d0 ! [mg CHL/mmol N] Maximum CHL a : N ratio = 0.3 gCHL gN^-1 +Chl2N_max_d = 4.2d0 +Chl2N_max_c = 3.5d0 ! NEW +res_phy = 0.01d0 ! [1/day] Maintenance respiration rate constant +res_phy_d = 0.01d0 +res_phy_c = 0.01d0 ! NEW +biosynth = 2.33d0 ! [mmol C/mmol N] Cost of biosynthesis +biosynthSi = 0.d0 +/ + +&pairon_chem +totalligand = 1.d0 ! [mumol/m3] order 1. Total free ligand +ligandStabConst = 100.d0 ! [m3/mumol] order 100. Ligand-free iron stability constant +/ + +&pazooplankton +graz_max = 0.31d0 ! [mmol N/(m3 * day)] Maximum grazing loss parameter +epsilonr = 0.09d0 ! [(mmol N)2 /m6] Half saturation constant for grazing loss +res_het = 0.028d0 ! [1/day] Respiration by heterotrophs and mortality (loss to detritus) +Redfield = 6.625 ! [mmol C/mmol N] Redfield ratio of C:N = 106:16 +loss_het = 0.04d0 ! [1/day] Temperature dependent N degradation of extracellular organic N (EON) +pzDia = 1.0d0 !0.5d0 ! Maximum diatom preference +sDiaNsq = 0.d0 +pzPhy = 0.5d0 !0.25d0 !1.0d0 ! Maximum nano-phytoplankton preference (NEW: 3/12) +sPhyNsq = 0.d0 +pzCocco = 0.666d0 ! NEW (8/12) +sCoccoNsq = 0.d0 ! NEW +pzMicZoo = 1.0d0 ! NEW 3Zoo Maximum nano-phytoplankton preference +sMicZooNsq = 0.d0 ! NEW 3Zoo +tiny_het = 1.d-5 ! for more stable computation of HetRespFlux (_plus). Value can be > tiny because HetRespFlux ~ hetC**2. +/ + +&pasecondzooplankton +graz_max2 = 0.1d0 ! [mmol N/(m3 * day)] Maximum grazing loss parameter +epsilon2 = 0.0144d0 ! [(mmol N)2 /m6] Half saturation constant for grazing loss +res_zoo2 = 0.0107d0 ! [1/day] Respiration by heterotrophs and mortality (loss to detritus) +loss_zoo2 = 0.003d0 ! [1/day] Temperature dependent N degradation of extracellular organic N (EON) +fecal_rate_n = 0.104d0 ! [1/day] Temperature dependent N degradation of \ +fecal_rate_c = 0.236d0 +fecal_rate_n_mes = 0.25d0 ! NEW 3Zoo +fecal_rate_c_mes = 0.32d0 ! NEW 3Zoo +pzDia2 = 1.5d0 !1.d0 ! Maximum diatom preference +sDiaNsq2 = 0.d0 +pzPhy2 = 0.5d0 ! Maximum diatom preference +sPhyNsq2 = 0.d0 +pzCocco2 = 0.5d0 ! NEW +sCoccoNsq2 = 0.d0 ! NEW +pzHet = 1.5d0 !0.8d0 ! Maximum diatom preference +sHetNsq = 0.d0 +pzMicZoo2 = 1.0d0 ! NEW 3Zoo Maximum nano-phytoplankton preference +sMicZooNsq2 = 0.d0 +t1_zoo2 = 28145.d0 ! Krill temp. function constant1 +t2_zoo2 = 272.5d0 ! Krill temp. function constant2 +t3_zoo2 = 105234.d0 ! Krill temp. function constant3 +t4_zoo2 = 274.15d0 ! Krill temp. function constant3 +/ + +&pathirdzooplankton +graz_max3 = 0.46d0 ! NEW 3Zoo [mmol N/(m3 * day)] Maximum grazing loss parameter +epsilon3 = 0.64d0 ! NEW 3Zoo [(mmol N)2 /m6] Half saturation constant for grazing loss +loss_miczoo = 0.01d0 ! NEW 3Zoo [1/day] Temperature dependent N degradation of extracellular organic N (EON) +res_miczoo = 0.01d0 ! NEW 3Zoo [1/day] Respiration by heterotrophs and mortality (loss to detritus) +pzDia3 = 0.5d0 ! NEW 3Zoo Maximum diatom preference +sDiaNsq3 = 0.d0 ! NEW 3Zoo +pzPhy3 = 1.0d0 ! NEW 3Zoo Maximum nano-phytoplankton preference +sPhyNsq3 = 0.d0 ! NEW 3Zoo +pzCocco3 = 0.d0 ! NEW 3Zoo Maximum coccolithophore preference ! ATTENTION: This value needs to be tuned; I start with zero preference! +sCoccoNsq3 = 0.d0 ! NEW 3Zoo +/ + +&pagrazingdetritus +pzDet = 0.5d0 ! Maximum small detritus prefence by first zooplankton +sDetNsq = 0.d0 +pzDetZ2 = 0.5d0 ! Maximum large detritus preference by first zooplankton +sDetZ2Nsq = 0.d0 +pzDet2 = 0.5d0 ! Maximum small detritus prefence by second zooplankton +sDetNsq2 = 0.d0 +pzDetZ22 = 0.5d0 ! Maximum large detritus preference by second zooplankton +sDetZ2Nsq2 = 0.d0 +/ + +&paaggregation +agg_PD = 0.165d0 ! [m3/(mmol N * day)] Maximum aggregation loss parameter for DetN +agg_PP = 0.015d0 ! [m3/(mmol N * day)] Maximum aggregation loss parameter for PhyN and DiaN (plankton) +/ + +&padin_rho_N +rho_N = 0.11d0 ! [1/day] Temperature dependent N degradation of extracellular organic N (EON) (Remineralization of DON) +/ + +&padic_rho_C1 +rho_C1 = 0.1d0 ! [1/day] Temperature dependent C degradation of extracellular organic C (EOC) +/ + +&paphytoplankton_N +lossN = 0.05d0 ! [1/day] Phytoplankton loss of organic N compounds +lossN_d = 0.05d0 +lossN_c = 0.05d0 ! NEW +/ + +&paphytoplankton_C +lossC = 0.10d0 ! [1/day] Phytoplankton loss of carbon +lossC_d = 0.10d0 +lossC_c = 0.10d0 ! NEW +/ + +&paphytoplankton_ChlA +deg_Chl = 0.25d0 !0.2d0 !0.25d0 ! [1/day] +deg_Chl_d = 0.15d0 !0.2d0 !0.15d0 +deg_Chl_c = 0.2d0 ! NEW (has been 0.5) +/ + +&padetritus_N +gfin = 0.3d0 ! NEW 3Zoo [] Grazing efficiency (fraction of grazing flux into zooplankton pool) +grazEff2 = 0.8d0 ! [] Grazing efficiency (fraction of grazing flux into second zooplankton pool) +grazEff3 = 0.8d0 ! NEW 3Zoo [] Grazing efficiency (fraction of grazing flux into microzooplankton pool) +reminN = 0.165d0 ! [1/day] Temperature dependent remineralisation rate of detritus +/ + +&padetritus_C +reminC = 0.15d0 ! [1/day] Temperature dependent remineralisation rate of detritus +rho_c2 = 0.1d0 ! [1/day] Temperature dependent C degradation of TEP-C +/ + +&paheterotrophs +lossN_z = 0.1d0 +lossC_z = 0.1d0 +/ + +&paseczooloss +lossN_z2 = 0.02d0 +lossC_z2 = 0.02d0 +/ + +&pathirdzooloss +lossN_z3 = 0.05d0 ! NEW 3Zoo +lossC_z3 = 0.05d0 ! NEW 3Zoo +/ + +&paco2lim ! NEW +Cunits = 976.5625 ! Conversion factor between [mol/m3] (model) and [umol/kg] (function): (1000 * 1000) / 1024 +a_co2_phy = 1.162e+00 ! [dimensionless] +a_co2_dia = 1.040e+00 ! [dimensionless] +a_co2_cocco = 1.109e+00 ! [dimensionless] +a_co2_calc = 1.102e+00 ! [dimensionless] +b_co2_phy = 4.888e+01 ! [mol/kg] +b_co2_dia = 2.890e+01 ! [mol/kg] +b_co2_cocco = 3.767e+01 ! [mol/kg] +b_co2_calc = 4.238e+01 ! [mol/kg] +c_co2_phy = 2.255e-01 ! [kg/mol] +c_co2_dia = 8.778e-01 ! [kg/mol] +c_co2_cocco = 3.912e-01 ! [kg/mol] +c_co2_calc = 7.079e-01 ! [kg/mol] +d_co2_phy = 1.023e+07 ! [kg/mol] +d_co2_dia = 2.640e+06 ! [kg/mol] +d_co2_cocco = 9.450e+06 ! [kg/mol] +d_co2_calc = 1.343e+07 ! [kg/mol] +/ + +&pairon +Fe2N = 0.033d0 ! Fe2C * 6.625 +Fe2N_benthos = 0.15d0 ! test, default was 0.14 Fe2C_benthos * 6.625 - will have to be tuned. [umol/m2/day] +kScavFe = 0.07d0 +dust_sol = 0.02d0 ! Dissolution of Dust for bioavaliable +RiverFeConc = 100 +/ + +&pacalc +calc_prod_ratio = 0.02 +calc_diss_guts = 0.0d0 +calc_diss_rate = 0.005714 ! 20.d0/3500.d0 +calc_diss_rate2 = 0.005714d0 +calc_diss_omegac = 0.197d0 ! NEW DISS Value from Aumont et al. 2015, will be used with OmegaC_diss flag +calc_diss_exp = 1.d0 ! NEW DISS Exponent in the dissolution rate of calcite, will be used with OmegaC_diss flag +/ + +&pabenthos_decay_rate +decayRateBenN = 0.005d0 +decayRateBenC = 0.005d0 +decayRateBenSi = 0.005d0 +q_NC_Denit = 0.86d0 ! N:C quota of the denitrification process +/ + +&paco2_flux_param +permil = 0.000000976 ! 1.e-3/1024.5d0 ! Converting DIC from [mmol/m3] to [mol/kg] +permeg = 1.e-6 ! [atm/uatm] Changes units from uatm to atm +!X1 = exp(-5.d0*log(10.d0)) ! Lowest ph-value = 7.7 (phlo) +!X2 = exp(-9.d0*log(10.d0)) ! Highest ph-value = 9.5 (phhi) +Xacc = 1.e-12 ! Accuracy for ph-iteration (phacc) +CO2_for_spinup = 278.d0 ! [uatm] Atmospheric partial pressure of CO2 +/ + +&paalkalinity_restoring +surf_relax_Alk = 3.2e-07 !10.d0/31536000.d0 +/ + +&paballasting +rho_POC = 1033.d0 ! kg m-3; density of POC (see Table 1 in Cram et al., 2018) +rho_PON = 1033.d0 ! kg m-3; density of PON (see Table 1 in Cram et al., 2018) +rho_CaCO3 = 2830.d0 ! kg m-3; density of CaCO3 (see Table 1 in Cram et al., 2018) +rho_opal = 2090.d0 ! kg m-3; density of Opal (see Table 1 in Cram et al., 2018) +rho_ref_part = 1230.d0 ! kg m-3; reference particle density (see Cram et al., 2018) +rho_ref_water = 1027.d0 ! kg m-3; reference seawater density (see Cram et al., 2018) +visc_ref_water = 0.00158d0 ! kg m-1 s-1; reference seawater viscosity, at Temp=4 degC (see Cram et al., 2018) +w_ref1 = 10.d0 ! m s-1; reference sinking velocity of small detritus +w_ref2 = 200.d0 ! m s-1; reference sinking velocity of large detritus +depth_scaling1 = 0.015d0 ! s-1; factor to increase sinking speed of det1 with depth, set to 0 if not wanted +depth_scaling2 = 0.d0 ! s-1; factor to increase sinking speed of det2 with depth, set to 0 if not wanted +max_sinking_velocity = 250.d0 ! d-1; for numerical stability, set a maximum possible sinking velocity here (applies to both detritus classes) +/ + +&paciso +ciso_init = .false. ! initial fractionation of bulk organic matter +ciso_14 = .false. ! include inorganic radiocarbon +ciso_organic_14 = .false. ! include organic radiocarbon +lambda_14 = 3.8561e-12 ! corresponding to 1 year = 365.00 days +delta_CO2_13 = -6.61 ! atmospheric d13C (permil), global-mean value +big_delta_CO2_14(1) = 0. ! atmospheric D14C (permil), northern hemisphere polewards of 30°N +big_delta_CO2_14(2) = 0. ! atmospheric D14C (permil), (sub) tropical zone 30°N - 30°S +big_delta_CO2_14(3) = 0. ! atmospheric D14C (permil), southern hemisphere polewards of 30°S +atbox_spinup = .false. +cosmic_14_init = 2.0 +/ + diff --git a/config/namelist.tra.recom.2p3z2d b/config/namelist.tra.recom.2p3z2d new file mode 100644 index 000000000..d1733bb91 --- /dev/null +++ b/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. +/ From 2372cc6175c698794c3122ea4b620170f93cb542 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 23 Jul 2025 14:16:44 +0200 Subject: [PATCH 08/17] switch to sea ice masked co2 flux --- src/gen_forcing_couple.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 13b6c4e47..98a84a79e 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -99,7 +99,7 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) use gen_bulk use force_flux_consv_interface #if defined (__recom) - use recom_modules, only: x_co2atm, GloCO2flux + use REcoM_GloVar, only: x_co2atm, GloCO2flux_seaicemask #endif implicit none @@ -215,11 +215,11 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) #if defined (__recom) elseif (i.eq.8) then ! Convert CO2 flux from mmol/(m² day) to kg/(m² s) - ! GloCO2flux is in [mmol/m2/day], need [kg/m2/s] + ! GloCO2flux_seaicemask is in [mmol/m2/day], need [kg/m2/s] ! Conversion: mmol/day -> mol/s -> kg/s ! 1 mmol/day = 1e-3 mol / (86400 s) = 1.157407407407407e-8 mol/s ! 1 mol CO2 = 44.0095 g/mol = 0.0440095 kg/mol (NIST 2018) - exchange(:) = GloCO2flux(:) * 1.157407407407407e-8_WP * 0.0440095_WP ! [kg m⁻² s⁻¹] + exchange(:) = GloCO2flux_seaicemask(:) * 1.157407407407407e-8_WP * 0.0440095_WP ! [kg m⁻² s⁻¹] #endif else print *, 'not installed yet or error in cpl_oasis3mct_send', mype From 398831e7516d45bbcf736d11a0b6864e4c593f31 Mon Sep 17 00:00:00 2001 From: chrisdane Date: Wed, 23 Jul 2025 14:28:23 +0200 Subject: [PATCH 09/17] add xco2 from oifs --- src/int_recom/recom_init.F90 | 7 ++++++- src/int_recom/recom_main.F90 | 1 + 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/src/int_recom/recom_init.F90 b/src/int_recom/recom_init.F90 index 361980ef1..fef4b667e 100644 --- a/src/int_recom/recom_init.F90 +++ b/src/int_recom/recom_init.F90 @@ -96,7 +96,7 @@ subroutine recom_init(tracers, partit, mesh) allocate(LocBenthos ( benthos_num )) allocate(decayBenthos ( benthos_num )) ! [1/day] Decay rate of detritus in the benthic layer allocate(PAR3D ( nl-1, node_size )) - + GloFeDust = 0.d0 AtmFeInput = 0.d0 GloNDust = 0.d0 @@ -131,6 +131,11 @@ subroutine recom_init(tracers, partit, mesh) LocBenthos = 0.d0 decayBenthos = 0.d0 PAR3D = 0.d0 + + if (use_atbox) then + allocate(x_co2atm ( node_size )) + x_co2atm = 0.d0 + endif ! pco2surf = 0.d0 ! dflux = 0.d0 diff --git a/src/int_recom/recom_main.F90 b/src/int_recom/recom_main.F90 index 452cb1b6d..1a3803e3e 100755 --- a/src/int_recom/recom_main.F90 +++ b/src/int_recom/recom_main.F90 @@ -208,6 +208,7 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) LocAtmCO2_14 = AtmCO2_14(lat_zone(lat_val), month) end if end if + LocAtmCO2 = x_co2atm(n) ! from oifs end if ! use_atbox if (ciso) then From e885e92b3080937983045309a5c916a7dafcc274 Mon Sep 17 00:00:00 2001 From: chrisdane Date: Wed, 23 Jul 2025 14:36:01 +0200 Subject: [PATCH 10/17] corrected if atmbox case --- src/int_recom/recom_init.F90 | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/int_recom/recom_init.F90 b/src/int_recom/recom_init.F90 index fef4b667e..336247eb2 100644 --- a/src/int_recom/recom_init.F90 +++ b/src/int_recom/recom_init.F90 @@ -133,9 +133,12 @@ subroutine recom_init(tracers, partit, mesh) PAR3D = 0.d0 if (use_atbox) then + allocate(x_co2atm ( 1 )) + else allocate(x_co2atm ( node_size )) - x_co2atm = 0.d0 endif + x_co2atm = 0.d0 + ! pco2surf = 0.d0 ! dflux = 0.d0 From 038e78170aac903fa169d998d26b8c9ae12c35e1 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 6 Aug 2025 09:05:30 +0200 Subject: [PATCH 11/17] change dim of wf to match number of tracers depending on recom setting --- src/int_recom/recom_extra.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/int_recom/recom_extra.F90 b/src/int_recom/recom_extra.F90 index 441fbd724..363438443 100644 --- a/src/int_recom/recom_extra.F90 +++ b/src/int_recom/recom_extra.F90 @@ -23,7 +23,7 @@ subroutine Depth_calculations(n, nn, wf, zf, thick, recipthick, partit, mesh) integer , intent(in) :: n ! Current node ! Output arrays - real(kind=8), dimension(mesh%nl,5), intent(out) :: wf ! [m/day] Flux velocities at the border of the control volumes + real(kind=8), dimension(mesh%nl,6), intent(out) :: wf ! [m/day] Flux velocities at the border of the control volumes real(kind=8), dimension(mesh%nl), intent(out) :: zf ! [m] Depth of vertical fluxes real(kind=8), dimension(mesh%nl-1), intent(out) :: thick ! [m] Distance between two nodes = layer thickness real(kind=8), dimension(mesh%nl-1), intent(out) :: recipthick ! [1/m] Reciprocal thickness From b5c6725fe0842676f20da0b96a91c8d42cf9a1ef Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 6 Aug 2025 14:42:26 +0200 Subject: [PATCH 12/17] add xco2atm output option --- config/namelist.io.recom.2p3z2d | 95 +++++++++++++++++++++++++++++++++ src/io_meandata.F90 | 5 ++ 2 files changed, 100 insertions(+) create mode 100644 config/namelist.io.recom.2p3z2d diff --git a/config/namelist.io.recom.2p3z2d b/config/namelist.io.recom.2p3z2d new file mode 100644 index 000000000..8381d77c2 --- /dev/null +++ b/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, +/ diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 57492a64c..7efa667d6 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -504,6 +504,11 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) ! output RECOM 2D #if defined(__recom) +CASE ('xCO2atm ') + if (use_REcoM) then + call def_stream(nod2D, myDim_nod2D, 'xCO2atm', 'atmospheric CO2 mass mixing ratio', 'mole fraction', x_co2atm(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) + end if + CASE ('dpCO2s ') if (use_REcoM) then call def_stream(nod2D, myDim_nod2D, 'dpCO2s', 'Difference of oceanic pCO2 minus atmospheric pCO2', 'uatm', GlodPCO2surf(:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) From 6e19194eafc8a326c6bf2481e56db39123ce0f87 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 6 Aug 2025 14:43:15 +0200 Subject: [PATCH 13/17] multiply xco2 by 1e6 to go from mole fraction to ppm --- src/gen_forcing_couple.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index 98a84a79e..d8fce6af7 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -402,7 +402,7 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) ! 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) ! [mole fraction] + x_co2atm(:) = exchange(:) * ((28.9647_WP/44.0095_WP)*1e6_WP) ! [mole fraction] end if #endif #else ! oifs From 20e79e739e3b3c68e4e0b0670d972b62607fa8c4 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 6 Aug 2025 17:25:47 +0200 Subject: [PATCH 14/17] co2 correct flux unit conversion --- src/gen_forcing_couple.F90 | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index d8fce6af7..33c4e1d20 100755 --- a/src/gen_forcing_couple.F90 +++ b/src/gen_forcing_couple.F90 @@ -214,12 +214,10 @@ subroutine update_atm_forcing(istep, ice, tracers, dynamics, partit, mesh) end do #if defined (__recom) elseif (i.eq.8) then - ! Convert CO2 flux from mmol/(m² day) to kg/(m² s) - ! GloCO2flux_seaicemask is in [mmol/m2/day], need [kg/m2/s] - ! Conversion: mmol/day -> mol/s -> kg/s - ! 1 mmol/day = 1e-3 mol / (86400 s) = 1.157407407407407e-8 mol/s + ! 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.157407407407407e-8_WP * 0.0440095_WP ! [kg m⁻² s⁻¹] + 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 From f8224a83686fab8ec1d0d4e86c2016010f8c0e22 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Tue, 2 Sep 2025 09:55:02 +0200 Subject: [PATCH 15/17] merge main into branch --- src/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index bab416317..c2b7b067a 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -250,7 +250,7 @@ if(OPENMP_REPRODUCIBLE) endif() if(${RECOM_COUPLED}) - target_compile_definitions(${PROJECT_NAME} PRIVATE __recom USE_PRECISION=2 __3Zoo2Det __coccos)# __usetp) + target_compile_definitions(${PROJECT_NAME} PRIVATE __recom USE_PRECISION=2 __3Zoo2Det)# __usetp) endif() if(${CISO_COUPLED}) From 2fb0cb55b554ed6702b6c786a64a043292477e04 Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Tue, 18 Nov 2025 15:49:22 +0100 Subject: [PATCH 16/17] saving wip, recomp + cavity runs, judiths suggestions not included yet --- src/fesom_module.F90 | 1 + src/gen_ic3d.F90 | 28 ++++++++++++++++++--- src/int_recom/recom_forcing.F90 | 44 +++++++++++++++++++++++++++------ src/int_recom/recom_main.F90 | 7 ++++-- src/oce_ale.F90 | 18 ++++++++++++++ 5 files changed, 86 insertions(+), 12 deletions(-) diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 50e379133..7ecc72a3c 100755 --- a/src/fesom_module.F90 +++ b/src/fesom_module.F90 @@ -38,6 +38,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) diff --git a/src/gen_ic3d.F90 b/src/gen_ic3d.F90 index 4af214cd6..732777d0d 100644 --- a/src/gen_ic3d.F90 +++ b/src/gen_ic3d.F90 @@ -412,7 +412,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) @@ -423,7 +425,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 @@ -481,6 +485,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 @@ -497,6 +509,7 @@ SUBROUTINE do_ic3d(tracers, partit, mesh) !! ** Purpose : read 3D initial conditions for tracers from netcdf and interpolate on model grid !!---------------------------------------------------------------------- USE insitu2pot_interface + USE ieee_arithmetic IMPLICIT NONE type(t_mesh), intent(in), target :: mesh type(t_partit), intent(inout), target :: partit @@ -521,7 +534,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!" @@ -534,6 +551,11 @@ SUBROUTINE do_ic3d(tracers, partit, mesh) DEALLOCATE(bilin_indx_i, bilin_indx_j) 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) diff --git a/src/int_recom/recom_forcing.F90 b/src/int_recom/recom_forcing.F90 index 59bef0116..52d434e01 100644 --- a/src/int_recom/recom_forcing.F90 +++ b/src/int_recom/recom_forcing.F90 @@ -185,15 +185,38 @@ subroutine REcoM_Forcing(zNodes, n, Nn, state, SurfSW, Loc_slp, Temp, Sali, Sali stop endif - call flxco2(co2flux, co2ex, dpco2surf, & - ph, pco2surf, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis, K0, & - REcoM_T, REcoM_S, REcoM_Alk, REcoM_DIC, REcoM_Si, REcoM_Phos, kw660, LocAtmCO2, Patm, thick(One), Nmocsy, Lond,Latd, & - optCON='mol/m3',optT='Tpot ',optP='m ',optB='u74',optK1K2='l ',optKf='dg',optGAS='Pinsitu',optS='Sprc') + !!---- Skip CO2 flux calculation in ice shelf cavities (ulevels_nod2d > 1) + !!---- Cavity nodes have uninitialized biogeochemical tracers + !!---- NOTE: ulevels_nod2d is already indexed by local node numbers + if (mesh%ulevels_nod2d(n) > 1) then + !! Set CO2 flux to zero in cavity regions + co2flux(1) = 0.0d0 + co2ex(1) = 0.0d0 + dpco2surf(1) = 0.0d0 + ph(1) = 8.0d0 + pco2surf(1) = 0.0d0 + fco2(1) = 0.0d0 + co2(1) = 0.0d0 + hco3(1) = 0.0d0 + co3(1) = 0.0d0 + OmegaA(1) = 0.0d0 + OmegaC(1) = 0.0d0 + BetaD(1) = 0.0d0 + rhoSW(1) = 1027.0d0 + else + !! Calculate CO2 flux for non-cavity nodes + call flxco2(co2flux, co2ex, dpco2surf, & + ph, pco2surf, fco2, co2, hco3, co3, OmegaA, OmegaC, BetaD, rhoSW, p, tempis, K0, & + REcoM_T, REcoM_S, REcoM_Alk, REcoM_DIC, REcoM_Si, REcoM_Phos, kw660, LocAtmCO2, Patm, thick(One), Nmocsy, Lond,Latd, & + optCON='mol/m3',optT='Tpot ',optP='m ',optB='u74',optK1K2='l ',optKf='dg',optGAS='Pinsitu',optS='Sprc') + endif ! changed optK1K2='l ' to 'm10' if((co2flux(1)>1.e10) .or. (co2flux(1)<-1.e10)) then ! co2flux(1)=0.0 print*, 'ERROR: co2 flux !' + print*, 'ulevels_nod2d(n): ',mesh%ulevels_nod2d(n) + print*, 'node n (local): ',n print*, 'pco2surf: ',pco2surf print*, 'co2: ',co2 print*, 'rhoSW: ', rhoSW @@ -230,9 +253,16 @@ subroutine REcoM_Forcing(zNodes, n, Nn, state, SurfSW, Loc_slp, Temp, Sali, Sali ppo = Loc_slp/Pa2atm !1 !slp divided by 1 atm REcoM_O2 = max(tiny*1e-3,state(one,ioxy)*1e-3) ! convert from mmol/m3 to mol/m3 for mocsy - call o2flux(REcoM_T, REcoM_S, kw660, ppo, REcoM_O2, Nmocsy, o2ex) - oflux = o2ex * 1.e3 *SecondsPerDay !* (1.d0 - Loc_ice_conc) [mmol/m2/d] - o2flux_seaicemask = o2ex * 1.e3 ! back to mmol here [mmol/m2/s] + !!---- Skip O2 flux calculation in ice shelf cavities + if (mesh%ulevels_nod2d(n) > 1) then + o2ex(1) = 0.0d0 + oflux = 0.0d0 + o2flux_seaicemask = 0.0d0 + else + call o2flux(REcoM_T, REcoM_S, kw660, ppo, REcoM_O2, Nmocsy, o2ex) + oflux = o2ex * 1.e3 *SecondsPerDay !* (1.d0 - Loc_ice_conc) [mmol/m2/d] + o2flux_seaicemask = o2ex * 1.e3 ! back to mmol here [mmol/m2/s] + endif ! Source-Minus-Sinks diff --git a/src/int_recom/recom_main.F90 b/src/int_recom/recom_main.F90 index 1a3803e3e..6b72c9c2d 100755 --- a/src/int_recom/recom_main.F90 +++ b/src/int_recom/recom_main.F90 @@ -57,6 +57,7 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) use o_PARAM use MOD_PARTIT use MOD_PARSUP + use ieee_arithmetic use recom_declarations use bio_fluxes_interface @@ -153,7 +154,8 @@ subroutine recom(ice, dynamics, tracers, partit, mesh) !********************************* LOOP STARTS ***************************************** do n=1, myDim_nod2D ! needs exchange_nod in the end -! if (ulevels_nod2D(n)>1) cycle + !!---- Skip cavity nodes - REcoM not supported in ice shelf cavities + if (ulevels_nod2D(n)>1) cycle ! nzmin = ulevels_nod2D(n) !!---- Number of vertical layers @@ -636,7 +638,8 @@ subroutine bio_fluxes(tracers, partit, mesh) use recom_locvar use recom_glovar use recom_config - + use o_PARAM + use ieee_arithmetic use mod_mesh use MOD_PARTIT use MOD_PARSUP diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index fe0051854..96517248d 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -3339,6 +3339,7 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) use solve_tracers_ale_interface use write_step_info_interface use check_blowup_interface + use ieee_arithmetic use fer_solve_interface use impl_vert_visc_ale_vtransp_interface #if defined (FESOM_PROFILING) @@ -3356,6 +3357,7 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) real(kind=8) :: t0,t1, t2, t30, t3, t4, t5, t6, t7, t8, t9, t10, loc, glo integer :: node integer :: nz, elem, nzmin, nzmax !for KE diagnostic + integer :: tr_num ! for cavity NaN cleanup !___________________________________________________________________________ ! pointer on necessary derived types real(kind=WP), dimension(:), pointer :: eta_n @@ -3792,6 +3794,22 @@ subroutine oce_timestep_ale(n, ice, dynamics, tracers, partit, mesh) call fesom_profiler_end("oce_tracer_solve") call fesom_profiler_start("oce_thickness_update") #endif +#if defined (__recom) + ! CAVITY FIX: Clean up NaN created by tracer advection/diffusion at cavity nodes + ! Cavity nodes can create NaN during numerical operations with zero/near-zero thickness + if (use_cavity) then + do node=1, partit%myDim_nod2D+partit%eDim_nod2D + if (mesh%ulevels_nod2D(node) > 1) then + ! Clean NaN in ALL levels for cavity nodes (not just cavity layers) + do tr_num=1, tracers%num_tracers + where (.not. ieee_is_finite(tracers%data(tr_num)%values(:, node))) + tracers%data(tr_num)%values(:, node) = 0.0_WP + end where + enddo + endif + enddo + endif +#endif !___________________________________________________________________________ ! Update hnode=hnode_new, helem From b557ca854c193465fbf7468f8db577c359d57f5b Mon Sep 17 00:00:00 2001 From: Jan Streffing Date: Wed, 3 Jun 2026 16:23:21 +0200 Subject: [PATCH 17/17] mask netcdf output by wet-cell topography instead of value==0 PR #799 set any accumulated output value below 1e-30 to _FillValue. That also masks legitimate ocean zeros: ice-free sea-ice concentration (a_ice and all other 2D ice variables) and vanishing 3D fields such as IDEMIX energy were written as NaN instead of 0, contrary to the CMIP/CMOR convention. The masking was applied uniformly to 2D and 3D fields. Replace the value-based test with a proper wet-cell mask derived from the mesh level arrays (ulevels/nlevels for elements, ulevels_nod2D/ nlevels_nod2D for nodes). Only cells above a cavity or below the bottom topography are written as _FillValue; every wet cell keeps its averaged value, including a genuine zero. 2D fields and non-spatial vertical axes (density / ice classes) are never masked. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/io_meandata.F90 | 76 ++++++++++++++++++++++++++++++++++++--------- 1 file changed, 61 insertions(+), 15 deletions(-) diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index fba00c3d8..d7809e000 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -1977,7 +1977,46 @@ subroutine update_means ! ! !_______________________________________________________________________________ -! main output routine called at the end of each time step --> here is decided if +! Return the vertically valid (wet) index range [ul_loc, kmax_loc] for column j of +! the given output stream. Cells with level index outside this range are below the +! bottom topography or above a cavity and are the only ones written as _FillValue. +! For 2D fields (nlev_loc==1) and non-spatial vertical axes (density / ice classes) +! the whole column is valid, so a legitimate zero (e.g. ice-free a_ice or vanishing +! IDEMIX energy) is kept as zero rather than turned into a missing value. +subroutine get_wet_range(entry, mesh, nlev_loc, j, ul_loc, kmax_loc) + use mod_mesh + implicit none + type(Meandata), intent(in) :: entry + type(t_mesh), intent(in) :: mesh + integer, intent(in) :: nlev_loc, j + integer, intent(out) :: ul_loc, kmax_loc + integer :: nl_bot + + if (nlev_loc == mesh%nl .or. nlev_loc == mesh%nl-1) then + ! genuine 3D field on the vertical grid --> mask by topography / cavity + if (entry%is_elem_based) then + ul_loc = mesh%ulevels(j) + nl_bot = mesh%nlevels(j) + else + ul_loc = mesh%ulevels_nod2D(j) + nl_bot = mesh%nlevels_nod2D(j) + end if + ! full levels (nz, interfaces) keep nl_bot levels, mid layers (nz1) keep nl_bot-1 + if (nlev_loc == mesh%nl) then + kmax_loc = nl_bot + else + kmax_loc = nl_bot - 1 + end if + else + ! 2D fields and non-spatial vertical axes: every entry is valid + ul_loc = 1 + kmax_loc = nlev_loc + end if +end subroutine +! +! +!_______________________________________________________________________________ +! main output routine called at the end of each time step --> here is decided if ! output event is triggered subroutine output(istep, ice, dynamics, tracers, partit, mesh) use g_clock @@ -1999,6 +2038,7 @@ subroutine output(istep, ice, dynamics, tracers, partit, mesh) logical, save :: lfirst=.true. integer :: n, k integer :: i, j !for OMP loops + integer :: nlev_loc, ul_loc, kmax_loc !for wet-cell masking of output logical :: do_output type(Meandata), pointer :: entry type(t_mesh), intent(in) , target :: mesh @@ -2139,36 +2179,42 @@ subroutine output(istep, ice, dynamics, tracers, partit, mesh) !___________________________________________________________________ ! write double precision output if (entry%accuracy == i_real8) then -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J) + nlev_loc = size(entry%local_values_r8,dim=1) +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,ul_loc,kmax_loc) DO J=1, size(entry%local_values_r8,dim=2) - DO I=1, size(entry%local_values_r8,dim=1) - ! Check if point has valid data (non-zero accumulated value) - ! Use small epsilon to account for floating point precision - if (abs(entry%local_values_r8(I,J)) < 1.0e-30_real64) then - entry%local_values_r8_copy(I,J) = NC_FILL_DOUBLE ! No data - set to fill value + ! Determine the vertically valid (wet) range for this column. Cells + ! outside [ul_loc,kmax_loc] are below the bottom topography or above + ! the cavity and have never been touched -> they are the only ones set + ! to _FillValue. Wet cells keep their (averaged) value, even when it is + ! a legitimate zero (e.g. ice-free a_ice or vanishing IDEMIX energy). + call get_wet_range(entry, mesh, nlev_loc, J, ul_loc, kmax_loc) + DO I=1, nlev_loc + if (I < ul_loc .or. I > kmax_loc) then + entry%local_values_r8_copy(I,J) = NC_FILL_DOUBLE ! dry cell - set to fill value else entry%local_values_r8_copy(I,J) = entry%local_values_r8(I,J) /real(entry%addcounter,real64) ! compute_means end if entry%local_values_r8(I,J) = 0._real64 ! clean_meanarrays - reset to 0 for next accumulation - END DO ! --> DO I=1, size(entry%local_values_r8,dim=1) + END DO ! --> DO I=1, nlev_loc END DO ! --> DO J=1, size(entry%local_values_r8,dim=2) !$OMP END PARALLEL DO !___________________________________________________________________ ! write single precision output else if (entry%accuracy == i_real4) then -!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J) + nlev_loc = size(entry%local_values_r4,dim=1) +!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(I,J,ul_loc,kmax_loc) DO J=1, size(entry%local_values_r4,dim=2) - DO I=1, size(entry%local_values_r4,dim=1) - ! Check if point has valid data (non-zero accumulated value) - ! Use small epsilon to account for floating point precision - if (abs(entry%local_values_r4(I,J)) < 1.0e-30_real32) then - entry%local_values_r4_copy(I,J) = NC_FILL_FLOAT ! No data - set to fill value + ! see comment in the double precision branch above + call get_wet_range(entry, mesh, nlev_loc, J, ul_loc, kmax_loc) + DO I=1, nlev_loc + if (I < ul_loc .or. I > kmax_loc) then + entry%local_values_r4_copy(I,J) = NC_FILL_FLOAT ! dry cell - set to fill value else entry%local_values_r4_copy(I,J) = entry%local_values_r4(I,J) /real(entry%addcounter,real32) ! compute_means end if entry%local_values_r4(I,J) = 0._real32 ! clean_meanarrays - reset to 0 for next accumulation - END DO ! --> DO I=1, size(entry%local_values_r4,dim=1) + END DO ! --> DO I=1, nlev_loc END DO ! --> DO J=1, size(entry%local_values_r4,dim=2) !$OMP END PARALLEL DO end if ! --> if (entry%accuracy == i_real8) then