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/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. +/ diff --git a/src/cpl_driver.F90 b/src/cpl_driver.F90 index 0fb3532cf..24a2538c9 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 @@ -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' diff --git a/src/fesom_module.F90 b/src/fesom_module.F90 index 9ea20678a..b391292f4 100755 --- a/src/fesom_module.F90 +++ b/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) diff --git a/src/gen_forcing_couple.F90 b/src/gen_forcing_couple.F90 index d5c78ad11..fdfb14365 100755 --- a/src/gen_forcing_couple.F90 +++ b/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 diff --git a/src/gen_ic3d.F90 b/src/gen_ic3d.F90 index 445928bed..67a3f7d1e 100644 --- a/src/gen_ic3d.F90 +++ b/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 diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index 007cf2490..166607550 100644 --- a/src/io_meandata.F90 +++ b/src/io_meandata.F90 @@ -832,6 +832,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) @@ -1002,7 +1007,7 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh) #endif !___________________________________________________________________________________________________________________________________ -!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 3D streams <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< +! 3D streams !___________________________________________________________________________________________________________________________________ CASE ('temp ') call def_stream((/nl-1, nod2D/), (/nl-1, myDim_nod2D/), 'temp', 'temperature', 'C', tracers%data(1)%values(:,:), io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh) @@ -2748,7 +2753,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 @@ -2770,6 +2814,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 @@ -2942,40 +2987,44 @@ subroutine output(istep, ice, dynamics, tracers, partit, mesh) ! data range wrecks GRIB packing precision (real ice values get ! quantized to ~0). Multio gets the clean mean (zeros stay zeros). 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) -#if defined(__MULTIO) - entry%local_values_r8_copy(I,J) = entry%local_values_r8(I,J) /real(entry%addcounter,real64) ! compute_means -#else - 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 #endif 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) -#if defined(__MULTIO) - entry%local_values_r4_copy(I,J) = entry%local_values_r4(I,J) /real(entry%addcounter,real32) ! compute_means -#else - 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 #endif 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 diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 0f8e93b57..086d21c82 100644 --- a/src/oce_ale.F90 +++ b/src/oce_ale.F90 @@ -3343,6 +3343,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) @@ -3360,6 +3361,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 @@ -3799,6 +3801,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 diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index 9c4fb7a37..39bcf7b19 100644 --- a/src/oce_ale_tracer.F90 +++ b/src/oce_ale_tracer.F90 @@ -235,7 +235,7 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) endif SinkingVel1 = 0.0d0 ! OG 16.03.23 SinkingVel2 = 0.0d0 ! OG 16.03.23 -#endif +#endif ! do tracer AB (Adams-Bashfort) interpolation only for advectiv part ! needed @@ -267,6 +267,24 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh) call save_cmor_advection(tr_num, tracers%work%del_ttf(:, 1:myDim_nod2D), myDim_nod2D, nl-1) end if +!if (.FALSE.) then +! O:G - tra_diag +!#if defined (__recom) +! if (tracers%data(tr_num)%ltra_diag) then +! do n=1, myDim_nod2D+eDim_nod2D +! nu1 = ulevels_nod2D(n) +! nl1 = nlevels_nod2D(n) +! do nz = nu1, nl1-1 + ! Horizontal advection part +! tracers%work%tra_advhoriz(nz,n,tr_num) = tracers%work%del_ttf_advhoriz(nz,n) + ! Vertical advection part +! tracers%work%tra_advvert (nz,n,tr_num) = tracers%work%del_ttf_advvert(nz,n) +! end do +! end do +! end if +!#endif +!endif + !___________________________________________________________________________ ! diffuse tracers if (flag_debug .and. mype==0) print *, achar(27)//'[37m'//' --> call diff_tracers_ale'//achar(27)//'[0m'