module tracer_mass_fixer_dyn_mod implicit none logical, save :: call_mass_fixer_dyn = .false. ! flag for correcting mass non-conservations due to dynamics logical, save :: found_startfi_co = .false. ! flags for finding mass_predyn_ in startfi.nc logical, save :: found_startfi_o2 = .false. logical, save :: found_startfi_h2 = .false. logical, save :: found_startfi_ho2 = .false. logical, save :: found_startfi_h2o2 = .false. logical, save :: found_startfi_n2 = .false. logical, save :: found_startfi_ar = .false. logical, save :: found_startfi_he = .false. real, allocatable, save :: mass_predyn(:) ! tracer mass at the end of physiq before dynamics starts ! tracer mass at the end of the last physiq step for writing and reading startfi real, save :: mass_predyn_co, mass_predyn_o2, mass_predyn_h2, mass_predyn_ho2, & mass_predyn_h2o2, mass_predyn_n2, mass_predyn_ar, mass_predyn_he !$omp threadprivate(call_mass_fixer_dyn) !$omp threadprivate(found_startfi_co, found_startfi_o2, found_startfi_h2, found_startfi_ho2) !$omp threadprivate(found_startfi_h2o2, found_startfi_n2,found_startfi_ar,found_startfi_he) !$omp threadprivate(mass_predyn) !$omp threadprivate(mass_predyn_co, mass_predyn_o2, mass_predyn_h2, mass_predyn_ho2) !$omp threadprivate(mass_predyn_h2o2, mass_predyn_n2, mass_predyn_ar, mass_predyn_he) contains subroutine tracer_mass_fixer_dyn(ngrid, nlayer, nq, ptimestep, pq, zplev, firstcall, pdq_corrdyn) !-------------------------------------------------------------------------------- ! Correction of mass mixing ratios for tracers that represent a significant fraction of their constituent elements ! at the beginning of physiq, to ensure mass conservation during dynamics. ! Y. LUO 03/2026 !--------------------------------------------------------------------------------- use tracer_mod, only: noms, igcm_co, igcm_o2, igcm_h2, igcm_ho2, igcm_h2o2, igcm_n2, igcm_ar, igcm_he implicit none integer, intent(in) :: ngrid ! number of atmospheric columns integer, intent(in) :: nlayer ! number of atmospheric layers integer, intent(in) :: nq ! number of tracers in traceur.def real, intent(in) :: ptimestep ! physics timestep (s) real, intent(in) :: pq(ngrid, nlayer, nq) ! tracer mass mixing ratios (kg per 1 kg of air) real, intent(in) :: zplev(ngrid, nlayer + 1) ! inter-layer pressure (Pa) logical, intent(in) :: firstcall ! True at the first call real, intent(out) :: pdq_corrdyn(ngrid, nlayer, nq) ! tracer mass mixing ratio tendencies in order to correct ! mass non-conservations in dynamics real mass_postdyn(nq) ! tracer mass when dynamics ends and physics starts real masscorrfac(nq) ! correction factor for mass non-conservations in dynamics real pq_corrdyn(ngrid, nlayer, nq) ! tracer mass mixing ratio after correcting mass non-conservations in dynamics call compute_tracer_mass_global(ngrid, nlayer, nq, pq, zplev, mass_postdyn) pq_corrdyn = pq if (igcm_co .ne. 0 .and. (.not. firstcall .or. found_startfi_co)) then masscorrfac(igcm_co) = mass_predyn(igcm_co) / mass_postdyn(igcm_co) pq_corrdyn(:, :, igcm_co) = pq(:, :, igcm_co) * masscorrfac(igcm_co) endif if (igcm_o2 .ne. 0 .and. (.not. firstcall .or. found_startfi_o2)) then masscorrfac(igcm_o2) = mass_predyn(igcm_o2) / mass_postdyn(igcm_o2) pq_corrdyn(:, :, igcm_o2) = pq(:, :, igcm_o2) * masscorrfac(igcm_o2) endif if (igcm_h2 .ne. 0 .and. (.not. firstcall .or. found_startfi_h2)) then masscorrfac(igcm_h2) = mass_predyn(igcm_h2) / mass_postdyn(igcm_h2) pq_corrdyn(:, :, igcm_h2) = pq(:, :, igcm_h2) * masscorrfac(igcm_h2) endif if (igcm_ho2 .ne. 0 .and. (.not. firstcall .or. found_startfi_ho2)) then masscorrfac(igcm_ho2) = mass_predyn(igcm_ho2) / mass_postdyn(igcm_ho2) pq_corrdyn(:, :, igcm_ho2) = pq(:, :, igcm_ho2) * masscorrfac(igcm_ho2) endif if (igcm_h2o2 .ne. 0 .and. (.not. firstcall .or. found_startfi_h2o2)) then masscorrfac(igcm_h2o2) = mass_predyn(igcm_h2o2) / mass_postdyn(igcm_h2o2) pq_corrdyn(:, :, igcm_h2o2) = pq(:, :, igcm_h2o2) * masscorrfac(igcm_h2o2) endif if (igcm_n2 .ne. 0 .and. (.not. firstcall .or. found_startfi_n2)) then masscorrfac(igcm_n2) = mass_predyn(igcm_n2) / mass_postdyn(igcm_n2) pq_corrdyn(:, :, igcm_n2) = pq(:, :, igcm_n2) * masscorrfac(igcm_n2) endif if (igcm_ar .ne. 0 .and. (.not. firstcall .or. found_startfi_ar)) then masscorrfac(igcm_ar) = mass_predyn(igcm_ar) / mass_postdyn(igcm_ar) pq_corrdyn(:, :, igcm_ar) = pq(:, :, igcm_ar) * masscorrfac(igcm_ar) endif if (igcm_he .ne. 0 .and. (.not. firstcall .or. found_startfi_he)) then masscorrfac(igcm_he) = mass_predyn(igcm_he) / mass_postdyn(igcm_he) pq_corrdyn(:, :, igcm_he) = pq(:, :, igcm_he) * masscorrfac(igcm_he) endif pdq_corrdyn = (pq_corrdyn - pq) / ptimestep end subroutine tracer_mass_fixer_dyn subroutine compute_tracer_mass_global(ngrid, nlayer, nq, pq, zplev, mass) use tracer_mod, only: mmol use comcstfi_h, only: g use geometry_mod, only: cell_area use tracer_mod, only: noms use planetwide_mod, only: planetwide_sumval implicit none integer, intent(in) :: ngrid ! number of atmospheric columns integer, intent(in) :: nlayer ! number of atmospheric layers integer, intent(in) :: nq ! number of tracers in traceur.def real, intent(in) :: pq(ngrid, nlayer, nq) ! tracer mass mixing ratios (kg per 1 kg of air) real, intent(in) :: zplev(ngrid, nlayer + 1) ! inter-layer pressure (Pa) real, intent(out) :: mass(nq) ! tracer mass in the whole atmosphere real colden(ngrid, nq) ! tracer column density (molecules cm-2) integer iq, ig, l do iq = 1, nq if (noms(iq) .eq. "co" .or. & noms(iq) .eq. "o2" .or. & noms(iq) .eq. "h2" .or. & noms(iq) .eq. "ho2" .or. & noms(iq) .eq. "h2o2" .or. & noms(iq) .eq. "n2" .or. & noms(iq) .eq. "ar" .or. & noms(iq) .eq. "he") then do ig = 1, ngrid colden(ig, iq) = 0. enddo do l = 1, nlayer do ig = 1, ngrid colden(ig, iq) = colden(ig, iq) + & pq(ig, l, iq) * & (zplev(ig, l) - zplev(ig, l+1)) * & 6.022e22 / (mmol(iq) * g) enddo enddo call planetwide_sumval( & colden(:, iq) / 6.022e23 * mmol(iq) * 1.e4 * cell_area, & mass(iq) ) endif enddo end subroutine compute_tracer_mass_global subroutine set_mass_predyn_from_startfi use tracer_mod, only: igcm_co, igcm_o2, igcm_h2, igcm_ho2, & igcm_h2o2, igcm_n2, igcm_ar, igcm_he implicit none if (igcm_co.ne.0 .and. found_startfi_co) then mass_predyn(igcm_co) = mass_predyn_co endif if (igcm_o2.ne.0 .and. found_startfi_o2) then mass_predyn(igcm_o2) = mass_predyn_o2 endif if (igcm_h2.ne.0 .and. found_startfi_h2) then mass_predyn(igcm_h2) = mass_predyn_h2 endif if (igcm_ho2.ne.0 .and. found_startfi_ho2) then mass_predyn(igcm_ho2) = mass_predyn_ho2 endif if (igcm_h2o2.ne.0 .and. found_startfi_h2o2) then mass_predyn(igcm_h2o2) = mass_predyn_h2o2 endif if (igcm_n2.ne.0 .and. found_startfi_n2) then mass_predyn(igcm_n2) = mass_predyn_n2 endif if (igcm_ar.ne.0 .and. found_startfi_ar) then mass_predyn(igcm_ar) = mass_predyn_ar endif if (igcm_he.ne.0 .and. found_startfi_he) then mass_predyn(igcm_he) = mass_predyn_he endif end subroutine set_mass_predyn_from_startfi subroutine ini_tracer_mass_fixer_dyn(nq) implicit none integer, intent(in) :: nq ! number of tracers in traceur.def allocate(mass_predyn(nq)) mass_predyn(1:nq) = -999. ! newstart writes mass_predyn_ in startfi.nc with this value, ! ensuring mass_predyn_ <= 0. As a result, found_startfi_ ! remains .false., and no correction is applied at the first physics step ! when initialized from newstart output. end subroutine ini_tracer_mass_fixer_dyn subroutine end_tracer_mass_fixer_dyn implicit none if (allocated(mass_predyn)) deallocate(mass_predyn) end subroutine end_tracer_mass_fixer_dyn end module tracer_mass_fixer_dyn_mod