Changeset 4150 for trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

Timestamp:

Mar 24, 2026, 5:27:57 PM (3 weeks ago)

Author:

yluo

Message:

Mars PCM:
Refactored the tracer mass-conservation fixer for dynamics introduced in r4142:
Added a runtime flag to enable/disable the correction of tracer mass non-conservation in dynamics.
Reorganized the implementation into a modular structure.
Restricted the correction to a subset of tracers: CO, O2, H2, HO2, H2O2, N2, Ar, and He.
Enabled correction for dynamics timesteps preceding the first physics timestep of a simulation.
YCL

File:

• trunk/LMDZ.MARS/libf/phymars/physiq_mod.F (modified) (7 diffs)

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -275 +275 @@
       REAL,INTENT(in) :: pt(ngrid,nlayer) ! temperature (K)
       REAL,INTENT(in) :: pq(ngrid,nlayer,nq) ! tracers (.../kg_of_air)
-      REAL pq_corrdyn(ngrid,nlayer,nq) ! tracer mass mixing ratio after correcting
-                                       ! mass non-conservations in dynamics
       REAL,INTENT(in) :: flxw(ngrid,nlayer) ! vertical mass flux (ks/s) at lower boundary of layer
 
@@ -286 +284 @@
       REAL,INTENT(out) :: pdt(ngrid,nlayer) ! temperature tendency (K/s)
       REAL,INTENT(out) :: pdq(ngrid,nlayer,nq) ! tracer tendencies (../kg/s)
-      REAL pdq_corrdyn(ngrid,nlayer,nq) ! tracer mass mixing ratio tendencies in order to correct
-                                        ! mass non-conservations in dynamics
       REAL,INTENT(out) :: pdpsrf(ngrid) ! surface pressure tendency (Pa/s)
 
@@ -451 +447 @@
       real colden(ngrid,nq)     ! vertical column of tracers
       real mass(nq)             ! global mass of tracers (g)
-      real, allocatable, save :: mass_predyn(:) ! tracer mass in the whole atmosphere
-                                                ! before dynamics
-      real mass_postdyn(nq) ! tracer mass in the whole atmosphere after dynamics
-
-!$OMP THREADPRIVATE(mass_predyn)
-
-      real masscorrfac(nq) ! mass correction factor for non-conservations in dynamics
       REAL mtot(ngrid)          ! Total mass of water vapor (kg/m2)
       REAL mstormdtot(ngrid)    ! Total mass of stormdust tracer (kg/m2)
@@ -826 +815 @@
      &                              presnivs,pseudoalt,mlayer)
 #endif
-
-         allocate(mass_predyn(nq))
-
       ENDIF        !  (end of "if firstcall")
 
@@ -892 +878 @@
       zplay(:,:) = pplay(:,:)
       ps(:) = pplev(:,1)
-
-c     Correct mass non-conservation in dynamics
-c     ~~~~~~~~~~~~~~~~~~
-      call compute_global_mass(pq, ngrid, nlayer, nq, zplev,
-     &                         mass_postdyn)
-
-      if (.not.firstcall) then ! At first call, mass_predyn does not exist
-         do iq = 1, nq
-            if (noms(iq) .ne. "dust_mass" .and.
-     &          noms(iq) .ne. "dust_number" .and.
-     &          noms(iq) .ne. "ccn_mass" .and.
-     &          noms(iq) .ne. "ccn_number" .and.
-     &          noms(iq) .ne. "ccnco2_mass" .and.
-     &          noms(iq) .ne. "ccnco2_number" .and.
-     &          noms(iq) .ne. "stormdust_mass" .and.
-     &          noms(iq) .ne. "stormdust_number" .and.
-     &          noms(iq) .ne. "topdust_mass" .and.
-     &          noms(iq) .ne. "topdust_number") then
-               masscorrfac(iq) = mass_predyn(iq)/mass_postdyn(iq)
-            endif
-         enddo
-
-         pq_corrdyn = pq
-         do ig = 1, ngrid
-            do l = 1, min(nlayer, 40) ! Correction only applied to the lowest 40 layers
-               do iq = 1, nq
-                  if (noms(iq) .ne. "dust_mass" .and.
-     &                noms(iq) .ne. "dust_number" .and.
-     &                noms(iq) .ne. "ccn_mass" .and.
-     &                noms(iq) .ne. "ccn_number" .and.
-     &                noms(iq) .ne. "ccnco2_mass" .and.
-     &                noms(iq) .ne. "ccnco2_number" .and.
-     &                noms(iq) .ne. "stormdust_mass" .and.
-     &                noms(iq) .ne. "stormdust_number" .and.
-     &                noms(iq) .ne. "topdust_mass" .and.
-     &                noms(iq) .ne. "topdust_number" .and.
-     &                iq .ne. igcm_co2) then
-                     ! Correcting species other than CO2
-                     pq_corrdyn(ig, l, iq) =
-     &                    pq(ig, l, iq) * masscorrfac(iq)
-                     ! Adjust CO2 to compensate for mass creation/loss
-                     pq_corrdyn(ig, l, igcm_co2) = 
-     &                    pq_corrdyn(ig, l, igcm_co2) -
-     &                    pq(ig, l, iq) * (masscorrfac(iq) - 1.)
-                  endif
-               enddo
-            enddo
-         enddo
-
-         pdq_corrdyn = (pq_corrdyn - pq)/ptimestep ! Mass fixer pdq
-         pdq = pdq + pdq_corrdyn
-      endif
 
 !#ifndef MESOSCALE
@@ -2634 +2568 @@
         ENDDO
       ENDDO
-
-      call compute_global_mass(zq, ngrid, nlayer, nq, zplev, 
-     &                         mass_predyn)
 
       ! Density
@@ -4199 +4130 @@
       END SUBROUTINE physiq
 
-      subroutine compute_global_mass(zq, ngrid, nlayer, nq, 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
-
-         real, intent(in)    :: zq(ngrid, nlayer, nq)
-         integer, intent(in) :: ngrid
-         integer, intent(in) :: nlayer
-         integer, intent(in) :: nq
-         real, intent(in)    :: zplev(ngrid, nlayer + 1)
-         real, intent(out)   :: mass(nq)
-         real    colden(ngrid, nq)
-         integer iq, ig, l
-
-         do iq = 1, nq
-            if (noms(iq) .ne. "dust_mass" .and.
-     &          noms(iq) .ne. "dust_number" .and.
-     &          noms(iq) .ne. "ccn_mass" .and.
-     &          noms(iq) .ne. "ccn_number" .and.
-     &          noms(iq) .ne. "ccnco2_mass" .and.
-     &          noms(iq) .ne. "ccnco2_number" .and.
-     &          noms(iq) .ne. "stormdust_mass" .and.
-     &          noms(iq) .ne. "stormdust_number" .and.
-     &          noms(iq) .ne. "topdust_mass" .and.
-     &          noms(iq) .ne. "topdust_number") then
-
-               do ig = 1, ngrid
-                   colden(ig, iq) = 0.
-               end do
-               do l = 1, nlayer
-                  do ig = 1, ngrid
-                     colden(ig, iq) = colden(ig, iq) + zq(ig, l, iq)*
-     &                                (zplev(ig, l) - zplev(ig, l+1))*
-     &                                6.022e22/(mmol(iq)*g)
-                  end do
-               end do
-
-               call planetwide_sumval(colden(:, iq)/6.022e23
-     &                    *mmol(iq)*1.e4*cell_area, mass(iq))
-
-            end if
-         end do
-
-      end subroutine
-
       END MODULE physiq_mod