MODULE tracco2i_mod ! ! This module does the work for the interactive CO2 tracers ! Authors: Patricia Cadule and Olivier Boucher ! ! Purpose and description: ! ----------------------- ! Main routine for the interactive carbon cycle ! Gather all carbon fluxes and emissions from ORCHIDEE, PISCES and fossil fuel ! Compute the net flux in source field which is used in phytrac ! Compute global CO2 mixing ratio for radiation scheme if option is activated ! Redistribute CO2 evenly over the atmosphere if transport is desactivated ! ! Purpose: ! Interactive CO2 tracer module for LMDZ atmospheric model ! Supports emission-driven AND concentration-driven modes ! Supports daily emission files, AND monthly with SZ98 interpolation ! Supports 2D (lonxlat) correction fields for land and ocean ! ! Authors: ! Patricia Cadule (IPSL) and Olivier Boucher (IPSL/LMD) ! ! History: ! 2020: Initial version (P. Cadule and O. Boucher) ! 2026-01: Added concentration-driven mode + SZ98 interpolation (P. Cadule) ! 2026-02: Added 2D (lonxlat) correction fields for land and ocean (P. Cadule) ! 2026-08: the daily relaxation pulls the mass-weighted global mean ! towards its pre-industrial target by a uniform increment, so the carbon ! it supplies measures the disequilibrium the free run would otherwise ! drift by (P. Cadule) ! ! ! Operating Modes: ! 1. EMISSION-DRIVEN (carbon_cycle_conc_driven = .FALSE.) ! - CO2 evolves from surface fluxes and emissions ! - source(:,id_CO2) contains the net flux ! - Requires CO2 in tracer.def (co2_tracer_available = .TRUE.) ! ! 2. CONCENTRATION-DRIVEN (carbon_cycle_conc_driven = .TRUE.) ! - Atmospheric CO2 is prescribed from co2_ppm ! - source(:,id_CO2) = 0 (no flux feedback) ! - co2_send transmits prescribed CO2 to land/ocean models ! - Can operate WITHOUT CO2 in tracer.def (co2_tracer_available = .FALSE.) ! ! ! IMPORTANT: Transport Control ! carbon_cycle_tr does NOT control tracer transport in LMDZ dynamics ! Transport is controlled by the presence of CO2 in traceur.def and tracer.def ! ! For true concentration-driven mode (no transport CPU cost): ! - Remove CO2 from traceur.def and tracer.def ! - co2_tracer_available will be .FALSE. ! - co2_send will use co2_ppm directly ! ! Emission Reading Modes: ! (A) read_daily_co2ff = .TRUE. ==> Files have time = year_len (360/365/366) ! (B) read_daily_co2ff = .FALSE. ==> Files have time = 12, use SZ98 interpolation ! ! ! !=============================================================================== ! IMPLICIT NONE PRIVATE PUBLIC :: tracco2i_init, tracco2i !----------------------------------------------------------------------------- ! Module-level variables !----------------------------------------------------------------------------- INTEGER, SAVE :: day_pre = -1 !$OMP THREADPRIVATE(day_pre) LOGICAL, SAVE :: first_call = .TRUE. !$OMP THREADPRIVATE(first_call) LOGICAL, SAVE :: check_fCO2_nbp_in_cfname = .FALSE. !$OMP THREADPRIVATE(check_fCO2_nbp_in_cfname) !--------------------------------------------------------------------------- ! CO2 tracer availability flag !--------------------------------------------------------------------------- ! This flag indicates whether the CO2 tracer exists in the tracer arrays. ! When carbon_cycle_tr = .FALSE. or CO2 is not in tracer.def, the tracer ! arrays may not include CO2. Accessing tr_seri(:,:,id_CO2) with invalid ! id_CO2 would cause segmentation faults. ! ! The flag is set ONCE at first call based on: ! id_CO2 >= 1 .AND. id_CO2 <= nbtr !--------------------------------------------------------------------------- LOGICAL, SAVE :: co2_tracer_available = .FALSE. !$OMP THREADPRIVATE(co2_tracer_available) CONTAINS !=============================================================================== ! SUBROUTINE: tracco2i_init ! ! Purpose: ! Initialise the carbon cycle module. Called from phytrac_init. ! Must be called BEFORE first call to phys_output_write. !=============================================================================== ! << PC !SUBROUTINE tracco2i_init() ! ! USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, & ! carbon_cycle_tr, carbon_cycle_rad, & ! carbon_cycle_conc_driven ! ! IMPLICIT NONE ! ! ! Call carbon_cycle_init if any carbon cycle pathway is active ! ! This allocates all necessary arrays ! IF (carbon_cycle_cpl .OR. carbon_cycle_tr .OR. & ! carbon_cycle_rad .OR. carbon_cycle_conc_driven) THEN ! CALL carbon_cycle_init() ! END IF ! !END SUBROUTINE tracco2i_init SUBROUTINE tracco2i_init() USE infotrac_phy, ONLY: tracers USE strings_mod, ONLY: strIdx, maxlen USE carbon_cycle_mod, ONLY: carbon_cycle_init, carbon_cycle_cpl, & carbon_cycle_tr, carbon_cycle_rad, & carbon_cycle_conc_driven USE carbon_cycle_mod, ONLY: id_CO2 IMPLICIT NONE !--- Names of the tracers actually carried by the physics, in tr_seri order CHARACTER(LEN=maxlen), ALLOCATABLE :: tnames(:) !--- Resolve the CO2 slot on the packed physics list. tr_seri holds only ! the tracers flagged isInPhysics by init_infotrac_phy, in table order, ! so the index must be taken on that packed list, never on the full ! tracers table. strIdx returns 0 when the name is absent, which is the ! value the safety check further down expects for a run without a CO2 ! tracer. tnames = PACK(tracers(:)%name, MASK=tracers(:)%isInPhysics) id_CO2 = strIdx(tnames, 'CO2') ! Call carbon_cycle_init if any carbon cycle pathway is active ! This allocates all necessary arrays IF (carbon_cycle_cpl .OR. carbon_cycle_tr .OR. & carbon_cycle_rad .OR. carbon_cycle_conc_driven) THEN CALL carbon_cycle_init() END IF END SUBROUTINE tracco2i_init ! >> PC !=============================================================================== ! SUBROUTINE: tracco2i ! ! Purpose: ! Main routine for interactive carbon cycle processing ! Handles both emission-driven and concentration-driven modes ! ! Description: ! 1. Gathers surface CO2 fluxes from Land (ORCHIDEE), Ocean (PISCES), ! Fossil Fuel (and Biomass Burning). ! 2. Builds the tracer source term [kg CO2 m-2 s-1]. ! 3. Computes global mean CO2 mixing ratio and atmospheric mass. ! 4. Handles scenario-based target CO2 (diagnostic or emission-comp). ! ! Authors: ! 2026: Patricia Cadule ! 2020: Patricia Cadule, Olivier Boucher (initial) ! ! ! Called from: ! phytrac ! !=============================================================================== SUBROUTINE tracco2i(pdtphys, debutphy, & xlat, xlon, pphis, pphi, & t_seri, pplay, paprs, tr_seri, source) USE dimphy USE infotrac_phy, ONLY: nbtr USE geometry_mod, ONLY: cell_area ! Explicit import of needed variables from carbon_cycle_mod ! << PC USE carbon_cycle_mod, ONLY: id_CO2 !--- Global-mean relaxation controls USE carbon_cycle_mod, ONLY: carbon_cycle_nudge USE carbon_cycle_mod, ONLY: carbon_cycle_nudge_target USE carbon_cycle_mod, ONLY: carbon_cycle_nudge_tau USE carbon_cycle_mod, ONLY: carbon_cycle_nudge_warn USE carbon_cycle_mod, ONLY: nudge_carbon_cumul USE carbon_cycle_mod, ONLY: rco2_nudge_mean, rco2_nudge_target USE carbon_cycle_mod, ONLY: rco2_nudge_incr, rco2_nudge_flux USE carbon_cycle_mod, ONLY: nbcf_in, fields_in, cfname_in, cfunits_in USE carbon_cycle_mod, ONLY : nbcf_out, fields_out, yfields_out, cfname_out USE carbon_cycle_mod, ONLY: fco2_ocn_day, fco2_ff, fco2_bb USE carbon_cycle_mod, ONLY: fco2_land, fco2_ocean USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor, var_fco2_ocean_cor, fco2_ocean_cor USE carbon_cycle_mod, ONLY: carbon_cycle_ocean_flux_units USE carbon_cycle_mod, ONLY: read_fco2_land_cor, var_fco2_land_cor, fco2_land_cor USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor_2d, read_fco2_land_cor_2d USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_base, fco2_land_cor_2d_base USE carbon_cycle_mod, ONLY: co2_send USE carbon_cycle_mod, ONLY: fco2_land_nbp, fco2_land_nep, fco2_land_fLuc USE carbon_cycle_mod, ONLY: fco2_land_fwoodharvest, fco2_land_fHarvest USE carbon_cycle_mod, ONLY: carbon_cycle_cpl, carbon_cycle_tr USE carbon_cycle_mod, ONLY: carbon_cycle_rad, RCO2_glo, RCO2_tot USE carbon_cycle_mod, ONLY: carbon_cycle_conc_driven USE carbon_cycle_mod, ONLY: ocean_area_tot_glo, land_area_tot_glo USE carbon_cycle_mod, ONLY: ocean_area_tot, land_area_tot USE carbon_cycle_mod, ONLY: co2mass_diag ! >> PC USE mod_grid_phy_lmdz USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root ! << PC USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter, reduce_sum, reduce_max ! >> PC USE mod_phys_lmdz_omp_data, ONLY: is_omp_root USE phys_cal_mod USE phys_state_var_mod, ONLY: pctsrf USE indice_sol_mod, ONLY: nbsrf, is_ter, is_lic, is_oce, is_sic USE yomcst_mod_h USE clesphys_mod_h USE print_control_mod, ONLY: lunout IMPLICIT NONE !----------------------------------------------------------------------------- ! Arguments !----------------------------------------------------------------------------- REAL, INTENT(IN) :: pdtphys ! Physics time step [s] LOGICAL, INTENT(IN) :: debutphy ! TRUE at first physics call REAL, DIMENSION(klon), INTENT(IN) :: xlat, xlon, pphis REAL, DIMENSION(klon,klev), INTENT(IN) :: pphi, t_seri, pplay REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs REAL, DIMENSION(klon,nbtr), INTENT(INOUT) :: source ! Surface flux [U/m2/s] REAL, DIMENSION(klon,klev,nbtr), INTENT(INOUT) :: tr_seri ! Tracer conc [U/kgA] !----------------------------------------------------------------------------- ! Local variables !----------------------------------------------------------------------------- CHARACTER(LEN=20), PARAMETER :: modname = 'tracco2i' INTEGER :: it, k, i, nb REAL, DIMENSION(klon,klev) :: m_air ! Mass of air in each grid box [kg] ! Global arrays for gather REAL, DIMENSION(klon_glo,klev) :: co2_glo, m_air_glo REAL, DIMENSION(klon_glo,nbsrf) :: pctsrf_glo REAL, DIMENSION(klon_glo) :: pctsrf_ter_glo, pctsrf_oce_glo, pctsrf_sic_glo REAL, DIMENSION(klon_glo) :: cell_area_glo ! <> PC ! Local sums for MPI reduction REAL :: co2_mass_local, air_mass_local REAL :: co2_mass_global, air_mass_global REAL, PARAMETER :: secinday = 86400.0 ! << PC !--------------------------------------------------------------------------- ! Global test on the initial state of the CO2 field !--------------------------------------------------------------------------- REAL :: co2_absmax_local, co2_absmax_global LOGICAL :: co2_field_is_empty !--------------------------------------------------------------------------- ! Global-mean relaxation of the transported CO2 ! ! rco2_exact is the unrounded mass-weighted mean, kept apart from RCO2_glo on ! purpose. RCO2_glo carries a legacy rounding to 1.0e-8 kg/kg, near 0.0066 ! ppm. That rounding is harmless for the concentration handed to ORCHIDEE and ! PISCES, and ruinous for the relaxation !--------------------------------------------------------------------------- REAL :: rco2_exact, rco2_target, d_rco2_nudge REAL :: nudge_weight, nudge_flux_pgc, nudge_incr_ppm, nudge_carbon_pgc !--- Conversion from mixing ratio [kg CO2 / kg dry air] to dry mole fraction ! [ppm]. Written once here rather than repeated inline. REAL :: kgkg_to_ppm ! >> PC !=========================================================================== ! 0. SAFETY CHECK: Validate tracer index bounds ! ! When CO2 is not in tracer.def, id_CO2 may be invalid. ! We MUST check before any access to tr_seri or source with id_CO2. !=========================================================================== IF (first_call) THEN ! Check if id_CO2 is a valid tracer index IF (id_CO2 >= 1 .AND. id_CO2 <= nbtr) THEN co2_tracer_available = .TRUE. ELSE co2_tracer_available = .FALSE. END IF IF (is_omp_root .AND. is_mpi_root) THEN WRITE(lunout,*) '==============================================' WRITE(lunout,*) modname, ': SAFETY CHECK' WRITE(lunout,*) ' nbtr = ', nbtr WRITE(lunout,*) ' id_CO2 = ', id_CO2 WRITE(lunout,*) ' carbon_cycle_tr = ', carbon_cycle_tr WRITE(lunout,*) ' carbon_cycle_conc_driven = ', carbon_cycle_conc_driven WRITE(lunout,*) ' co2_tracer_available = ', co2_tracer_available IF (.NOT. co2_tracer_available) THEN WRITE(lunout,*) ' WARNING: CO2 tracer NOT in tracer list!' WRITE(lunout,*) ' Tracer operations will be SKIPPED.' WRITE(lunout,*) ' co2_send will use co2_ppm directly.' END IF WRITE(lunout,*) '==============================================' END IF first_call = .FALSE. END IF !=========================================================================== ! 1. INITIALISATION AT FIRST PHYSICS STEP !=========================================================================== IF (debutphy) THEN ! Initialise tracer field ONLY if CO2 tracer is available IF (co2_tracer_available) THEN IF (carbon_cycle_conc_driven) THEN ! Concentration-driven: force tracer to prescribed value tr_seri(:,:,id_CO2) = co2_ppm * 1.0e-6 * (RMCO2 / RMD) ELSE ! << PC !--- Emission-driven: initialise only if the field is ! empty, and take the test over the whole globe. ! co2_absmax_local = MAXVAL(ABS(tr_seri(:,:,id_CO2))) co2_field_is_empty = .FALSE. CALL reduce_max(co2_absmax_local, co2_absmax_global) !$OMP MASTER IF (is_mpi_root) THEN co2_field_is_empty = (co2_absmax_global < 1.0e-15) END IF !$OMP END MASTER CALL bcast(co2_field_is_empty) IF (co2_field_is_empty) THEN tr_seri(:,:,id_CO2) = co2_ppm0 * 1.0e-6 * (RMCO2 / RMD) IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) modname, ': CO2 field empty on the whole', & ' globe, initialised to co2_ppm0 =', co2_ppm0, ' ppm' END IF ELSE IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) modname, ': CO2 field inherited from the', & ' restart, global max =', & co2_absmax_global * 1.0e6 * RMD / RMCO2, ' ppm' END IF END IF ! >> PC END IF END IF ! << PC !--- The relaxation accumulator is NOT reset here. It is ! carried across restarts and its initial value is set once, by ! phyetat0, from the restart or from carbon_cycle_nudge_reset. Resetting ! it at every segment start would wipe the restored value. ! >> PC ! Check if fCO2_nbp is in coupling fields check_fCO2_nbp_in_cfname = .FALSE. DO nb = 1, nbcf_in IF (cfname_in(nb) == "fCO2_nbp") check_fCO2_nbp_in_cfname = .TRUE. END DO ! Gather surface fractions for area calculations CALL gather(pctsrf, pctsrf_glo) CALL gather(pctsrf(:,is_ter), pctsrf_ter_glo) CALL gather(pctsrf(:,is_oce), pctsrf_oce_glo) CALL gather(pctsrf(:,is_sic), pctsrf_sic_glo) CALL gather(cell_area(:), cell_area_glo) ! Compute total areas for flux corrections ! IMPORTANT: areas are computed whenever ANY correction mode is ! active (scalar OR 2D), not only in scalar mode. This ensures: ! (a) Mixed modes (e.g. scalar ocean + 2D land) work correctly ! (b) Diagnostic output is meaningful in all modes ! (c) Global integral verification of 2D fields is possible IF (read_fco2_ocean_cor .OR. read_fco2_ocean_cor_2d) THEN !$OMP MASTER IF (is_mpi_root .AND. is_omp_root) THEN ocean_area_tot = 0.0 DO i = 1, klon_glo ocean_area_tot = ocean_area_tot + & (pctsrf_oce_glo(i) + pctsrf_sic_glo(i)) * cell_area_glo(i) END DO END IF !$OMP END MASTER CALL bcast(ocean_area_tot) END IF IF (read_fco2_land_cor .OR. read_fco2_land_cor_2d) THEN !$OMP MASTER IF (is_mpi_root .AND. is_omp_root) THEN land_area_tot = 0.0 DO i = 1, klon_glo land_area_tot = land_area_tot + pctsrf_ter_glo(i) * cell_area_glo(i) END DO END IF !$OMP END MASTER CALL bcast(land_area_tot) END IF ! Diagnostic output: show correction mode and areas IF (is_omp_root .AND. is_mpi_root) THEN WRITE(lunout,*) modname, ': --- Flux correction configuration ---' WRITE(lunout,*) modname, ': read_fco2_ocean_cor = ', read_fco2_ocean_cor WRITE(lunout,*) modname, ': read_fco2_ocean_cor_2d = ', read_fco2_ocean_cor_2d WRITE(lunout,*) modname, ': read_fco2_land_cor = ', read_fco2_land_cor WRITE(lunout,*) modname, ': read_fco2_land_cor_2d = ', read_fco2_land_cor_2d WRITE(lunout,*) modname, ': ocean_area_tot =', ocean_area_tot, ' m2' WRITE(lunout,*) modname, ': land_area_tot =', land_area_tot, ' m2' IF (read_fco2_ocean_cor) THEN WRITE(lunout,*) modname, ': var_fco2_ocean_cor =', & var_fco2_ocean_cor, ' PgC/yr (scalar mode)' END IF IF (read_fco2_land_cor) THEN WRITE(lunout,*) modname, ': var_fco2_land_cor =', & var_fco2_land_cor, ' PgC/yr (scalar mode)' END IF END IF END IF ! debutphy !=========================================================================== ! 2. COMPUTE MASS OF AIR !=========================================================================== DO k = 1, klev m_air(:,k) = (paprs(:,k) - paprs(:,k+1)) / RG * cell_area(:) END DO !=========================================================================== ! 3. READ/UPDATE EMISSIONS !=========================================================================== CALL co2_emissions(debutphy) !=========================================================================== ! 3b. READ 2D FLUX CORRECTIONS (once at debutphy) ! Populates fco2_ocean_cor / fco2_land_cor ! Only active when read_fco2_ocean_cor_2d or read_fco2_land_cor_2d ! is set to .TRUE. in physiq.def. !=========================================================================== CALL co2_flux_corrections(debutphy) !=========================================================================== ! 4. PROCESS COUPLING FIELDS !=========================================================================== fco2_land(:) = 0.0 fco2_ocean(:) = 0.0 fco2_land_nbp(:) = 0.0 fco2_land_nep(:) = 0.0 fco2_land_fLuc(:) = 0.0 fco2_land_fwoodharvest(:) = 0.0 fco2_land_fHarvest(:) = 0.0 ! << PC ! ! Map coupling fields from ORCHIDEE/PISCES ! DO nb = 1, nbcf_in ! SELECT CASE (cfname_in(nb)) ! CASE ("fCO2_nep") ! fco2_land_nep(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter) ! CASE ("fCO2_fLuc") ! fco2_land_fLuc(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter) ! CASE ("fCO2_fwoodharvest") ! fco2_land_fwoodharvest(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter) ! CASE ("fCO2_fHarvest") ! fco2_land_fHarvest(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter) ! CASE ("fCO2_nbp") ! fco2_land_nbp(:) = fields_in(:,nb) * RMCO2/RMC * pctsrf(:,is_ter) ! CASE ("fCO2_fgco2") ! fco2_ocean(:) = -1.0 * fco2_ocn_day(:) * RMCO2/1.0e3 * & ! (pctsrf(:,is_oce) + pctsrf(:,is_sic)) ! END SELECT ! END DO ! Map coupling fields from ORCHIDEE and PISCES, honouring the declared unit DO nb = 1, nbcf_in SELECT CASE (cfname_in(nb)) CASE ("fCO2_nep") fco2_land_nep(:) = & orc_flux_to_kgco2_per_m2_s(fields_in(:,nb), cfunits_in(nb)) & * pctsrf(:,is_ter) CASE ("fCO2_fLuc") fco2_land_fLuc(:) = & orc_flux_to_kgco2_per_m2_s(fields_in(:,nb), cfunits_in(nb)) & * pctsrf(:,is_ter) ! << PC !--- ORCHIDEE sends this field as fCO2_fWoodharvest ! CASE ("fCO2_fwoodharvest") CASE ("fCO2_fWoodharvest") ! >> PC fco2_land_fwoodharvest(:) = & orc_flux_to_kgco2_per_m2_s(fields_in(:,nb), cfunits_in(nb)) & * pctsrf(:,is_ter) CASE ("fCO2_fHarvest") fco2_land_fHarvest(:) = & orc_flux_to_kgco2_per_m2_s(fields_in(:,nb), cfunits_in(nb)) & * pctsrf(:,is_ter) CASE ("fCO2_nbp") fco2_land_nbp(:) = & orc_flux_to_kgco2_per_m2_s(fields_in(:,nb), cfunits_in(nb)) & * pctsrf(:,is_ter) CASE ("fCO2_fgco2") !-- The ocean flux arrives through fco2_ocn_day, filled by cpl_mod. ! Its unit is the one declared for this field. The sign is negative: ! a positive ocean carbon flux is an uptake, a sink for the ! atmosphere. ! << PC !--- coherence barrier: LMDZ cannot read the NEMO namelist, so it ! checks the observable consequence: if the carbon coupling is requested ! yet the first ocean flux arrives at the coupler fill value, the two ! sides disagree and the run is stopped at once rather than left to ! drift IF (debutphy .AND. carbon_cycle_cpl) THEN IF (ANY(ABS(fco2_ocn_day(:)) > 1.0e-2)) THEN IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) modname, ': FATAL, ocean CO2 flux at fill value at first step' WRITE(lunout,*) modname, ': carbon_cycle_cpl is active but NEMO sends no CO2 flux' WRITE(lunout,*) modname, ': set sn_snd_co2 = coupled in the NEMO namelist' END IF CALL abort_physic(modname, 'ocean CO2 coupling requested but not sent', 1) END IF END IF ! >> PC ! << PC !-- Reject the coupler fill value before it enters the source. ! A correct air-sea carbon flux never exceeds a few 1e-6 in ! magnitude. When NEMO runs with sn_snd_co2 = none this field is ! never sent and arrives at its fill value, about -99999, which the ! mapping would turn into a constant source of order 1e5 ! kgCO2 m-2 s-1, the cause of the tracer blow-up IF (ANY(ABS(fco2_ocn_day(:)) > 1.0e-2)) THEN IF (is_mpi_root .AND. is_omp_root) WRITE(lunout,*) & 'tracco2i WARNING: ocean CO2 flux non-physical, ', & 'check sn_snd_co2 in the NEMO namelist' WHERE (ABS(fco2_ocn_day(:)) > 1.0e-2) fco2_ocn_day(:) = 0.0 END IF ! >> PC ! << PC ! fco2_ocean(:) = -1.0 & ! * orc_flux_to_kgco2_per_m2_s(fco2_ocn_day(:), cfunits_in(nb)) & ! * (pctsrf(:,is_oce) + pctsrf(:,is_sic)) ! >> PC ! << PC !--- Convert the ocean flux from the unit PISCES sends, moles of carbon ! per square metre per second, to the kilograms of CO2 per square ! metre per second the source term expects. The unit is taken from ! carbon_cycle_ocean_flux_units and not from the coupling channel, ! because the flux arrives on fco2_ocn_day while that channel stands ! empty fco2_ocean(:) = -1.0 & * orc_flux_to_kgco2_per_m2_s(fco2_ocn_day(:), carbon_cycle_ocean_flux_units) & * (pctsrf(:,is_oce) + pctsrf(:,is_sic)) ! >> PC ! << PC ! an unmatched coupling-field name is reported rather than dropped in silence CASE DEFAULT IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) modname, ': WARNING, received coupling field [', & TRIM(cfname_in(nb)), '] matches no known name and was ignored;', & ' check coupling_fields.def against the CASE names above' END IF ! >> PC END SELECT END DO ! << PC IF (is_mpi_root .AND. is_omp_root) THEN DO nb = 1, nbcf_in WRITE(lunout,*) 'RECV fields_in ', TRIM(cfname_in(nb)), & ' unit=[', TRIM(cfunits_in(nb)), '] min/max=', & MINVAL(fields_in(:,nb)), MAXVAL(fields_in(:,nb)) END DO WRITE(lunout,*) 'RECV fco2_ocn_day min/max=', & MINVAL(fco2_ocn_day), MAXVAL(fco2_ocn_day) END IF ! >> PC !=========================================================================== ! 5. APPLY FLUX CORRECTIONS ! ! The 2D base rate field is multiplied by the CURRENT ! pctsrf at EVERY time step, so the correction tracks sea-ice ! evolution dynamically. This makes the 2D path mathematically ! identical to the scalar path (both use runtime pctsrf). ! ! Use cases: ! - 2D alone: base_rate * pctsrf (from piControl correction file) ! - Scalar alone: scalar_rate * pctsrf (from physiq.def parameter) ! - Both: additive combination (fine-tuning) ! ! SIGN CONVENTION (source term, Section 7): ! source = fco2_ff + fco2_bb + fco2_land + fco2_ocean ! - fco2_ocean_cor - fco2_land_cor ! ! Positive correction REDUCES atmospheric CO2 (net sink). ! Negative correction INCREASES atmospheric CO2 (net source). !=========================================================================== ! --- Ocean correction --- ! Always start from zero, then add active components fco2_ocean_cor(:) = 0.0 ! (a) 2D base rate * current pctsrf IF (read_fco2_ocean_cor_2d) THEN fco2_ocean_cor(:) = fco2_ocean_cor_2d_base(:) * & (pctsrf(:,is_oce) + pctsrf(:,is_sic)) END IF ! (b) Scalar uniform correction (additive if 2D also active) IF (read_fco2_ocean_cor .AND. ocean_area_tot > 0.0) THEN fco2_ocean_cor(:) = fco2_ocean_cor(:) + & var_fco2_ocean_cor * 1.0e12 * RMCO2/RMC / & ocean_area_tot / secinday / REAL(year_len) * & (pctsrf(:,is_oce) + pctsrf(:,is_sic)) END IF ! --- Land correction --- fco2_land_cor(:) = 0.0 ! (a) 2D base rate * current pctsrf IF (read_fco2_land_cor_2d) THEN fco2_land_cor(:) = fco2_land_cor_2d_base(:) * pctsrf(:,is_ter) END IF ! (b) Scalar uniform correction (additive if 2D also active) IF (read_fco2_land_cor .AND. land_area_tot > 0.0) THEN fco2_land_cor(:) = fco2_land_cor(:) + & var_fco2_land_cor * 1.0e12 * RMCO2/RMC / & land_area_tot / secinday / REAL(year_len) * & pctsrf(:,is_ter) END IF !=========================================================================== ! 6. AGGREGATE LAND FLUXES !=========================================================================== IF (check_fCO2_nbp_in_cfname) THEN fco2_land(:) = fco2_land_nbp(:) ELSE fco2_land(:) = fco2_land_nep(:) + fco2_land_fLuc(:) + & fco2_land_fwoodharvest(:) + fco2_land_fHarvest(:) END IF ! << PC CALL gather(fco2_ff, fco2_ff_glo) CALL gather(fco2_bb, fco2_bb_glo) CALL gather(fco2_land, fco2_land_glo) CALL gather(fco2_ocean, fco2_ocean_glo) IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) 'tracco2i DIAG: fco2_ff min/max =', MINVAL(fco2_ff_glo), MAXVAL(fco2_ff_glo) WRITE(lunout,*) 'tracco2i DIAG: fco2_bb min/max =', MINVAL(fco2_bb_glo), MAXVAL(fco2_bb_glo) WRITE(lunout,*) 'tracco2i DIAG: fco2_land min/max =', MINVAL(fco2_land_glo), MAXVAL(fco2_land_glo) WRITE(lunout,*) 'tracco2i DIAG: fco2_ocean min/max =', MINVAL(fco2_ocean_glo), MAXVAL(fco2_ocean_glo) END IF !! << PC restart probe ! !--- Global checksums of the ocean flux as received (fco2_ocn_day) and as ! ! used (fco2_ocean), and of the land flux, at every step: a serial sum ! ! over the gathered global vector on the MPI root, printed with sixteen ! ! significant digits. A difference anywhere on the globe shows here, ! ! where the root-subdomain min/max of the RECV lines cannot see it. ! ! *Diagnostic only, to be removed once the restart test is closed* ! CALL gather(fco2_ocn_day, field_in_glo) ! IF (is_mpi_root .AND. is_omp_root) THEN ! WRITE(lunout,'(A,4(1X,ES24.16))') 'RSTPROBE fco2_ocn_day sum sumsq min max', & ! SUM(field_in_glo), SUM(field_in_glo * field_in_glo), & ! MINVAL(field_in_glo), MAXVAL(field_in_glo) ! WRITE(lunout,'(A,4(1X,ES24.16))') 'RSTPROBE fco2_ocean sum sumsq min max', & ! SUM(fco2_ocean_glo), SUM(fco2_ocean_glo * fco2_ocean_glo), & ! MINVAL(fco2_ocean_glo), MAXVAL(fco2_ocean_glo) ! WRITE(lunout,'(A,4(1X,ES24.16))') 'RSTPROBE fco2_land sum sumsq min max', & ! SUM(fco2_land_glo), SUM(fco2_land_glo * fco2_land_glo), & ! MINVAL(fco2_land_glo), MAXVAL(fco2_land_glo) ! END IF ! >> PC restart probe DO nb = 1, nbcf_in CALL gather(fields_in(:,nb), field_in_glo) IF (is_mpi_root .AND. is_omp_root) & WRITE(lunout,*) 'tracco2i DIAG: field ', TRIM(cfname_in(nb)), & ' units=[', TRIM(cfunits_in(nb)), '] min/max=', MINVAL(field_in_glo), MAXVAL(field_in_glo) END DO ! >> PC !=========================================================================== ! 7. BUILD FINAL SOURCE TERM !=========================================================================== IF (co2_tracer_available) THEN ! Tracer exists - can write to source(:,id_CO2) IF (carbon_cycle_conc_driven) THEN !--------------------------------------------------------------------- ! CONCENTRATION-DRIVEN MODE ! Source term is zero (no flux feedback to atmosphere) ! Fluxes computed above are for diagnostic outputs only !--------------------------------------------------------------------- source(:,id_CO2) = 0.0 ELSE !--------------------------------------------------------------------- ! EMISSION-DRIVEN MODE ! Sum all flux components !--------------------------------------------------------------------- source(:,id_CO2) = fco2_ff(:) + fco2_bb(:) + fco2_land(:) + fco2_ocean(:) & - fco2_ocean_cor(:) - fco2_land_cor(:) ! << PC diagnostic source, a retirer apres usage IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) 'SRC pctsrf_ter min/max=', & MINVAL(pctsrf(:,is_ter)), MAXVAL(pctsrf(:,is_ter)) WRITE(lunout,*) 'SRC fco2_ff min/max=', & MINVAL(fco2_ff), MAXVAL(fco2_ff) WRITE(lunout,*) 'SRC fco2_bb min/max=', & MINVAL(fco2_bb), MAXVAL(fco2_bb) WRITE(lunout,*) 'SRC fco2_land min/max=', & MINVAL(fco2_land), MAXVAL(fco2_land) WRITE(lunout,*) 'SRC fco2_ocean min/max=', & MINVAL(fco2_ocean), MAXVAL(fco2_ocean) WRITE(lunout,*) 'SRC source_CO2 min/max=', & MINVAL(source(:,id_CO2)), MAXVAL(source(:,id_CO2)) END IF ! >> PC END IF END IF ! If co2_tracer_available = .FALSE., we do NOT touch source(:,id_CO2) !=========================================================================== ! 8. COMPUTE GLOBAL MEAN CO2 (daily update) !=========================================================================== IF (debutphy .OR. day_cur /= day_pre) THEN IF (carbon_cycle_conc_driven) THEN !--------------------------------------------------------------------- ! Concentration-driven: prescribe from co2_ppm !--------------------------------------------------------------------- RCO2_glo = co2_ppm * 1.0e-6 * (RMCO2/RMD) RCO2_glo = REAL(INT(RCO2_glo * 1.0e8)) / 1.0e8 ! Force tracer to prescribed value (ONLY if available) IF (co2_tracer_available) THEN tr_seri(:,:,id_CO2) = RCO2_glo END IF IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) modname, ': Concentration-driven CO2 =', co2_ppm, ' ppm' END IF ELSE IF (co2_tracer_available) THEN !--------------------------------------------------------------------- ! Emission-driven with tracer available: compute from tracer field !--------------------------------------------------------------------- ! << PC ! co2_mass_local = SUM(tr_seri(:,:,id_CO2) * m_air(:,:)) ! air_mass_local = SUM(m_air(:,:)) ! ! co2_mass_global = 0.0 ! air_mass_global = 0.0 ! CALL reduce_sum(co2_mass_local, co2_mass_global) ! CALL reduce_sum(air_mass_local, air_mass_global) ! !!$OMP BARRIER ! IF (air_mass_global > 0.0) THEN ! RCO2_glo = co2_mass_global / air_mass_global ! ELSE ! RCO2_glo = co2_ppm0 * 1.0e-6 * (RMCO2/RMD) ! END IF ! RCO2_glo = REAL(INT(RCO2_glo * 1.0e8)) / 1.0e8 !!$OMP BARRIER ! ! IF (is_mpi_root .AND. is_omp_root) THEN ! WRITE(lunout,*) modname, ': global CO2 =', RCO2_glo*1.0e6*RMD/RMCO2, ' ppm' ! END IF ! ! ! If no transport flag, homogenise tracer field ! IF (.NOT. carbon_cycle_tr) THEN ! tr_seri(:,:,id_CO2) = RCO2_glo ! END IF ! ! ELSE ! !--------------------------------------------------------------------- ! ! Tracer not available: use prescribed value ! !--------------------------------------------------------------------- ! RCO2_glo = co2_ppm * 1.0e-6 * (RMCO2/RMD) ! END IF !!!--- !--- Atmospheric CO2 burden divided by the mass of the atmosphere. ! Weighting by the air mass of each layer is what makes the result ! a conserved quantity. An arithmetic mean over grid points would ! give a thin polar column the same say as a deep tropical one, ! and would drift as soon as the surface pressure field evolves. ! This is the concentration that ORCHIDEE and PISCES are told to ! breathe, and, when carbon_cycle_rad is on, the one the radiation ! sees. co2_mass_local = SUM(tr_seri(:,:,id_CO2) * m_air(:,:)) air_mass_local = SUM(m_air(:,:)) co2_mass_global = 0.0 air_mass_global = 0.0 CALL reduce_sum(co2_mass_local, co2_mass_global) CALL reduce_sum(air_mass_local, air_mass_global) ! << PC !--- Conversion factor, mixing ratio to dry mole fraction kgkg_to_ppm = 1.0e6 * RMD / RMCO2 rco2_exact = 0.0 rco2_target = 0.0 d_rco2_nudge = 0.0 nudge_flux_pgc = 0.0 ! >> PC !$OMP MASTER IF (is_mpi_root) THEN IF (air_mass_global > 0.0) THEN rco2_exact = co2_mass_global / air_mass_global ELSE !--- An atmosphere of zero mass means the pressure field has ! been lost. Fall back on the prescribed pre-industrial ! value, and say so, rather than carry a meaningless ratio. rco2_exact = co2_ppm0 * 1.0e-6 * (RMCO2/RMD) WRITE(lunout,*) modname, ': WARNING, null global air mass,', & ' global CO2 fallen back on co2_ppm0' END IF ! << PC !--- Pull the global mean towards its target ! ! The increment is uniform in mixing ratio. Adding a constant c ! to q shifts the mass-weighted mean, the ratio of the sum of q ! times m to the sum of m, by c times the sum of m over the sum ! of m, which is c exactly. No iteration and no residual. The ! same constant cancels in any difference between two points, so ! the inter-hemispheric gradient, the seasonal cycle and the ! vertical structure all pass through untouched. What is held ! here is the single degree of freedom that the surface fluxes ! cannot set on their own ! ! Everything below is formed from rco2_exact. The rounding is ! applied afterwards, to RCO2_glo alone. Computing the increment ! from the rounded mean would surrender up to half a rounding ! step every day, that is about 2.6 PgC per year on the diagnosed ! flux, and would destroy the only number this device exists to ! produce IF (carbon_cycle_nudge) THEN IF (carbon_cycle_nudge_target > 0.0) THEN rco2_target = carbon_cycle_nudge_target * 1.0e-6 * (RMCO2/RMD) ELSE rco2_target = co2_ppm * 1.0e-6 * (RMCO2/RMD) END IF !--- The update is daily, so the weight is one day over tau, ! capped at unity. A tau of one day resets the mean exactly. nudge_weight = 1.0 / MAX(carbon_cycle_nudge_tau, 1.0) d_rco2_nudge = nudge_weight * (rco2_target - rco2_exact) !--- Carbon handed to the atmosphere by the relaxation, put in the ! units of the carbon budget. This number measures how far the ! coupled system sits from the target = a spin-up has converged ! when it settles near zero. The daily amount is annualised by ! year_len, valid because this block runs once per day. nudge_carbon_pgc = d_rco2_nudge * air_mass_global & * (RMC/RMCO2) * 1.0e-12 nudge_flux_pgc = nudge_carbon_pgc * REAL(year_len) END IF ! >> PC !--- Round to 1.0e-8 kg/kg, near 0.0066 ppm, so that the diagnosed ! concentration is one number rather than a family of values ! differing in their last bits. Rounding to nearest, not ! truncation: truncation would place the diagnostic systematically ! below the true mean. RCO2_glo carries the post-nudge mean, so it ! stays consistent with the tracer that is shifted below. RCO2_glo = rco2_exact + d_rco2_nudge IF (RCO2_glo > 0.0 .AND. RCO2_glo < 1.0) THEN RCO2_glo = REAL(NINT(RCO2_glo * 1.0e8)) / 1.0e8 END IF END IF !! << PC ! !--- Total atmospheric CO2 mass ! ! * Emission-driven: the mass the tracer carries ! ! * Concentration-driven, or no tracer: the mass implied by the prescribed ! ! global mixing ratio (same molar masses as the transport) ! IF (carbon_cycle_conc_driven .OR. .NOT. co2_tracer_available) THEN ! co2mass_diag = co2_ppm * 1.0e-6 * (RMCO2 / RMD) * air_mass_global ! ELSE ! co2mass_diag = co2_mass_global ! END IF !! >> PC !$OMP END MASTER CALL bcast(RCO2_glo) ! << PC CALL bcast(rco2_exact) CALL bcast(rco2_target) CALL bcast(d_rco2_nudge) CALL bcast(nudge_flux_pgc) CALL bcast(co2mass_diag) !--- Every thread applies the same increment to its own points. The ! uniform offset shifts the burden without touching any gradient. IF (carbon_cycle_nudge) THEN tr_seri(:,:,id_CO2) = tr_seri(:,:,id_CO2) + d_rco2_nudge END IF ! >> PC IF (is_mpi_root .AND. is_omp_root) THEN WRITE(lunout,*) modname, ': global CO2 =', RCO2_glo*1.0e6*RMD/RMCO2, ' ppm' ! << PC IF (carbon_cycle_nudge) THEN nudge_incr_ppm = d_rco2_nudge * kgkg_to_ppm nudge_carbon_cumul = nudge_carbon_cumul & + nudge_flux_pgc / REAL(year_len) !--- One line per day, fixed fields, single writer, meant to be ! read back by a script rather than by eye. The tag is there so ! that a grep gives a table and nothing else. Fields: year, day ! of month, month, exact global mean [ppm], target [ppm], ! applied increment [ppm], annualised flux [PgC/yr], carbon ! injected since the last reset [PgC]. WRITE(lunout,'(A,1X,I6,1X,I3,1X,I3,5(1X,ES16.8))') & 'NUDGE_DIAG', year_cur, day_cur, mth_cur, & rco2_exact * kgkg_to_ppm, & rco2_target * kgkg_to_ppm, & nudge_incr_ppm, & nudge_flux_pgc, & nudge_carbon_cumul !--- A large increment in a single day means the tracer was not ! where it was believed to be. Say so loudly, and carry on: ! aborting halfway through a spin-up costs more than a noisy ! log, and the run remains diagnosable afterwards. IF (ABS(nudge_incr_ppm) > carbon_cycle_nudge_warn) THEN WRITE(lunout,*) modname, ': WARNING, relaxation increment', & nudge_incr_ppm, ' ppm exceeds', carbon_cycle_nudge_warn, & ' ppm. Check the initial state of the CO2 tracer.' END IF END IF ! >> PC END IF ! << PC !--- Publish the relaxation diagnostics for the history rco2_nudge_mean = rco2_exact * kgkg_to_ppm rco2_nudge_target = rco2_target * kgkg_to_ppm rco2_nudge_incr = d_rco2_nudge * kgkg_to_ppm rco2_nudge_flux = nudge_flux_pgc CALL bcast(nudge_carbon_cumul) ! >> PC ! << PC !--- Give co2_ppm global-mean CO2 the run handed to the radiation (RCO2_glo) IF (carbon_cycle_tr .AND. carbon_cycle_rad & .AND. .NOT. carbon_cycle_conc_driven & .AND. .NOT. carbon_cycle_nudge) THEN co2_ppm = RCO2_glo * kgkg_to_ppm END IF ! >> PC !--- With transport switched off the atmosphere is a single well-mixed ! box: every grid point is reset to the global mean, so that surface ! fluxes change the burden without building any gradient. This branch ! is unreachable while the relaxation is active, because conf_phys ! refuses carbon_cycle_nudge without carbon_cycle_tr. IF (.NOT. carbon_cycle_tr) THEN tr_seri(:,:,id_CO2) = RCO2_glo END IF ELSE !--------------------------------------------------------------------- ! No CO2 tracer: the concentration is the prescribed one !--------------------------------------------------------------------- RCO2_glo = co2_ppm * 1.0e-6 * (RMCO2/RMD) END IF ! >> PC ! << PC !--- Total atmospheric CO2 mass: ! with a transported tracer the mass is the true three-dimensional integral ! co2_mass_global. Otherwise the atmosphere carries the prescribed ! global mean, and the mass is that mean weighted by the total air ! mass IF (carbon_cycle_conc_driven .OR. .NOT. co2_tracer_available) THEN air_mass_local = SUM(m_air(:,:)) air_mass_global = 0.0 CALL reduce_sum(air_mass_local, air_mass_global) END IF !$OMP MASTER IF (is_mpi_root) THEN IF (carbon_cycle_conc_driven .OR. .NOT. co2_tracer_available) THEN co2mass_diag = co2_ppm * 1.0e-6 * (RMCO2 / RMD) * air_mass_global ELSE co2mass_diag = co2_mass_global END IF END IF !$OMP END MASTER CALL bcast(co2mass_diag) ! >> PC day_pre = day_cur END IF !=========================================================================== ! 9. UPDATE CO2 TO SEND TO SURFACE MODELS ! ! In concentration-driven mode OR when tracer not available: ! ==> send co2_ppm (prescribed) ! In emission-driven mode with tracer: ! ==> send computed atmospheric value !=========================================================================== IF (ALLOCATED(co2_send)) THEN IF (carbon_cycle_conc_driven .OR. .NOT. co2_tracer_available) THEN co2_send(:) = co2_ppm ELSE co2_send(:) = tr_seri(:,1,id_CO2) * 1.0e6 * RMD / RMCO2 END IF END IF !! << PC restart probe ! !--- Global checksum of the concentration handed to the surface models, ! ! same convention as the flux probes above ! IF (ALLOCATED(co2_send)) THEN ! CALL gather(co2_send, field_in_glo) ! IF (is_mpi_root .AND. is_omp_root) THEN ! WRITE(lunout,'(A,4(1X,ES24.16))') 'RSTPROBE co2_send sum sumsq min max', & ! SUM(field_in_glo), SUM(field_in_glo * field_in_glo), & ! MINVAL(field_in_glo), MAXVAL(field_in_glo) ! END IF ! END IF !! >> PC restart probe IF (is_mpi_root .AND. is_omp_root) THEN IF (co2_tracer_available) THEN WRITE(lunout,*) modname, ': tr_seri L1 min/max (ppm) =', & MINVAL(tr_seri(:,1,id_CO2)*1.0e6*RMD/RMCO2), & MAXVAL(tr_seri(:,1,id_CO2)*1.0e6*RMD/RMCO2) END IF IF (ALLOCATED(co2_send)) THEN WRITE(lunout,*) modname, ': co2_send min/max =', & MINVAL(co2_send), MAXVAL(co2_send) END IF END IF END SUBROUTINE tracco2i !=============================================================================== ! FUNCTION: to_lower !=============================================================================== PURE FUNCTION to_lower(s) RESULT(t) IMPLICIT NONE CHARACTER(*), INTENT(IN) :: s CHARACTER(LEN(s)) :: t INTEGER :: i, ia DO i = 1, LEN(s) ia = IACHAR(s(i:i)) IF (ia >= IACHAR('A') .AND. ia <= IACHAR('Z')) THEN t(i:i) = ACHAR(ia + 32) ELSE t(i:i) = s(i:i) END IF END DO END FUNCTION to_lower ! << PC !=========================================================================== ! FUNCTION: orc_flux_to_kgco2_per_m2_s ! ! Purpose: ! Convert a carbon flux from the unit declared in coupling_fields.def to ! kgCO2 per m2 per second. Mass and time conversion only. The caller ! applies the sign convention and the surface fraction. A unit not in the ! table stops the run, so a silent mis-scaling can no longer happen. ! ! Author: ! 2026: Patricia Cadule (IPSL) !=========================================================================== FUNCTION orc_flux_to_kgco2_per_m2_s(flux_in, units_in) RESULT(flux_out) USE yomcst_mod_h, ONLY: RMCO2, RMC USE print_control_mod, ONLY: lunout IMPLICIT NONE REAL, INTENT(IN) :: flux_in(:) CHARACTER(LEN=*), INTENT(IN) :: units_in REAL :: flux_out(SIZE(flux_in)) REAL :: factor CHARACTER(LEN=64) :: u REAL, PARAMETER :: secinday = 86400.0 CHARACTER(LEN=24), PARAMETER :: modname = 'orc_flux_to_kgco2' u = to_lower(ADJUSTL(TRIM(units_in))) SELECT CASE (u) CASE ('gc m-2 day-1', 'gc/m2/day', 'gc m-2 d-1', 'gc.m-2.day-1', 'gc.m-2.d-1') factor = 1.0e-3 / secinday * (RMCO2/RMC) ! gC/m2/day -> kgCO2/m2/s CASE ('gc m-2 s-1', 'gc/m2/s', 'gc.m-2.s-1') factor = 1.0e-3 * (RMCO2/RMC) ! gC/m2/s -> kgCO2/m2/s CASE ('kgc m-2 s-1', 'kgc/m2/s', 'kgc.m-2.s-1') factor = (RMCO2/RMC) ! kgC/m2/s -> kgCO2/m2/s CASE ('molc m-2 s-1', 'molc/m2/s', 'molc.m-2.s-1') factor = (RMCO2/1000.0) ! molC/m2/s -> kgCO2/m2/s CASE DEFAULT WRITE(lunout,*) modname, ': unsupported coupling unit [', & TRIM(units_in), ']' CALL abort_physic(modname, & 'unsupported carbon flux unit: '//TRIM(units_in), 1) END SELECT flux_out(:) = flux_in(:) * factor END FUNCTION orc_flux_to_kgco2_per_m2_s ! >> PC !=============================================================================== ! SUBROUTINE: ocean_flux_to_atm_source !=============================================================================== SUBROUTINE ocean_flux_to_atm_source(flux_in, units_in, frac_ocean, flux_out) USE yomcst_mod_h, ONLY: RMCO2, RMC USE clesphys_mod_h IMPLICIT NONE REAL, INTENT(IN) :: flux_in(:) CHARACTER(*), INTENT(IN) :: units_in REAL, INTENT(IN) :: frac_ocean(:) REAL, INTENT(OUT) :: flux_out(:) REAL :: factor CHARACTER(LEN=64) :: u REAL, PARAMETER :: secinday = 86400.0 CHARACTER(LEN=20), PARAMETER :: modname = 'ocean_flux_conv' u = to_lower(ADJUSTL(TRIM(units_in))) factor = -1.0 SELECT CASE (u) CASE ('gc m-2 day-1', 'gc/m2/day', 'gc m-2 d-1', 'gc.m-2.day-1', 'gc.m-2.d-1') flux_out(:) = factor * flux_in(:) * 1.0e-3 / secinday * (RMCO2/RMC) * frac_ocean(:) CASE ('gc m-2 s-1', 'gc/m2/s', 'gc.m-2.s-1') flux_out(:) = factor * flux_in(:) * 1.0e-3 * (RMCO2/RMC) * frac_ocean(:) CASE ('kgc m-2 s-1', 'kgc/m2/s', 'kgc.m-2.s-1') flux_out(:) = factor * flux_in(:) * (RMCO2/RMC) * frac_ocean(:) CASE ('molc.m-2.s-1', 'molc/m2/s', 'molc m-2 s-1') flux_out(:) = factor * flux_in(:) * (RMCO2 / 1000.0) * frac_ocean(:) CASE DEFAULT CALL abort_physic(modname, 'Unsupported ocean flux units: '//TRIM(units_in), 1) END SELECT END SUBROUTINE ocean_flux_to_atm_source !=============================================================================== ! "SZ98" ROUTINES ! ! These routines implement the Sheng & Zwiers (1998) interpolation method ! for converting monthly mean emissions to smooth daily values whilst ! exactly preserving the monthly means. ! ! Reference: ! Sheng, J., & Zwiers, F. (1998). An improved scheme for time-dependent ! boundary conditions in atmospheric general circulation models. ! Climate Dynamics, 14, 609-613. !=============================================================================== !=============================================================================== ! SUBROUTINE: co2_emissions ! ! Purpose: ! Read CO2 fossil-fuel (and optionally biomass-burning) emissions from ! NetCDF files and provide daily flux values to the atmospheric model. ! ! Authors: ! 2020: Olivier Boucher (IPSL/LMD) - initial ! 2026: Patricia Cadule (IPSL) - daily mode + SZ98 interpolation ! ! ! Operating Modes: ! (A) read_daily_co2ff = .TRUE. ! Input file has time = year_len (360/365/366) daily values. ! Direct indexing by day-of-year. ! ! (B) read_daily_co2ff = .FALSE. ! Input file has time = 12 monthly mean values. ! SZ98 interpolation provides smooth daily values that preserve ! the monthly means exactly. !! ! SZ98 Interpolation: ! Sheng & Zwiers (1998) method ensures monthly means are preserved ! when interpolating to daily values. The tridiagonal system is solved ! ONCE per year. ! Sheng, J. & Zwiers, F.W. (1998). An improved scheme for time-dependent ! boundary conditions in atmospheric general circulation models. ! Climate Dynamics, 14, 609-613. ! DOI: 10.1007/s003820050244 ! ! Units: ! flx_co2ff/bb : kg CO2 m-2 s-1 ! fco2_ff/bb : kg CO2 m-2 s-1 (output to source term) ! !=============================================================================== SUBROUTINE co2_emissions(debutphy) USE dimphy USE geometry_mod, ONLY: cell_area USE mod_grid_phy_lmdz USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root USE mod_phys_lmdz_omp_data, ONLY: is_omp_root, omp_rank USE mod_phys_lmdz_para, ONLY: scatter USE print_control_mod, ONLY: lunout USE phys_cal_mod, ONLY: year_len, day_cur, mth_cur USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open USE netcdf, ONLY: nf90_get_var, nf90_nowrite, nf90_noerr, nf90_strerror, & nf90_inq_dimid, nf90_inquire_dimension ! Import variables from module USE carbon_cycle_mod, ONLY: fco2_ff, fco2_bb USE carbon_cycle_mod, ONLY: read_daily_co2ff USE clesphys_mod_h USE yomcst_mod_h IMPLICIT NONE !----------------------------------------------------------------------------- ! Arguments !----------------------------------------------------------------------------- LOGICAL, INTENT(IN) :: debutphy ! .TRUE. at first physics call !----------------------------------------------------------------------------- ! ! In monthly mode: (klon, 12) - 12 monthly means ! In daily mode: (klon, year_len) - one value per day ! ! These are allocated once at debutphy and persist across time steps. !----------------------------------------------------------------------------- REAL, ALLOCATABLE, SAVE :: flx_co2ff(:,:) ! Fossil-fuel CO2 [kgCO2/m2/s] REAL, ALLOCATABLE, SAVE :: flx_co2bb(:,:) ! Biomass-burning CO2 [kgCO2/m2/s] !$OMP THREADPRIVATE(flx_co2ff, flx_co2bb) !-- SZ98 interpolation targets (monthly mode only, THREADPRIVATE) REAL, ALLOCATABLE, SAVE :: target_ff(:,:) ! (klon, 12) mid-month targets REAL, ALLOCATABLE, SAVE :: target_bb(:,:) ! (klon, 12) mid-month targets !$OMP THREADPRIVATE(target_ff, target_bb) !-- Calendar information for SZ98 interpolation INTEGER, SAVE :: days_in_mth(12) ! Days per month for current year REAL, SAVE :: month_midday(12) ! Day-of-year of month mid-points !$OMP THREADPRIVATE(days_in_mth, month_midday) !-- Daily update guard: avoid recalculating every time step INTEGER, SAVE :: day_pre_emis = -1 !$OMP THREADPRIVATE(day_pre_emis) !----------------------------------------------------------------------------- ! Local variables for NetCDF reading !----------------------------------------------------------------------------- REAL, ALLOCATABLE :: flx_co2ff_glo(:,:) ! Global field (MPI root only) REAL, ALLOCATABLE :: flx_co2bb_glo(:,:) ! Global field (MPI root only) INTEGER :: ncid_in, varid, ncerr INTEGER :: n_glo, n_time_file, m, n_times INTEGER :: current_year_len, doy REAL, ALLOCATABLE :: vector_coord(:), time_coord(:) !-- Control flags LOGICAL, PARAMETER :: readco2ff = .TRUE. LOGICAL, PARAMETER :: readco2bb = .FALSE. CHARACTER(LEN=30), PARAMETER :: modname = 'co2_emissions' CHARACTER(LEN=80) :: abort_message !============================================================================= ! 1. INITIALISATION !============================================================================= IF (debutphy) THEN current_year_len = INT(year_len) !-- Determine number of time steps in the file IF (read_daily_co2ff) THEN n_times = current_year_len ! 360, 365, or 366 ELSE n_times = 12 ! Monthly END IF !------------------------------------------------------------------ ! A. Allocate arrays !------------------------------------------------------------------ ALLOCATE(flx_co2ff(klon, n_times)) ALLOCATE(flx_co2bb(klon, n_times)) flx_co2ff(:,:) = 0.0 flx_co2bb(:,:) = 0.0 IF (.NOT. read_daily_co2ff) THEN ALLOCATE(target_ff(klon, 12)) ALLOCATE(target_bb(klon, 12)) target_ff(:,:) = 0.0 target_bb(:,:) = 0.0 END IF !------------------------------------------------------------------ ! B. NetCDF reading !------------------------------------------------------------------ !$OMP MASTER IF (is_mpi_root) THEN !-- Allocate global buffer on MPI root IF (.NOT. ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(klon_glo, n_times)) IF (.NOT. ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(klon_glo, n_times)) flx_co2ff_glo(:,:) = 0.0 flx_co2bb_glo(:,:) = 0.0 !--------------------------------------------------------------- ! B.1 Read fossil-fuel emissions !--------------------------------------------------------------- IF (readco2ff) THEN CALL nf95_open('sflx_lmdz_co2_ff.nc', nf90_nowrite, ncid_in) !-- Validate spatial dimension CALL nf95_inq_varid(ncid_in, 'vector', varid) CALL nf95_gw_var(ncid_in, varid, vector_coord) n_glo = SIZE(vector_coord) DEALLOCATE(vector_coord) IF (n_glo /= klon_glo) THEN abort_message = 'sflx_lmdz_co2_ff: vector /= klon_glo' CALL abort_physic(modname, abort_message, 1) END IF !-- Validate time dimension CALL nf95_inq_varid(ncid_in, 'time', varid) CALL nf95_gw_var(ncid_in, varid, time_coord) n_time_file = SIZE(time_coord) DEALLOCATE(time_coord) IF (n_time_file /= n_times) THEN WRITE(abort_message, '(A,I4,A,I4)') & 'sflx_lmdz_co2_ff: time=', n_time_file, ' exp=', n_times CALL abort_physic(modname, abort_message, 1) END IF !-- Read emission field directly into (klon_glo, n_times) array CALL nf95_inq_varid(ncid_in, 'flx_co2', varid) ncerr = nf90_get_var(ncid_in, varid, flx_co2ff_glo) CALL nf95_close(ncid_in) !-- Reorder from file convention to LMDZ internal ordering ! ! Longitude convention: ! Determined by LMDZ grid definition file used in preprocessing. ! If grid uses [0,360]: shift_lon = 0 ! If grid uses [-180,+180]: shift_lon = nbp_lon / 2 (default) ! ! If a different grid definition is used (e.g. xfirst=0), this value ! must be adjusted accordingly: shift_lon = 0 for xfirst=0. ! CSHIFT by nbp_lon/2 = 72 positions rotates each band by 180 degrees, ! aligning the Python [0,360] convention with the LMDZ [-180,+180] grid. !--------------------------------------------------------------- DO m = 1, n_times CALL reorder_global_1d(flx_co2ff_glo(:,m), nbp_lon, nbp_lat, & flip_lat=.FALSE., shift_lon=nbp_lon/2) END DO WRITE(lunout,*) modname, ': Reorder checked: flip_lat=F, shift_lon=', nbp_lon/2 WRITE(*,*) modname, ': Read FF emissions OK, n_times=', n_times IF (read_daily_co2ff) THEN WRITE(*,*) ' Mode: DAILY (direct indexing by day-of-year)' ELSE WRITE(*,*) ' Mode: MONTHLY (SZ98 interpolation to daily)' END IF ELSE flx_co2ff_glo(:,:) = 0.0 END IF !--------------------------------------------------------------- ! B.2 Read biomass-burning emissions (disabled by default) !--------------------------------------------------------------- IF (readco2bb) THEN CALL nf95_open('sflx_lmdz_co2_bb.nc', nf90_nowrite, ncid_in) CALL nf95_inq_varid(ncid_in, 'vector', varid) CALL nf95_gw_var(ncid_in, varid, vector_coord) n_glo = SIZE(vector_coord) DEALLOCATE(vector_coord) IF (n_glo /= klon_glo) THEN abort_message = 'sflx_lmdz_co2_bb: vector /= klon_glo' CALL abort_physic(modname, abort_message, 1) END IF CALL nf95_inq_varid(ncid_in, 'time', varid) CALL nf95_gw_var(ncid_in, varid, time_coord) n_time_file = SIZE(time_coord) DEALLOCATE(time_coord) IF (n_time_file /= n_times) THEN abort_message = 'sflx_lmdz_co2_bb: time mismatch' CALL abort_physic(modname, abort_message, 1) END IF CALL nf95_inq_varid(ncid_in, 'flx_co2', varid) ncerr = nf90_get_var(ncid_in, varid, flx_co2bb_glo) CALL nf95_close(ncid_in) DO m = 1, n_times CALL reorder_global_1d(flx_co2ff_glo(:,m), nbp_lon, nbp_lat, & flip_lat=.FALSE., shift_lon=nbp_lon/2) END DO ELSE flx_co2bb_glo(:,:) = 0.0 END IF END IF ! is_mpi_root !$OMP END MASTER !------------------------------------------------------------------ ! C. Distribute global field !------------------------------------------------------------------ !-- Dummy allocation for non-MPI-root ranks IF (.NOT. ALLOCATED(flx_co2ff_glo)) ALLOCATE(flx_co2ff_glo(0, 0)) IF (.NOT. ALLOCATED(flx_co2bb_glo)) ALLOCATE(flx_co2bb_glo(0, 0)) !-- Scatter: distributes first dimension (klon_glo ==> klon) CALL scatter(flx_co2ff_glo, flx_co2ff) CALL scatter(flx_co2bb_glo, flx_co2bb) !-- Release global buffers IF (ALLOCATED(flx_co2ff_glo)) DEALLOCATE(flx_co2ff_glo) IF (ALLOCATED(flx_co2bb_glo)) DEALLOCATE(flx_co2bb_glo) !------------------------------------------------------------------ ! D. SZ98 TARGETS: compute mid-month interpolation targets ! (monthly mode only, each thread processes its own data) ! ! The SZ98 tridiagonal solve finds target values T_j such that ! piecewise-linear interpolation between mid-month points ! preserves the original monthly means M_j exactly: ! M_j = a_j * T_{j-1} + b_j * T_j + c_j * T_{j+1} ! ! This is done ONCE per year. Daily interpolation is O(1). !------------------------------------------------------------------ IF (.NOT. read_daily_co2ff) THEN CALL get_days_in_month(current_year_len, days_in_mth) CALL compute_month_middays(days_in_mth, month_midday) CALL sz98_solve_tridiagonal(flx_co2ff, days_in_mth, 12, klon, target_ff) CALL sz98_solve_tridiagonal(flx_co2bb, days_in_mth, 12, klon, target_bb) IF (is_mpi_root) THEN WRITE(*,*) modname, ': SZ98 targets computed, year_len=', current_year_len WRITE(*,*) ' days_in_month =', days_in_mth WRITE(*,*) ' month_midday =', month_midday END IF END IF END IF ! debutphy !============================================================================= ! 2. DAILY UPDATE: set fco2_ff and fco2_bb for the current day ! ! Guard: only recalculate when the day changes (day_cur /= day_pre_emis). ! Within a single day, the emissions remain constant regardless of how ! many physics time steps occur. This saves unnecessary computation. ! ! IMPORTANT: day_cur from phys_cal_mod is the day of the MONTH (1-31), ! NOT the day of year. For daily-file indexing and SZ98 interpolation, ! we need the day-of-year, computed from mth_cur and day_cur. !============================================================================= IF (day_cur /= day_pre_emis) THEN !-- Compute day-of-year from current month and day doy = compute_day_of_year(mth_cur, day_cur, INT(year_len)) IF (read_daily_co2ff) THEN !--------------------------------------------------------------- ! (A) DAILY MODE: direct indexing by day-of-year !--------------------------------------------------------------- IF (doy >= 1 .AND. doy <= SIZE(flx_co2ff, 2)) THEN fco2_ff(:) = flx_co2ff(:, doy) fco2_bb(:) = flx_co2bb(:, doy) ELSE WRITE(*,*) modname, ': WARNING: doy out of range =', doy fco2_ff(:) = 0.0 fco2_bb(:) = 0.0 END IF ELSE !--------------------------------------------------------------- ! (B) MONTHLY MODE: SZ98 piecewise-linear interpolation ! Interpolates between mid-month target values to give a ! smooth daily field that preserves monthly means. !--------------------------------------------------------------- CALL sz98_interpolate_daily(target_ff, month_midday, doy, & klon, 12, fco2_ff) CALL sz98_interpolate_daily(target_bb, month_midday, doy, & klon, 12, fco2_bb) END IF day_pre_emis = day_cur END IF END SUBROUTINE co2_emissions !=============================================================================== ! FUNCTION: compute_day_of_year ! ! Purpose: ! Convert (month, day_of_month) to day-of-year (1 to year_len). ! Uses the current calendar's month lengths. ! ! Arguments: ! month : Current month (1-12) ! day : Day of the month (1-31) ! year_length : Total number of days in the year (360, 365, or 366) ! ! Returns: ! Day-of-year (integer, 1 to year_length) !=============================================================================== INTEGER FUNCTION compute_day_of_year(month, day, year_length) IMPLICIT NONE INTEGER, INTENT(IN) :: month, day, year_length INTEGER :: dim(12), im CALL get_days_in_month(year_length, dim) compute_day_of_year = day DO im = 1, month - 1 compute_day_of_year = compute_day_of_year + dim(im) END DO END FUNCTION compute_day_of_year !=============================================================================== ! SUBROUTINE: get_days_in_month ! ! Purpose: ! Return the number of days in each month for a given year length. ! Supports 360-day, 365-day (noleap), and 366-day (leap) calendars. ! ! Arguments: ! year_length : Total days in the year (360, 365, or 366) ! days_in_month : Output array of 12 month lengths !=============================================================================== SUBROUTINE get_days_in_month(year_length, days_in_month) IMPLICIT NONE INTEGER, INTENT(IN) :: year_length INTEGER, INTENT(OUT) :: days_in_month(12) SELECT CASE (year_length) CASE (360) !-- 360-day calendar: uniform 30-day months days_in_month(:) = 30 CASE (365) !-- Standard non-leap year (noleap/365_day calendar) days_in_month = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) CASE (366) !-- Leap year (gregorian/standard calendar) days_in_month = (/ 31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) CASE DEFAULT !-- Fallback: assume non-leap year days_in_month = (/ 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 /) END SELECT END SUBROUTINE get_days_in_month !=============================================================================== ! SUBROUTINE: compute_month_middays ! ! Purpose: ! Compute the day-of-year corresponding to the mid-point of each month. ! Used as interpolation nodes for the SZ98 piecewise-linear scheme. ! ! Example (non-leap year): ! January: mid-point = 0 + 31/2 = 15.5 (day 15-16 boundary) ! February: mid-point = 31 + 28/2 = 45.0 (day 45) ! ... ! ! Arguments: ! days_in_month : Number of days in each month (12) ! month_midday : Output: day-of-year of each month mid-point (12) !=============================================================================== SUBROUTINE compute_month_middays(days_in_month, month_midday) IMPLICIT NONE INTEGER, INTENT(IN) :: days_in_month(12) REAL, INTENT(OUT) :: month_midday(12) INTEGER :: cumul(13), im !-- Cumulative days at the START of each month cumul(1) = 0 DO im = 1, 12 cumul(im+1) = cumul(im) + days_in_month(im) END DO !-- Mid-point = start + half the month length DO im = 1, 12 month_midday(im) = REAL(cumul(im)) + REAL(days_in_month(im)) / 2.0 END DO END SUBROUTINE compute_month_middays !=============================================================================== ! SUBROUTINE: sz98_solve_tridiagonal ! ! Purpose: ! Solve the Sheng & Zwiers (1998) tridiagonal system to find mid-month ! target values T that preserve monthly means M under piecewise-linear ! interpolation: ! ! M_j = a_j * T_{j-1} + b_j * T_j + c_j * T_{j+1} ! ! Mathematical Background: ! The coefficients depend on month lengths d_j: ! Interior months (j = 2, ..., 11): ! a_j = d_j / (4 * (d_{j-1} + d_j)) ! c_j = d_j / (4 * (d_j + d_{j+1})) ! b_j = 1 - a_j - c_j ! ! Boundary months: ! j=1: b_1 = (3*d_1+4*d_2)/(4*(d_1+d_2)), c_1 = d_1/(4*(d_1+d_2)) ! j=12: a_N = d_N/(4*(d_{N-1}+d_N)), b_N = (4*d_{N-1}+3*d_N)/(4*(d_{N-1}+d_N)) ! ! When all d_j = 30 (360-day calendar): a=c=1/8, b=3/4 (interior) ! ! Algorithm: ! Thomas algorithm (tridiagonal matrix solver) ! ! Arguments: ! monthly_means : Input monthly mean values (klon_loc, 12) ! days_in_month : Number of days in each month (12) ! n_months : Number of months (always 12) ! klon_loc : Number of local grid points ! targets : Output mid-month target values (klon_loc, 12) ! ! Reference: ! Sheng, J. & Zwiers, F.W. (1998), Climate Dynamics, 14, 609-613. !=============================================================================== SUBROUTINE sz98_solve_tridiagonal(monthly_means, days_in_month, & n_months, klon_loc, targets) IMPLICIT NONE INTEGER, INTENT(IN) :: n_months INTEGER, INTENT(IN) :: klon_loc REAL, INTENT(IN) :: monthly_means(klon_loc, n_months) INTEGER, INTENT(IN) :: days_in_month(n_months) REAL, INTENT(OUT) :: targets(klon_loc, n_months) !-- Tridiagonal coefficients REAL :: a(n_months) ! Sub-diagonal REAL :: b(n_months) ! Main diagonal REAL :: c(n_months) ! Super-diagonal REAL :: d(n_months) ! Month lengths as real !-- Thomas algorithm working arrays REAL :: cp(n_months) ! Modified super-diagonal REAL :: dp(klon_loc, n_months) INTEGER :: j, i REAL :: denom !--------------------------------------------------------------------------- ! Convert month lengths to real for coefficient computation !--------------------------------------------------------------------------- DO j = 1, n_months d(j) = REAL(days_in_month(j)) END DO !--------------------------------------------------------------------------- ! Build tridiagonal matrix (SZ98 variable-length-month formulation) !--------------------------------------------------------------------------- !-- First row (boundary condition: no month before January) a(1) = 0.0 b(1) = (3.0*d(1) + 4.0*d(2)) / (4.0*(d(1) + d(2))) c(1) = d(1) / (4.0*(d(1) + d(2))) !-- Interior rows (February through November) DO j = 2, n_months - 1 a(j) = d(j) / (4.0*(d(j-1) + d(j))) c(j) = d(j) / (4.0*(d(j) + d(j+1))) b(j) = 1.0 - a(j) - c(j) END DO !-- Last row (boundary condition: no month after December) a(n_months) = d(n_months) / (4.0*(d(n_months-1) + d(n_months))) b(n_months) = (4.0*d(n_months-1) + 3.0*d(n_months)) / & (4.0*(d(n_months-1) + d(n_months))) c(n_months) = 0.0 !--------------------------------------------------------------------------- ! Thomas Algorithm: Forward elimination !--------------------------------------------------------------------------- !-- First row cp(1) = c(1) / b(1) DO i = 1, klon_loc dp(i, 1) = monthly_means(i, 1) / b(1) END DO !-- Rows 2 to n_months DO j = 2, n_months denom = b(j) - a(j) * cp(j-1) IF (j < n_months) THEN cp(j) = c(j) / denom END IF DO i = 1, klon_loc dp(i, j) = (monthly_means(i, j) - a(j) * dp(i, j-1)) / denom END DO END DO !--------------------------------------------------------------------------- ! Thomas Algorithm: Back substitution !--------------------------------------------------------------------------- !-- Last row is already solved DO i = 1, klon_loc targets(i, n_months) = dp(i, n_months) END DO !-- Remaining rows (backwards) DO j = n_months - 1, 1, -1 DO i = 1, klon_loc targets(i, j) = dp(i, j) - cp(j) * targets(i, j+1) END DO END DO END SUBROUTINE sz98_solve_tridiagonal !=============================================================================== ! SUBROUTINE: sz98_interpolate_daily ! ! Purpose: ! Interpolate mid-month target values to a daily value using ! piecewise-linear interpolation between month mid-points. ! ! Method: ! 1. Find which two month mid-points bracket the current day-of-year ! 2. Linear interpolation between those targets ! 3. Constant extension before first / after last mid-point ! ! ! Arguments: ! targets : Mid-month target values from sz98_solve_tridiagonal ! month_midday : Day-of-year of each month mid-point (12) ! day_of_year : Current day of year (1 to year_len) ! klon_loc : Number of local grid points ! n_months : Number of months (12) ! daily_value : Output: interpolated emission for current day (klon_loc) !=============================================================================== SUBROUTINE sz98_interpolate_daily(targets, month_midday, day_of_year, & klon_loc, n_months, daily_value) IMPLICIT NONE INTEGER, INTENT(IN) :: klon_loc INTEGER, INTENT(IN) :: n_months REAL, INTENT(IN) :: targets(klon_loc, n_months) REAL, INTENT(IN) :: month_midday(n_months) INTEGER, INTENT(IN) :: day_of_year REAL, INTENT(OUT) :: daily_value(klon_loc) REAL :: day_real, weight INTEGER :: m_before, m_after, im, i day_real = REAL(day_of_year) !--------------------------------------------------------------------------- ! Find bracketing month mid-points !--------------------------------------------------------------------------- IF (day_real <= month_midday(1)) THEN !-- Before first mid-point: constant extension (use January target) DO i = 1, klon_loc daily_value(i) = targets(i, 1) END DO RETURN ELSE IF (day_real >= month_midday(n_months)) THEN !-- After last mid-point: constant extension (use December target) DO i = 1, klon_loc daily_value(i) = targets(i, n_months) END DO RETURN ELSE !-- Find the interval [month_midday(m_before), month_midday(m_after)] m_before = 1 m_after = 2 DO im = 1, n_months - 1 IF (day_real >= month_midday(im) .AND. & day_real < month_midday(im+1)) THEN m_before = im m_after = im + 1 EXIT END IF END DO !-- Linear interpolation between the two bracketing targets weight = (day_real - month_midday(m_before)) / & (month_midday(m_after) - month_midday(m_before)) DO i = 1, klon_loc daily_value(i) = (1.0 - weight) * targets(i, m_before) + & weight * targets(i, m_after) END DO END IF END SUBROUTINE sz98_interpolate_daily !=============================================================================== ! SUBROUTINE: co2_flux_corrections ! ! Purpose: ! Read 2D spatially-varying flux correction fields from NetCDF files and ! distribute them to all MPI ranks and OMP threads via scatter. ! ! These corrections allow the user to apply gridpoint-level adjustments ! to the ocean and/or land net CO2 flux, rather than a single global ! scalar ! ! NetCDF file format: ! dimensions: ! vector = klon_glo (LMDZ 1D grid, same as emission files) ! variables: ! float vector(vector) (coordinate variable) ! float (vector) (correction field) ! units: "kgCO2/m2/s" (same unit as fco2_ocean_cor / fco2_land_cor) ! ! The variable name is configurable via fco2_ocean_cor_2d_var / fco2_land_cor_2d_var ! in physiq.def (default: "fco2_cor"). ! ! Unit convention: ! The 2D field must be in kgCO2/m2/s per unit TOTAL grid-cell area ! (i.e., already weighted by the appropriate surface-type fraction if needed). ! ! Sign convention: ! Positive values REDUCE the source term (same as scalar mode): ! source = fco2_ff + fco2_bb + fco2_land + fco2_ocean ! - fco2_ocean_cor - fco2_land_cor ! ! Called from: ! tracco2i, once at debutphy ! ! Authors: ! 2026-02: P. Cadule ! !=============================================================================== SUBROUTINE co2_flux_corrections(debutphy) USE dimphy USE mod_grid_phy_lmdz, ONLY: klon_glo, nbp_lon, nbp_lat USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root USE mod_phys_lmdz_omp_data, ONLY: is_omp_root USE mod_phys_lmdz_para, ONLY: scatter USE print_control_mod, ONLY: lunout USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_varid, nf95_open USE netcdf, ONLY: nf90_get_var, nf90_nowrite, nf90_noerr, nf90_strerror ! Import 2D correction flags and filenames from carbon_cycle_mod USE carbon_cycle_mod, ONLY: read_fco2_ocean_cor_2d USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_file USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_var USE carbon_cycle_mod, ONLY: read_fco2_land_cor_2d USE carbon_cycle_mod, ONLY: fco2_land_cor_2d_file USE carbon_cycle_mod, ONLY: fco2_land_cor_2d_var ! Import target THREADPRIVATE arrays (already allocated in carbon_cycle_init) ! Section 5 computes: fco2_ocean_cor = base * pctsrf at every time step. USE carbon_cycle_mod, ONLY: fco2_ocean_cor_2d_base USE carbon_cycle_mod, ONLY: fco2_land_cor_2d_base IMPLICIT NONE !----------------------------------------------------------------------------- ! Arguments !----------------------------------------------------------------------------- LOGICAL, INTENT(IN) :: debutphy ! .TRUE. at first physics call !----------------------------------------------------------------------------- ! Local variables for NetCDF reading !----------------------------------------------------------------------------- REAL, ALLOCATABLE :: cor_glo(:) ! Global correction field (klon_glo) REAL, ALLOCATABLE :: vector_coord(:) ! Coordinate variable for validation INTEGER :: ncid_in, varid, ncerr INTEGER :: n_glo CHARACTER(LEN=30), PARAMETER :: modname = 'co2_flux_corrections' CHARACTER(LEN=80) :: abort_message !----------------------------------------------------------------------------- ! first call (debutphy) !----------------------------------------------------------------------------- IF (.NOT. debutphy) RETURN IF (.NOT. read_fco2_ocean_cor_2d .AND. .NOT. read_fco2_land_cor_2d) RETURN !============================================================================= ! A. OCEAN 2D CORRECTION ! ! Read the ocean correction field !============================================================================= IF (read_fco2_ocean_cor_2d) THEN !------------------------------------------------------------------ ! A.1 Read NetCDF !------------------------------------------------------------------ !$OMP MASTER IF (is_mpi_root) THEN !-- Allocate global buffer on MPI root ALLOCATE(cor_glo(klon_glo)) cor_glo(:) = 0.0 !-- Open NetCDF file CALL nf95_open(TRIM(fco2_ocean_cor_2d_file), nf90_nowrite, ncid_in) !-- Validate spatial dimension CALL nf95_inq_varid(ncid_in, 'vector', varid) CALL nf95_gw_var(ncid_in, varid, vector_coord) n_glo = SIZE(vector_coord) DEALLOCATE(vector_coord) IF (n_glo /= klon_glo) THEN WRITE(abort_message, '(A,I8,A,I8)') & 'ocean_cor_2d: vector=', n_glo, ' /= klon_glo=', klon_glo CALL abort_physic(modname, abort_message, 1) END IF !-- Read correction field CALL nf95_inq_varid(ncid_in, TRIM(fco2_ocean_cor_2d_var), varid) ncerr = nf90_get_var(ncid_in, varid, cor_glo) IF (ncerr /= nf90_noerr) THEN abort_message = 'ocean_cor_2d: '//TRIM(nf90_strerror(ncerr)) CALL abort_physic(modname, abort_message, 1) END IF CALL nf95_close(ncid_in) !-- Diagnostic output WRITE(lunout,*) modname, ': --- Ocean 2D correction ---' WRITE(lunout,*) modname, ': File = ', TRIM(fco2_ocean_cor_2d_file) WRITE(lunout,*) modname, ': Var = ', TRIM(fco2_ocean_cor_2d_var) WRITE(lunout,*) modname, ': Min = ', MINVAL(cor_glo), ' kgCO2/m2/s' WRITE(lunout,*) modname, ': Max = ', MAXVAL(cor_glo), ' kgCO2/m2/s' WRITE(lunout,*) modname, ': Mean = ', SUM(cor_glo)/REAL(klon_glo), ' kgCO2/m2/s' END IF ! is_mpi_root !$OMP END MASTER !------------------------------------------------------------------ ! A.2 Dummy allocation for non-root ranks !------------------------------------------------------------------ IF (.NOT. ALLOCATED(cor_glo)) ALLOCATE(cor_glo(0)) !------------------------------------------------------------------ ! A.3 Scatter: distribute global field to all ranks/threads ! Section 5 computes fco2_ocean_cor = base * pctsrf every step. !------------------------------------------------------------------ CALL scatter(cor_glo, fco2_ocean_cor_2d_base) !-- Release global buffer immediately IF (ALLOCATED(cor_glo)) DEALLOCATE(cor_glo) END IF ! read_fco2_ocean_cor_2d !============================================================================= ! B. LAND 2D CORRECTION ! ! Same pattern as ocean, targeting fco2_land_cor(:). !============================================================================= IF (read_fco2_land_cor_2d) THEN !------------------------------------------------------------------ ! B.1 Read NetCDF !------------------------------------------------------------------ !$OMP MASTER IF (is_mpi_root) THEN !-- Allocate global buffer on MPI root ALLOCATE(cor_glo(klon_glo)) cor_glo(:) = 0.0 !-- Open NetCDF file CALL nf95_open(TRIM(fco2_land_cor_2d_file), nf90_nowrite, ncid_in) !-- Validate spatial dimension CALL nf95_inq_varid(ncid_in, 'vector', varid) CALL nf95_gw_var(ncid_in, varid, vector_coord) n_glo = SIZE(vector_coord) DEALLOCATE(vector_coord) IF (n_glo /= klon_glo) THEN WRITE(abort_message, '(A,I8,A,I8)') & 'land_cor_2d: vector=', n_glo, ' /= klon_glo=', klon_glo CALL abort_physic(modname, abort_message, 1) END IF !-- Read correction field CALL nf95_inq_varid(ncid_in, TRIM(fco2_land_cor_2d_var), varid) ncerr = nf90_get_var(ncid_in, varid, cor_glo) IF (ncerr /= nf90_noerr) THEN abort_message = 'land_cor_2d: '//TRIM(nf90_strerror(ncerr)) CALL abort_physic(modname, abort_message, 1) END IF CALL nf95_close(ncid_in) !-- Diagnostic output WRITE(lunout,*) modname, ': --- Land 2D correction ---' WRITE(lunout,*) modname, ': File = ', TRIM(fco2_land_cor_2d_file) WRITE(lunout,*) modname, ': Var = ', TRIM(fco2_land_cor_2d_var) WRITE(lunout,*) modname, ': Min = ', MINVAL(cor_glo), ' kgCO2/m2/s' WRITE(lunout,*) modname, ': Max = ', MAXVAL(cor_glo), ' kgCO2/m2/s' WRITE(lunout,*) modname, ': Mean = ', SUM(cor_glo)/REAL(klon_glo), ' kgCO2/m2/s' END IF ! is_mpi_root !$OMP END MASTER !------------------------------------------------------------------ ! B.2 Dummy allocation !------------------------------------------------------------------ IF (.NOT. ALLOCATED(cor_glo)) ALLOCATE(cor_glo(0)) !------------------------------------------------------------------ ! B.3 Scatter to fco2_land_cor_2d_base !------------------------------------------------------------------ CALL scatter(cor_glo, fco2_land_cor_2d_base) !-- Release global buffer immediately IF (ALLOCATED(cor_glo)) DEALLOCATE(cor_glo) END IF ! read_fco2_land_cor_2d END SUBROUTINE co2_flux_corrections !=============================================================================== ! SUBROUTINE: reorder_global_1d ! ! Purpose: ! Reorder a 1D global physics field (klon_glo) to correct grid ordering ! convention mismatches between external NetCDF files and the LMDZ internal ! physics grid. ! ! LMDZ Physics Grid Structure: ! The LMDZ physics grid is NOT a regular (nlon x nlat) rectangle. ! It has SINGULAR POLE POINTS: ! ! field(1) = North Pole (1 point) ! field(2 : nlon+1) = 1st interior band (N ==> S) ! field(nlon+2 : 2*nlon+1) = 2nd interior band ! ... ! field((nlat_int-1)*nlon+2 : nlat_int*nlon+1) = last interior band ! field(klon_glo) = South Pole (1 point) ! ! where: ! nlat_int = nlat - 2 (number of interior latitude bands) ! klon_glo = nlon * nlat_int + 2 = nlon * (nlat - 2) + 2 ! ! THIS IS NOT nlon * nlat. ! ! ! Arguments: ! field(klon) : INOUT - flat 1D field in LMDZ physics layout ! klon = nlon * (nlat - 2) + 2 ! nlon : IN - longitude points per interior band (nbp_lon) ! nlat : IN - total latitude count including poles (nbp_lat) ! flip_lat : IN - .TRUE. to reverse latitude band ordering ! (S=>N becomes N=>S, or vice versa) ! Also swaps pole values. ! shift_lon : IN - circular shift in longitude (grid points, interior only) ! Typical: nlon/2 for [-180,+180] => [0,360] conversion ! Positive = left rotation (eastward shift of origin) ! Poles are scalar averages: unaffected by longitude shift. ! ! Algorithm: ! Step 1: If flip_lat - swap pole values (indices 1 and klon_glo) ! Step 2: Copy interior bands to 2D work array ! Step 3: Write back with reversed band order (if flip_lat) ! and/or CSHIFT (if shift_lon /= 0) ! Step 4: the field is now in LMDZ internal ordering ! ! ! Authors: ! 2026-02: P. Cadule - IPSL !=============================================================================== SUBROUTINE reorder_global_1d(field, nlon, nlat, flip_lat, shift_lon) IMPLICIT NONE !-- Arguments INTEGER, INTENT(IN) :: nlon ! Number of longitudes (nbp_lon) INTEGER, INTENT(IN) :: nlat ! Total latitudes incl. poles (nbp_lat) REAL, INTENT(INOUT) :: field(nlon * (nlat - 2) + 2) ! LMDZ physics 1D field LOGICAL, INTENT(IN) :: flip_lat ! Flip latitude band ordering INTEGER, INTENT(IN) :: shift_lon ! Longitude circular shift (grid points) !-- Local variables INTEGER :: nlat_int ! Number of interior latitude bands INTEGER :: klon_glo ! Total field size (for clarity) REAL :: work(nlon, nlat - 2) ! 2D work array for interior bands REAL :: temp_pole ! Temporary for pole swap INTEGER :: j, idx_src, idx_dst !-- Quick exit if no reordering requested IF (.NOT. flip_lat .AND. shift_lon == 0) RETURN nlat_int = nlat - 2 klon_glo = nlon * nlat_int + 2 !--------------------------------------------------------------------------- ! Step 1: If flipping latitudes, swap pole points ! ! LMDZ convention: ! field(1) = North Pole ! field(klon_glo) = South Pole ! ! File convention (CF, S=>N): ! field(1) = South Pole (needs to become North Pole) ! field(klon_glo) = North Pole (needs to become South Pole) !--------------------------------------------------------------------------- IF (flip_lat) THEN temp_pole = field(1) field(1) = field(klon_glo) field(klon_glo) = temp_pole END IF !--------------------------------------------------------------------------- ! Step 2: Copy interior bands to 2D work array ! ! Interior points span indices 2 to klon_glo-1, arranged as: ! band j: field(2 + (j-1)*nlon) .. field(1 + j*nlon) ! for j = 1 .. nlat_int !--------------------------------------------------------------------------- DO j = 1, nlat_int idx_src = 2 + (j - 1) * nlon work(1:nlon, j) = field(idx_src : idx_src + nlon - 1) END DO !--------------------------------------------------------------------------- ! Step 3: Write back with reversed band order and/or longitude shift ! ! If flip_lat: output band j reads from work band (nlat_int - j + 1) ! If shift_lon: each band is circularly shifted by shift_lon positions ! using Fortran intrinsic CSHIFT (left rotation) !--------------------------------------------------------------------------- DO j = 1, nlat_int idx_dst = 2 + (j - 1) * nlon IF (flip_lat) THEN IF (shift_lon /= 0) THEN field(idx_dst : idx_dst + nlon - 1) = & CSHIFT(work(1:nlon, nlat_int - j + 1), shift_lon) ELSE field(idx_dst : idx_dst + nlon - 1) = work(1:nlon, nlat_int - j + 1) END IF ELSE ! No latitude flip, longitude shift only field(idx_dst : idx_dst + nlon - 1) = CSHIFT(work(1:nlon, j), shift_lon) END IF END DO END SUBROUTINE reorder_global_1d END MODULE tracco2i_mod