Ignore:
Timestamp:
Mar 27, 2026, 12:04:04 PM (7 weeks ago)
Author:
yluo
Message:

Mars PCM:
Add a mass-conservation fixer for moldiff_MPF to correct tracer mass changes introduced by molecular diffusion.
The tendencies above the 1 Pa level remain unchanged. Below 1 Pa, tendencies are adjusted such that the column-integrated mass of each species after moldiff_MPF matches that before moldiff_MPF.
Special treatment is applied to H, O, and N due to their very low abundances below 1 Pa: the mixing ratios of H2, O2, and N2 below 1 Pa are adjusted to conserve the total abundances of H+H2, O+O2, and N+N2 in the whole column, respectively.
GCMGRID_P2 is revised to account for the fact that the mixing ratios of all diffused species do not sum to one.
At each grid point, the CO2 tendency is set to the negative sum of the tendencies of all other diffused species to ensure molecular diffusion does not create/destroy mass.
Add a flag call_mass_fixer_moldiff_MPF to enable/disable the above corrections (default: .false.).
YCL

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90

    r4035 r4159  
    44  real*8,parameter :: tdiffmin=5d0
    55  real*8,parameter :: dzres=2d0 ! grid resolution (km) for diffusion
     6  logical, save    :: call_mass_fixer_moldiff_MPF = .false. ! flag for correcting tracer mass non-conservations in moldiff_MPF
     7
     8!$OMP THREADPRIVATE(call_mass_fixer_moldiff_MPF)
    69
    710CONTAINS
    … …  
    5154
    5255      integer,dimension(nq) :: indic_diff
    53       integer ig,iq,nz,l,k,n,nn,p,ij0
     56      integer ig,iq,nz,l,k,n,nn,nn_o,nn_o2,nn_h,nn_h2,nn_n,nn_n2,nn_co2,p,ij0
    5457      integer istep,il,gcn,ntime,nlraf
    5558      real*8 masse
    … …  
    8083      real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie,alphaT
    8184!$OMP THREADPRIVATE(wi,Wad,Uthermal,Lambdaexo,Hspecie,alphaT)
    82       real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
    83 !$OMP THREADPRIVATE(Mtot1,Mtot2,Mraf1,Mraf2)
     85      real*8, dimension(:), allocatable, save :: Mtot1 ! mass of a species in whole column before molecular diffusion
     86      real*8, dimension(:), allocatable, save :: Mtot2 ! mass of a species in whole column after molecular diffusion which violates mass conservation
     87      real*8, dimension(:), allocatable, save :: Mtot3 ! mass of a species in whole column after mass non-conservation correction to species other than H, N, O
     88      real*8, dimension(:), allocatable, save :: Mhi1 ! mass of a species above 1 Pa level before molecular diffusion
     89      real*8, dimension(:), allocatable, save :: Mhi2 ! mass of a species above 1 Pa level after molecular diffusion before mass correction
     90!$OMP THREADPRIVATE(Mtot1, Mtot2, Mtot3, Mhi1, Mhi2)
     91      ! real*8, dimension(:), allocatable, save :: Mraf1, Mraf2
     92! !$OMP THREADPRIVATE(Mraf1, Mraf2)
    8493      integer,parameter :: ListeDiffNb=16
    8594      character(len=20),dimension(ListeDiffNb) :: ListeDiff
    … …  
    97106! vertical index limit for the molecular diffusion
    98107      integer,save :: il0 
    99 !$OMP THREADPRIVATE(i_h2,i_h,i_d,i_hd,il0)
     108! vertical index below which qnew is scaled to conserve mass
     109      integer,save :: il1
     110!$OMP THREADPRIVATE(i_h2,i_h,i_d,i_hd,il0,il1)
    100111
    101112! Tracer indexes in the GCM:
    … …  
    125136      integer ierr,compteur
    126137
     138      real*8, dimension(:), allocatable, save  :: masscorrfac  ! mass correction factors applied below 1 Pa level to conserve mass in whole column for species other than H, N, O
     139      real*8                                   :: masscorrfac1 ! mass correction factor for O2, N2, H2 below 1 Pa level to correct mass non-conservation in H, N, O
     140      real*8                                   :: sum_tend
     141
    127142      logical,save :: firstcall=.true.
    128143!      real abfac(ncompdiff)
    … …  
    130145      real,save :: step
    131146
    132 !$OMP THREADPRIVATE(ncompdiff,gcmind,firstcall,dij,step)     
     147!$OMP THREADPRIVATE(ncompdiff, gcmind, masscorrfac, firstcall, dij, step)
    133148
    134149! Initializations at first call
    … …  
    178193        print*,'number of diffused species:',ncompdiff
    179194        allocate(gcmind(ncompdiff))
     195        if (call_mass_fixer_moldiff_MPF) then
     196            allocate(masscorrfac(ncompdiff))
     197        endif
    180198
    181199! Store gcm indexes in gcmind
    … …  
    202220        enddo
    203221        il0=il0+1
     222            if (call_mass_fixer_moldiff_MPF) then
     223                ! find layer index below which scaling is applied to qnew in order to conserve mass
     224                do l = 1, nlayer
     225                    if (pplay(1, l) .gt. 1.) then ! threshold: 1 Pa
     226                        il1 = l
     227                    endif
     228                enddo
     229                il1 = il1 + 1
     230            endif ! if (call_mass_fixer_moldiff_MPF)
    204231        endif ! (is_master)
    205232        CALL bcast(il0)
     233        if (call_mass_fixer_moldiff_MPF) CALL bcast(il1)
    206234        print*,'vertical index for diffusion',il0,pplay(1,il0)
    207235
    … …  
    225253        allocate(lambdaexo(ncompdiff))
    226254        allocate(Hspecie(ncompdiff))
    227         allocate(Mtot1(ncompdiff))
    228         allocate(Mtot2(ncompdiff))
    229         allocate(Mraf1(ncompdiff))
    230         allocate(Mraf2(ncompdiff))
     255        ! allocate(Mraf1(ncompdiff))
     256        ! allocate(Mraf2(ncompdiff))
     257        if (call_mass_fixer_moldiff_MPF) then
     258            allocate(Mtot1(ncompdiff))
     259            allocate(Mtot2(ncompdiff))
     260            allocate(Mtot3(ncompdiff))
     261            allocate(Mhi1(ncompdiff))
     262            allocate(Mhi2(ncompdiff))
     263        endif
    231264
    232265        firstcall= .false.
    … …  
    317350! e.g. OH, O(1D)
    318351
    319         Mtot1(1:ncompdiff)=0d0
    320 
    321         do l=il0,nlayer
    322         do nn=1,ncompdiff
    323         Mtot1(nn)=Mtot1(nn)+1d0/g*qq(l,nn)*                             &
    324      &  (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
    325         enddo
    326         enddo
     352        if (call_mass_fixer_moldiff_MPF) then
     353            ! Compute vertically integrated mass for each diffused chemical species before diffusion
     354            ! From surface to TOA:
     355            Mtot1(:) = 0d0
     356
     357            do l = 1, nlayer
     358                Mtot1(:) = Mtot1(:) + 1d0 / g * qq(l, :) * &
     359                           (dble(pplev(ig, l)) - dble(pplev(ig, l + 1)))
     360            enddo
     361
     362            ! Above 1 Pa:
     363            Mhi1(:) = 0d0
     364
     365            do l = il1 + 1, nlayer
     366                Mhi1(:) = Mhi1(:) + 1d0 / g * qq(l, :) * &
     367                          (dble(pplev(ig, l)) - dble(pplev(ig, l + 1)))
     368            end do
     369        endif ! if (call_mass_fixer_moldiff_MPF)
    327370
    328371        Zmin=zz(il0)
    … …  
    398441! Total mass computed on the altitude grid
    399442
    400         Mraf1(1:ncompdiff)=0d0
    401         do nn=1,ncompdiff
    402         do l=1,nlraf
    403         Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf
    404         enddo
    405         enddo
     443!       Mraf1(1:ncompdiff)=0d0
     444!       do nn=1,ncompdiff
     445!       do l=1,nlraf
     446!       Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf
     447!       enddo
     448!       enddo
    406449
    407450! Reupdate values for mass conservation
    … …  
    637680! Compute the total mass of each species to check mass conservation     
    638681
    639         Mraf2(1:ncompdiff)=0d0
    640         do nn=1,ncompdiff
    641         do l=1,nlraf
    642         Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf
    643         enddo
    644         enddo
     682!        Mraf2(1:ncompdiff)=0d0
     683!        do nn=1,ncompdiff
     684!        do l=1,nlraf
     685!        Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf
     686!        enddo
     687!        enddo
    645688
    646689!        print*,'Mraf',Mraf2
    … …  
    679722! Compute total mass of each specie on the GCM vertical grid
    680723
    681         Mtot2(1:ncompdiff)=0d0
    682 
    683         do l=il0,nlayer
    684         do nn=1,ncompdiff
    685         Mtot2(nn)=Mtot2(nn)+1d0/g*qnew(l,nn)*                           &
    686      &  (dble(pplev(ig,l))-dble(pplev(ig,l+1)))
    687         enddo
    688         enddo
     724        if (call_mass_fixer_moldiff_MPF) then
     725            ! Compute vertically integrated mass for each diffused chemical species immediately after diffusion (before correction)
     726            ! From surface to TOA:
     727            Mtot2(:) = 0d0
     728
     729            do l = 1, nlayer
     730                Mtot2(:) = Mtot2(:) + 1d0 / g * qnew(l, :) * &
     731                           (dble(pplev(ig, l)) - dble(pplev(ig, l + 1)))
     732            enddo
     733
     734            ! Above 1 Pa:
     735            Mhi2(:) = 0d0
     736
     737            do l = il1 + 1, nlayer
     738                Mhi2(:) = Mhi2(:) + 1d0 / g * qnew(l, :) * &
     739                          (dble(pplev(ig, l)) - dble(pplev(ig, l + 1)))
     740            enddo
     741
     742            ! For species other than H, N, and O, scale their mixing ratios below 1 Pa to conserve mass from surface to TOA
     743            masscorrfac(:) = (Mtot1(:) - Mhi2(:)) / (Mtot2(:) - Mhi2(:))
     744
     745            do nn = 1, ncompdiff
     746                if (noms(gcmind(nn)) .eq. 'h' .or. noms(gcmind(nn)) .eq. 'n' .or. noms(gcmind(nn)) .eq. 'o') masscorrfac(nn) = 1d0
     747            enddo
     748
     749            do nn = 1, ncompdiff
     750                qnew(1:il1, nn) = qnew(1:il1, nn) * masscorrfac(nn)
     751            enddo
     752
     753            ! Compute vertically integrated mass for each diffused chemical species after diffusion and correction (to species other
     754            ! than H, N, and O)
     755            ! For species other than H, N, and O, Mtot3 should be equal to Mtot1, as correction has been applied
     756            ! For H, N, and O, Mtot3 should be equal to Mtot2, as correction has not been applied
     757            Mtot3(:) = 0d0
     758
     759            do l = 1, nlayer
     760                Mtot3(:) = Mtot3(:) + 1d0 / g * qnew(l, :) * &
     761                           (dble(pplev(ig, l)) - dble(pplev(ig, l + 1)))
     762            enddo
     763
     764            ! Treat H, O, and N separately by scaling the mixing ratios of H2, O2, and N2 below 1 Pa
     765            ! Cannot treat them like others because their abundances below 1 Pa are too small
     766            nn_h = 0
     767            nn_h2 = 0
     768            nn_o = 0
     769            nn_o2 = 0
     770            nn_n = 0
     771            nn_n2 = 0
     772
     773            do nn = 1, ncompdiff
     774                if (noms(gcmind(nn)) .eq. 'h')  nn_h  = nn
     775                if (noms(gcmind(nn)) .eq. 'h2') nn_h2 = nn
     776                if (noms(gcmind(nn)) .eq. 'o')  nn_o  = nn
     777                if (noms(gcmind(nn)) .eq. 'o2') nn_o2 = nn
     778                if (noms(gcmind(nn)) .eq. 'n')  nn_n  = nn
     779                if (noms(gcmind(nn)) .eq. 'n2') nn_n2 = nn
     780            end do
     781
     782            ! Correct O2 below 1 Pa to compensate for mass non-conservation of O
     783            if (nn_o2 .gt. 0 .and. nn_o .gt. 0) then
     784                masscorrfac1 = (Mtot3(nn_o2) + Mtot1(nn_o) - Mtot3(nn_o) - Mhi2(nn_o2))/(Mtot3(nn_o2) - Mhi2(nn_o2))
     785                do l = 1, il1
     786                    qnew(l, nn_o2) = qnew(l, nn_o2)*masscorrfac1
     787                end do
     788            endif
     789
     790            ! Correct H2 below 1 Pa to compensate for mass non-conservation of H
     791            if (nn_h2 .gt. 0 .and. nn_h .gt. 0) then
     792                masscorrfac1 = (Mtot3(nn_h2) + Mtot1(nn_h) - Mtot3(nn_h) - Mhi2(nn_h2))/(Mtot3(nn_h2) - Mhi2(nn_h2))
     793                do l = 1, il1
     794                    qnew(l, nn_h2) = qnew(l, nn_h2)*masscorrfac1
     795                end do
     796            endif
     797
     798            ! Correct N2 below 1 Pa to compensate for mass non-conservation of N
     799            if (nn_n2 .gt. 0 .and. nn_n .gt. 0) then
     800                masscorrfac1 = (Mtot3(nn_n2) + Mtot1(nn_n) - Mtot3(nn_n) - Mhi2(nn_n2))/(Mtot3(nn_n2) - Mhi2(nn_n2))
     801                do l = 1, il1
     802                    qnew(l, nn_n2) = qnew(l, nn_n2)*masscorrfac1
     803                end do
     804            endif
     805
     806        endif ! if (call_mass_fixer_moldiff_MPF)
    689807
    690808! Check mass  conservation of each specie on column
    … …  
    696814! Compute the diffusion trends du to diffusion
    697815
     816        if (call_mass_fixer_moldiff_MPF) then
     817            ! To ensure the sum of the tendencies of all species is zero, CO2 is used as the compensating species
     818            nn_co2 = 0
     819            do nn = 1, ncompdiff
     820                if (noms(gcmind(nn)) .eq. 'co2') then
     821                    nn_co2 = nn
     822                    exit
     823                endif
     824            enddo
     825
     826            if (nn_co2 .eq. 0) stop "Tracer CO2 does not exist or does not diffuse"
     827
     828            do l = 1, nlayer
     829                sum_tend = 0d0
     830
     831                do nn = 1, ncompdiff
     832                    if (nn .ne. nn_co2) then
     833                        pdqdiff(ig, l, gcmind(nn)) = (qnew(l, nn) - qq(l, nn)) / ptimestep
     834                        sum_tend = sum_tend + pdqdiff(ig, l, gcmind(nn))
     835                    endif
     836                enddo
     837
     838                ! CO2 tendency = 0 - sum(tendencies of all other tracers)
     839                pdqdiff(ig, l, gcmind(nn_co2)) = -sum_tend
     840            enddo ! l
     841        else
     842        ! The traditional algorithm: compute pdqdiff for all species (which may not add up to zero)
    698843        do l=1,nlayer
    699844        do nn=1,ncompdiff
    … …  
    701846        enddo
    702847        enddo
     848        endif ! if (call_mass_fixer_moldiff_MPF)
    703849
    704850! deallocation des tableaux
    … …  
    15691715        tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP
    15701716        rhonew(il)=sum(rhoknew(il,:))
     1717        if (call_mass_fixer_moldiff_MPF) then
     1718            do nn = 1, nq
     1719                ! To account for the fact that the mixing ratios of all diffused species do not sum to 1
     1720                qnew(il, nn) = rhoknew(il, nn) / rhonew(il) * sum(qq(il, :))
     1721            enddo
     1722        else
    15711723        do nn=1,nq
    15721724        qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
    15731725        enddo
     1726        endif ! if (call_mass_fixer_moldiff_MPF)
    15741727
    15751728        else ! pp < P(nl) need to extrapolate density of each specie
    … …  
    16071760        enddo
    16081761        rhonew(il)=sum(rhoknew(il,:))
     1762        if (call_mass_fixer_moldiff_MPF) then
     1763            do nn = 1, nq
     1764                ! To account for the fact that the mixing ratios of all diffused species do not sum to 1
     1765                qnew(il, nn) = rhoknew(il, nn) / rhonew(il) * sum(qq(il, :))
     1766            enddo
     1767        else
    16091768        do nn=1,nq
    16101769        qnew(il,nn)=rhoknew(il,nn)/rhonew(il)
    16111770        enddo
     1771        endif ! if (call_mass_fixer_moldiff_MPF)
    16121772        tnew(il)=T(nl)
    16131773
Note: See TracChangeset for help on using the changeset viewer.