| 1 | module tracer_mass_fixer_dyn_mod |
|---|
| 2 | |
|---|
| 3 | implicit none |
|---|
| 4 | |
|---|
| 5 | logical, save :: call_mass_fixer_dyn = .false. ! flag for correcting mass non-conservations due to dynamics |
|---|
| 6 | logical, save :: found_startfi_co = .false. ! flags for finding mass_predyn_<tracer name> in startfi.nc |
|---|
| 7 | logical, save :: found_startfi_o2 = .false. |
|---|
| 8 | logical, save :: found_startfi_h2 = .false. |
|---|
| 9 | logical, save :: found_startfi_ho2 = .false. |
|---|
| 10 | logical, save :: found_startfi_h2o2 = .false. |
|---|
| 11 | logical, save :: found_startfi_n2 = .false. |
|---|
| 12 | logical, save :: found_startfi_ar = .false. |
|---|
| 13 | logical, save :: found_startfi_he = .false. |
|---|
| 14 | real, allocatable, save :: mass_predyn(:) ! tracer mass at the end of physiq before dynamics starts |
|---|
| 15 | ! tracer mass at the end of the last physiq step for writing and reading startfi |
|---|
| 16 | real, save :: mass_predyn_co, mass_predyn_o2, mass_predyn_h2, mass_predyn_ho2, & |
|---|
| 17 | mass_predyn_h2o2, mass_predyn_n2, mass_predyn_ar, mass_predyn_he |
|---|
| 18 | |
|---|
| 19 | !$omp threadprivate(call_mass_fixer_dyn) |
|---|
| 20 | !$omp threadprivate(found_startfi_co, found_startfi_o2, found_startfi_h2, found_startfi_ho2) |
|---|
| 21 | !$omp threadprivate(found_startfi_h2o2, found_startfi_n2,found_startfi_ar,found_startfi_he) |
|---|
| 22 | !$omp threadprivate(mass_predyn) |
|---|
| 23 | !$omp threadprivate(mass_predyn_co, mass_predyn_o2, mass_predyn_h2, mass_predyn_ho2) |
|---|
| 24 | !$omp threadprivate(mass_predyn_h2o2, mass_predyn_n2, mass_predyn_ar, mass_predyn_he) |
|---|
| 25 | |
|---|
| 26 | contains |
|---|
| 27 | |
|---|
| 28 | subroutine tracer_mass_fixer_dyn(ngrid, nlayer, nq, ptimestep, pq, zplev, firstcall, pdq_corrdyn) |
|---|
| 29 | |
|---|
| 30 | !-------------------------------------------------------------------------------- |
|---|
| 31 | ! Correction of mass mixing ratios for tracers that represent a significant fraction of their constituent elements |
|---|
| 32 | ! at the beginning of physiq, to ensure mass conservation during dynamics. |
|---|
| 33 | ! Y. LUO 03/2026 |
|---|
| 34 | !--------------------------------------------------------------------------------- |
|---|
| 35 | |
|---|
| 36 | use tracer_mod, only: noms, igcm_co, igcm_o2, igcm_h2, igcm_ho2, igcm_h2o2, igcm_n2, igcm_ar, igcm_he |
|---|
| 37 | |
|---|
| 38 | implicit none |
|---|
| 39 | |
|---|
| 40 | integer, intent(in) :: ngrid ! number of atmospheric columns |
|---|
| 41 | integer, intent(in) :: nlayer ! number of atmospheric layers |
|---|
| 42 | integer, intent(in) :: nq ! number of tracers in traceur.def |
|---|
| 43 | real, intent(in) :: ptimestep ! physics timestep (s) |
|---|
| 44 | real, intent(in) :: pq(ngrid, nlayer, nq) ! tracer mass mixing ratios (kg per 1 kg of air) |
|---|
| 45 | real, intent(in) :: zplev(ngrid, nlayer + 1) ! inter-layer pressure (Pa) |
|---|
| 46 | logical, intent(in) :: firstcall ! True at the first call |
|---|
| 47 | real, intent(out) :: pdq_corrdyn(ngrid, nlayer, nq) ! tracer mass mixing ratio tendencies in order to correct |
|---|
| 48 | ! mass non-conservations in dynamics |
|---|
| 49 | |
|---|
| 50 | real mass_postdyn(nq) ! tracer mass when dynamics ends and physics starts |
|---|
| 51 | real masscorrfac(nq) ! correction factor for mass non-conservations in dynamics |
|---|
| 52 | real pq_corrdyn(ngrid, nlayer, nq) ! tracer mass mixing ratio after correcting mass non-conservations in dynamics |
|---|
| 53 | |
|---|
| 54 | call compute_tracer_mass_global(ngrid, nlayer, nq, pq, zplev, mass_postdyn) |
|---|
| 55 | |
|---|
| 56 | pq_corrdyn = pq |
|---|
| 57 | |
|---|
| 58 | if (igcm_co .ne. 0 .and. (.not. firstcall .or. found_startfi_co)) then |
|---|
| 59 | masscorrfac(igcm_co) = mass_predyn(igcm_co) / mass_postdyn(igcm_co) |
|---|
| 60 | pq_corrdyn(:, :, igcm_co) = pq(:, :, igcm_co) * masscorrfac(igcm_co) |
|---|
| 61 | endif |
|---|
| 62 | |
|---|
| 63 | if (igcm_o2 .ne. 0 .and. (.not. firstcall .or. found_startfi_o2)) then |
|---|
| 64 | masscorrfac(igcm_o2) = mass_predyn(igcm_o2) / mass_postdyn(igcm_o2) |
|---|
| 65 | pq_corrdyn(:, :, igcm_o2) = pq(:, :, igcm_o2) * masscorrfac(igcm_o2) |
|---|
| 66 | endif |
|---|
| 67 | |
|---|
| 68 | if (igcm_h2 .ne. 0 .and. (.not. firstcall .or. found_startfi_h2)) then |
|---|
| 69 | masscorrfac(igcm_h2) = mass_predyn(igcm_h2) / mass_postdyn(igcm_h2) |
|---|
| 70 | pq_corrdyn(:, :, igcm_h2) = pq(:, :, igcm_h2) * masscorrfac(igcm_h2) |
|---|
| 71 | endif |
|---|
| 72 | |
|---|
| 73 | if (igcm_ho2 .ne. 0 .and. (.not. firstcall .or. found_startfi_ho2)) then |
|---|
| 74 | masscorrfac(igcm_ho2) = mass_predyn(igcm_ho2) / mass_postdyn(igcm_ho2) |
|---|
| 75 | pq_corrdyn(:, :, igcm_ho2) = pq(:, :, igcm_ho2) * masscorrfac(igcm_ho2) |
|---|
| 76 | endif |
|---|
| 77 | |
|---|
| 78 | if (igcm_h2o2 .ne. 0 .and. (.not. firstcall .or. found_startfi_h2o2)) then |
|---|
| 79 | masscorrfac(igcm_h2o2) = mass_predyn(igcm_h2o2) / mass_postdyn(igcm_h2o2) |
|---|
| 80 | pq_corrdyn(:, :, igcm_h2o2) = pq(:, :, igcm_h2o2) * masscorrfac(igcm_h2o2) |
|---|
| 81 | endif |
|---|
| 82 | |
|---|
| 83 | if (igcm_n2 .ne. 0 .and. (.not. firstcall .or. found_startfi_n2)) then |
|---|
| 84 | masscorrfac(igcm_n2) = mass_predyn(igcm_n2) / mass_postdyn(igcm_n2) |
|---|
| 85 | pq_corrdyn(:, :, igcm_n2) = pq(:, :, igcm_n2) * masscorrfac(igcm_n2) |
|---|
| 86 | endif |
|---|
| 87 | |
|---|
| 88 | if (igcm_ar .ne. 0 .and. (.not. firstcall .or. found_startfi_ar)) then |
|---|
| 89 | masscorrfac(igcm_ar) = mass_predyn(igcm_ar) / mass_postdyn(igcm_ar) |
|---|
| 90 | pq_corrdyn(:, :, igcm_ar) = pq(:, :, igcm_ar) * masscorrfac(igcm_ar) |
|---|
| 91 | endif |
|---|
| 92 | |
|---|
| 93 | if (igcm_he .ne. 0 .and. (.not. firstcall .or. found_startfi_he)) then |
|---|
| 94 | masscorrfac(igcm_he) = mass_predyn(igcm_he) / mass_postdyn(igcm_he) |
|---|
| 95 | pq_corrdyn(:, :, igcm_he) = pq(:, :, igcm_he) * masscorrfac(igcm_he) |
|---|
| 96 | endif |
|---|
| 97 | |
|---|
| 98 | pdq_corrdyn = (pq_corrdyn - pq) / ptimestep |
|---|
| 99 | |
|---|
| 100 | end subroutine tracer_mass_fixer_dyn |
|---|
| 101 | |
|---|
| 102 | subroutine compute_tracer_mass_global(ngrid, nlayer, nq, pq, zplev, mass) |
|---|
| 103 | |
|---|
| 104 | use tracer_mod, only: mmol |
|---|
| 105 | use comcstfi_h, only: g |
|---|
| 106 | use geometry_mod, only: cell_area |
|---|
| 107 | use tracer_mod, only: noms |
|---|
| 108 | use planetwide_mod, only: planetwide_sumval |
|---|
| 109 | |
|---|
| 110 | implicit none |
|---|
| 111 | |
|---|
| 112 | integer, intent(in) :: ngrid ! number of atmospheric columns |
|---|
| 113 | integer, intent(in) :: nlayer ! number of atmospheric layers |
|---|
| 114 | integer, intent(in) :: nq ! number of tracers in traceur.def |
|---|
| 115 | real, intent(in) :: pq(ngrid, nlayer, nq) ! tracer mass mixing ratios (kg per 1 kg of air) |
|---|
| 116 | real, intent(in) :: zplev(ngrid, nlayer + 1) ! inter-layer pressure (Pa) |
|---|
| 117 | real, intent(out) :: mass(nq) ! tracer mass in the whole atmosphere |
|---|
| 118 | |
|---|
| 119 | real colden(ngrid, nq) ! tracer column density (molecules cm-2) |
|---|
| 120 | integer iq, ig, l |
|---|
| 121 | |
|---|
| 122 | do iq = 1, nq |
|---|
| 123 | if (noms(iq) .eq. "co" .or. & |
|---|
| 124 | noms(iq) .eq. "o2" .or. & |
|---|
| 125 | noms(iq) .eq. "h2" .or. & |
|---|
| 126 | noms(iq) .eq. "ho2" .or. & |
|---|
| 127 | noms(iq) .eq. "h2o2" .or. & |
|---|
| 128 | noms(iq) .eq. "n2" .or. & |
|---|
| 129 | noms(iq) .eq. "ar" .or. & |
|---|
| 130 | noms(iq) .eq. "he") then |
|---|
| 131 | |
|---|
| 132 | do ig = 1, ngrid |
|---|
| 133 | colden(ig, iq) = 0. |
|---|
| 134 | enddo |
|---|
| 135 | |
|---|
| 136 | do l = 1, nlayer |
|---|
| 137 | do ig = 1, ngrid |
|---|
| 138 | colden(ig, iq) = colden(ig, iq) + & |
|---|
| 139 | pq(ig, l, iq) * & |
|---|
| 140 | (zplev(ig, l) - zplev(ig, l+1)) * & |
|---|
| 141 | 6.022e22 / (mmol(iq) * g) |
|---|
| 142 | enddo |
|---|
| 143 | enddo |
|---|
| 144 | |
|---|
| 145 | call planetwide_sumval( & |
|---|
| 146 | colden(:, iq) / 6.022e23 * mmol(iq) * 1.e4 * cell_area, & |
|---|
| 147 | mass(iq) ) |
|---|
| 148 | |
|---|
| 149 | endif |
|---|
| 150 | enddo |
|---|
| 151 | |
|---|
| 152 | end subroutine compute_tracer_mass_global |
|---|
| 153 | |
|---|
| 154 | subroutine set_mass_predyn_from_startfi |
|---|
| 155 | |
|---|
| 156 | use tracer_mod, only: igcm_co, igcm_o2, igcm_h2, igcm_ho2, & |
|---|
| 157 | igcm_h2o2, igcm_n2, igcm_ar, igcm_he |
|---|
| 158 | |
|---|
| 159 | implicit none |
|---|
| 160 | |
|---|
| 161 | if (igcm_co.ne.0 .and. found_startfi_co) then |
|---|
| 162 | mass_predyn(igcm_co) = mass_predyn_co |
|---|
| 163 | endif |
|---|
| 164 | |
|---|
| 165 | if (igcm_o2.ne.0 .and. found_startfi_o2) then |
|---|
| 166 | mass_predyn(igcm_o2) = mass_predyn_o2 |
|---|
| 167 | endif |
|---|
| 168 | |
|---|
| 169 | if (igcm_h2.ne.0 .and. found_startfi_h2) then |
|---|
| 170 | mass_predyn(igcm_h2) = mass_predyn_h2 |
|---|
| 171 | endif |
|---|
| 172 | |
|---|
| 173 | if (igcm_ho2.ne.0 .and. found_startfi_ho2) then |
|---|
| 174 | mass_predyn(igcm_ho2) = mass_predyn_ho2 |
|---|
| 175 | endif |
|---|
| 176 | |
|---|
| 177 | if (igcm_h2o2.ne.0 .and. found_startfi_h2o2) then |
|---|
| 178 | mass_predyn(igcm_h2o2) = mass_predyn_h2o2 |
|---|
| 179 | endif |
|---|
| 180 | |
|---|
| 181 | if (igcm_n2.ne.0 .and. found_startfi_n2) then |
|---|
| 182 | mass_predyn(igcm_n2) = mass_predyn_n2 |
|---|
| 183 | endif |
|---|
| 184 | |
|---|
| 185 | if (igcm_ar.ne.0 .and. found_startfi_ar) then |
|---|
| 186 | mass_predyn(igcm_ar) = mass_predyn_ar |
|---|
| 187 | endif |
|---|
| 188 | |
|---|
| 189 | if (igcm_he.ne.0 .and. found_startfi_he) then |
|---|
| 190 | mass_predyn(igcm_he) = mass_predyn_he |
|---|
| 191 | endif |
|---|
| 192 | |
|---|
| 193 | end subroutine set_mass_predyn_from_startfi |
|---|
| 194 | |
|---|
| 195 | subroutine ini_tracer_mass_fixer_dyn(nq) |
|---|
| 196 | |
|---|
| 197 | implicit none |
|---|
| 198 | |
|---|
| 199 | integer, intent(in) :: nq ! number of tracers in traceur.def |
|---|
| 200 | |
|---|
| 201 | allocate(mass_predyn(nq)) |
|---|
| 202 | mass_predyn(1:nq) = -999. ! newstart writes mass_predyn_<tracer> in startfi.nc with this value, |
|---|
| 203 | ! ensuring mass_predyn_<tracer> <= 0. As a result, found_startfi_<tracer> |
|---|
| 204 | ! remains .false., and no correction is applied at the first physics step |
|---|
| 205 | ! when initialized from newstart output. |
|---|
| 206 | |
|---|
| 207 | end subroutine ini_tracer_mass_fixer_dyn |
|---|
| 208 | |
|---|
| 209 | subroutine end_tracer_mass_fixer_dyn |
|---|
| 210 | |
|---|
| 211 | implicit none |
|---|
| 212 | |
|---|
| 213 | if (allocated(mass_predyn)) deallocate(mass_predyn) |
|---|
| 214 | |
|---|
| 215 | end subroutine end_tracer_mass_fixer_dyn |
|---|
| 216 | |
|---|
| 217 | end module tracer_mass_fixer_dyn_mod |
|---|
