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) ! ! ! 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. !=============================================================================== 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() ENDIF END SUBROUTINE tracco2i_init !=============================================================================== ! 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 USE carbon_cycle_mod, ONLY: id_CO2 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: 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 mod_grid_phy_lmdz USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root USE mod_phys_lmdz_para, ONLY: gather, bcast, scatter, reduce_sum 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 ! Local sums for MPI reduction REAL :: co2_mass_local, air_mass_local REAL :: co2_mass_global, air_mass_global REAL, PARAMETER :: secinday = 86400.0 !=========================================================================== ! 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,*) '==============================================' ENDIF 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 ! Emission-driven: initialise only if essentially zero IF (MAXVAL(ABS(tr_seri(:,:,id_CO2))) < 1.0e-15) THEN tr_seri(:,:,id_CO2) = co2_ppm0 * 1.0e-6 * (RMCO2 / RMD) END IF END IF ENDIF ! 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. ENDDO ! 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) ENDIF 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) ENDDO 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(:) ENDDO !=========================================================================== ! 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 ! 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 ENDDO !=========================================================================== ! 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)) ENDIF ! (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)) ENDIF ! --- 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 !=========================================================================== ! 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(:) 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' ENDIF ELSE IF (co2_tracer_available) THEN !--------------------------------------------------------------------- ! Emission-driven with tracer available: compute from tracer field !--------------------------------------------------------------------- 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 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 ENDIF 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) ENDIF END DO END FUNCTION to_lower !=============================================================================== ! 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) ENDIF !-- 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) ENDIF !-- 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 ENDIF !--------------------------------------------------------------- ! 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) ENDIF 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) ENDIF 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 ENDIF 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 ENDIF 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