Index: /trunk/LMDZ.MARS/changelog.txt =================================================================== --- /trunk/LMDZ.MARS/changelog.txt (revision 4158) +++ /trunk/LMDZ.MARS/changelog.txt (revision 4159) @@ -5130,2 +5130,10 @@ == 26/03/2026 == YCL Add a flag call_mass_fixer_dyn for correcting tracer mass non-conservation in dynamics to reference def files in deftank. + +== 27/03/2026 == YCL +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.). Index: /trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90 =================================================================== --- /trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90 (revision 4158) +++ /trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90 (revision 4159) @@ -4,4 +4,7 @@ real*8,parameter :: tdiffmin=5d0 real*8,parameter :: dzres=2d0 ! grid resolution (km) for diffusion + logical, save :: call_mass_fixer_moldiff_MPF = .false. ! flag for correcting tracer mass non-conservations in moldiff_MPF + +!$OMP THREADPRIVATE(call_mass_fixer_moldiff_MPF) CONTAINS @@ -51,5 +54,5 @@ integer,dimension(nq) :: indic_diff - integer ig,iq,nz,l,k,n,nn,p,ij0 + integer ig,iq,nz,l,k,n,nn,nn_o,nn_o2,nn_h,nn_h2,nn_n,nn_n2,nn_co2,p,ij0 integer istep,il,gcn,ntime,nlraf real*8 masse @@ -80,6 +83,12 @@ real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie,alphaT !$OMP THREADPRIVATE(wi,Wad,Uthermal,Lambdaexo,Hspecie,alphaT) - real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2 -!$OMP THREADPRIVATE(Mtot1,Mtot2,Mraf1,Mraf2) + real*8, dimension(:), allocatable, save :: Mtot1 ! mass of a species in whole column before molecular diffusion + real*8, dimension(:), allocatable, save :: Mtot2 ! mass of a species in whole column after molecular diffusion which violates mass conservation + 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 + real*8, dimension(:), allocatable, save :: Mhi1 ! mass of a species above 1 Pa level before molecular diffusion + real*8, dimension(:), allocatable, save :: Mhi2 ! mass of a species above 1 Pa level after molecular diffusion before mass correction +!$OMP THREADPRIVATE(Mtot1, Mtot2, Mtot3, Mhi1, Mhi2) + ! real*8, dimension(:), allocatable, save :: Mraf1, Mraf2 +! !$OMP THREADPRIVATE(Mraf1, Mraf2) integer,parameter :: ListeDiffNb=16 character(len=20),dimension(ListeDiffNb) :: ListeDiff @@ -97,5 +106,7 @@ ! vertical index limit for the molecular diffusion integer,save :: il0 -!$OMP THREADPRIVATE(i_h2,i_h,i_d,i_hd,il0) +! vertical index below which qnew is scaled to conserve mass + integer,save :: il1 +!$OMP THREADPRIVATE(i_h2,i_h,i_d,i_hd,il0,il1) ! Tracer indexes in the GCM: @@ -125,4 +136,8 @@ integer ierr,compteur + 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 + real*8 :: masscorrfac1 ! mass correction factor for O2, N2, H2 below 1 Pa level to correct mass non-conservation in H, N, O + real*8 :: sum_tend + logical,save :: firstcall=.true. ! real abfac(ncompdiff) @@ -130,5 +145,5 @@ real,save :: step -!$OMP THREADPRIVATE(ncompdiff,gcmind,firstcall,dij,step) +!$OMP THREADPRIVATE(ncompdiff, gcmind, masscorrfac, firstcall, dij, step) ! Initializations at first call @@ -178,4 +193,7 @@ print*,'number of diffused species:',ncompdiff allocate(gcmind(ncompdiff)) + if (call_mass_fixer_moldiff_MPF) then + allocate(masscorrfac(ncompdiff)) + endif ! Store gcm indexes in gcmind @@ -202,6 +220,16 @@ enddo il0=il0+1 + if (call_mass_fixer_moldiff_MPF) then + ! find layer index below which scaling is applied to qnew in order to conserve mass + do l = 1, nlayer + if (pplay(1, l) .gt. 1.) then ! threshold: 1 Pa + il1 = l + endif + enddo + il1 = il1 + 1 + endif ! if (call_mass_fixer_moldiff_MPF) endif ! (is_master) CALL bcast(il0) + if (call_mass_fixer_moldiff_MPF) CALL bcast(il1) print*,'vertical index for diffusion',il0,pplay(1,il0) @@ -225,8 +253,13 @@ allocate(lambdaexo(ncompdiff)) allocate(Hspecie(ncompdiff)) - allocate(Mtot1(ncompdiff)) - allocate(Mtot2(ncompdiff)) - allocate(Mraf1(ncompdiff)) - allocate(Mraf2(ncompdiff)) + ! allocate(Mraf1(ncompdiff)) + ! allocate(Mraf2(ncompdiff)) + if (call_mass_fixer_moldiff_MPF) then + allocate(Mtot1(ncompdiff)) + allocate(Mtot2(ncompdiff)) + allocate(Mtot3(ncompdiff)) + allocate(Mhi1(ncompdiff)) + allocate(Mhi2(ncompdiff)) + endif firstcall= .false. @@ -317,12 +350,22 @@ ! e.g. OH, O(1D) - Mtot1(1:ncompdiff)=0d0 - - do l=il0,nlayer - do nn=1,ncompdiff - Mtot1(nn)=Mtot1(nn)+1d0/g*qq(l,nn)* & - & (dble(pplev(ig,l))-dble(pplev(ig,l+1))) - enddo - enddo + if (call_mass_fixer_moldiff_MPF) then + ! Compute vertically integrated mass for each diffused chemical species before diffusion + ! From surface to TOA: + Mtot1(:) = 0d0 + + do l = 1, nlayer + Mtot1(:) = Mtot1(:) + 1d0 / g * qq(l, :) * & + (dble(pplev(ig, l)) - dble(pplev(ig, l + 1))) + enddo + + ! Above 1 Pa: + Mhi1(:) = 0d0 + + do l = il1 + 1, nlayer + Mhi1(:) = Mhi1(:) + 1d0 / g * qq(l, :) * & + (dble(pplev(ig, l)) - dble(pplev(ig, l + 1))) + end do + endif ! if (call_mass_fixer_moldiff_MPF) Zmin=zz(il0) @@ -398,10 +441,10 @@ ! Total mass computed on the altitude grid - Mraf1(1:ncompdiff)=0d0 - do nn=1,ncompdiff - do l=1,nlraf - Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf - enddo - enddo +! Mraf1(1:ncompdiff)=0d0 +! do nn=1,ncompdiff +! do l=1,nlraf +! Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf +! enddo +! enddo ! Reupdate values for mass conservation @@ -637,10 +680,10 @@ ! Compute the total mass of each species to check mass conservation - Mraf2(1:ncompdiff)=0d0 - do nn=1,ncompdiff - do l=1,nlraf - Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf - enddo - enddo +! Mraf2(1:ncompdiff)=0d0 +! do nn=1,ncompdiff +! do l=1,nlraf +! Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf +! enddo +! enddo ! print*,'Mraf',Mraf2 @@ -679,12 +722,87 @@ ! Compute total mass of each specie on the GCM vertical grid - Mtot2(1:ncompdiff)=0d0 - - do l=il0,nlayer - do nn=1,ncompdiff - Mtot2(nn)=Mtot2(nn)+1d0/g*qnew(l,nn)* & - & (dble(pplev(ig,l))-dble(pplev(ig,l+1))) - enddo - enddo + if (call_mass_fixer_moldiff_MPF) then + ! Compute vertically integrated mass for each diffused chemical species immediately after diffusion (before correction) + ! From surface to TOA: + Mtot2(:) = 0d0 + + do l = 1, nlayer + Mtot2(:) = Mtot2(:) + 1d0 / g * qnew(l, :) * & + (dble(pplev(ig, l)) - dble(pplev(ig, l + 1))) + enddo + + ! Above 1 Pa: + Mhi2(:) = 0d0 + + do l = il1 + 1, nlayer + Mhi2(:) = Mhi2(:) + 1d0 / g * qnew(l, :) * & + (dble(pplev(ig, l)) - dble(pplev(ig, l + 1))) + enddo + + ! For species other than H, N, and O, scale their mixing ratios below 1 Pa to conserve mass from surface to TOA + masscorrfac(:) = (Mtot1(:) - Mhi2(:)) / (Mtot2(:) - Mhi2(:)) + + do nn = 1, ncompdiff + if (noms(gcmind(nn)) .eq. 'h' .or. noms(gcmind(nn)) .eq. 'n' .or. noms(gcmind(nn)) .eq. 'o') masscorrfac(nn) = 1d0 + enddo + + do nn = 1, ncompdiff + qnew(1:il1, nn) = qnew(1:il1, nn) * masscorrfac(nn) + enddo + + ! Compute vertically integrated mass for each diffused chemical species after diffusion and correction (to species other + ! than H, N, and O) + ! For species other than H, N, and O, Mtot3 should be equal to Mtot1, as correction has been applied + ! For H, N, and O, Mtot3 should be equal to Mtot2, as correction has not been applied + Mtot3(:) = 0d0 + + do l = 1, nlayer + Mtot3(:) = Mtot3(:) + 1d0 / g * qnew(l, :) * & + (dble(pplev(ig, l)) - dble(pplev(ig, l + 1))) + enddo + + ! Treat H, O, and N separately by scaling the mixing ratios of H2, O2, and N2 below 1 Pa + ! Cannot treat them like others because their abundances below 1 Pa are too small + nn_h = 0 + nn_h2 = 0 + nn_o = 0 + nn_o2 = 0 + nn_n = 0 + nn_n2 = 0 + + do nn = 1, ncompdiff + if (noms(gcmind(nn)) .eq. 'h') nn_h = nn + if (noms(gcmind(nn)) .eq. 'h2') nn_h2 = nn + if (noms(gcmind(nn)) .eq. 'o') nn_o = nn + if (noms(gcmind(nn)) .eq. 'o2') nn_o2 = nn + if (noms(gcmind(nn)) .eq. 'n') nn_n = nn + if (noms(gcmind(nn)) .eq. 'n2') nn_n2 = nn + end do + + ! Correct O2 below 1 Pa to compensate for mass non-conservation of O + if (nn_o2 .gt. 0 .and. nn_o .gt. 0) then + masscorrfac1 = (Mtot3(nn_o2) + Mtot1(nn_o) - Mtot3(nn_o) - Mhi2(nn_o2))/(Mtot3(nn_o2) - Mhi2(nn_o2)) + do l = 1, il1 + qnew(l, nn_o2) = qnew(l, nn_o2)*masscorrfac1 + end do + endif + + ! Correct H2 below 1 Pa to compensate for mass non-conservation of H + if (nn_h2 .gt. 0 .and. nn_h .gt. 0) then + masscorrfac1 = (Mtot3(nn_h2) + Mtot1(nn_h) - Mtot3(nn_h) - Mhi2(nn_h2))/(Mtot3(nn_h2) - Mhi2(nn_h2)) + do l = 1, il1 + qnew(l, nn_h2) = qnew(l, nn_h2)*masscorrfac1 + end do + endif + + ! Correct N2 below 1 Pa to compensate for mass non-conservation of N + if (nn_n2 .gt. 0 .and. nn_n .gt. 0) then + masscorrfac1 = (Mtot3(nn_n2) + Mtot1(nn_n) - Mtot3(nn_n) - Mhi2(nn_n2))/(Mtot3(nn_n2) - Mhi2(nn_n2)) + do l = 1, il1 + qnew(l, nn_n2) = qnew(l, nn_n2)*masscorrfac1 + end do + endif + + endif ! if (call_mass_fixer_moldiff_MPF) ! Check mass conservation of each specie on column @@ -696,4 +814,31 @@ ! Compute the diffusion trends du to diffusion + if (call_mass_fixer_moldiff_MPF) then + ! To ensure the sum of the tendencies of all species is zero, CO2 is used as the compensating species + nn_co2 = 0 + do nn = 1, ncompdiff + if (noms(gcmind(nn)) .eq. 'co2') then + nn_co2 = nn + exit + endif + enddo + + if (nn_co2 .eq. 0) stop "Tracer CO2 does not exist or does not diffuse" + + do l = 1, nlayer + sum_tend = 0d0 + + do nn = 1, ncompdiff + if (nn .ne. nn_co2) then + pdqdiff(ig, l, gcmind(nn)) = (qnew(l, nn) - qq(l, nn)) / ptimestep + sum_tend = sum_tend + pdqdiff(ig, l, gcmind(nn)) + endif + enddo + + ! CO2 tendency = 0 - sum(tendencies of all other tracers) + pdqdiff(ig, l, gcmind(nn_co2)) = -sum_tend + enddo ! l + else + ! The traditional algorithm: compute pdqdiff for all species (which may not add up to zero) do l=1,nlayer do nn=1,ncompdiff @@ -701,4 +846,5 @@ enddo enddo + endif ! if (call_mass_fixer_moldiff_MPF) ! deallocation des tableaux @@ -1569,7 +1715,14 @@ tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP rhonew(il)=sum(rhoknew(il,:)) + if (call_mass_fixer_moldiff_MPF) then + do nn = 1, nq + ! To account for the fact that the mixing ratios of all diffused species do not sum to 1 + qnew(il, nn) = rhoknew(il, nn) / rhonew(il) * sum(qq(il, :)) + enddo + else do nn=1,nq qnew(il,nn)=rhoknew(il,nn)/rhonew(il) enddo + endif ! if (call_mass_fixer_moldiff_MPF) else ! pp < P(nl) need to extrapolate density of each specie @@ -1607,7 +1760,14 @@ enddo rhonew(il)=sum(rhoknew(il,:)) + if (call_mass_fixer_moldiff_MPF) then + do nn = 1, nq + ! To account for the fact that the mixing ratios of all diffused species do not sum to 1 + qnew(il, nn) = rhoknew(il, nn) / rhonew(il) * sum(qq(il, :)) + enddo + else do nn=1,nq qnew(il,nn)=rhoknew(il,nn)/rhonew(il) enddo + endif ! if (call_mass_fixer_moldiff_MPF) tnew(il)=T(nl) Index: /trunk/LMDZ.MARS/libf/phymars/conf_phys.F =================================================================== --- /trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 4158) +++ /trunk/LMDZ.MARS/libf/phymars/conf_phys.F (revision 4159) @@ -51,4 +51,5 @@ use lwxn_mod, only: lwxn_linear, lwxn_alphan, lwxn_ncouche use tracer_mass_fixer_dyn_mod, only: call_mass_fixer_dyn + use moldiff_MPF_mod, only: call_mass_fixer_moldiff_MPF use callkeys_mod, only: startphy_file, activice, activeco2ice, & calladj, callatke, callcond, @@ -1255,4 +1256,12 @@ write(*,*) "call_mass_fixer_dyn = ", call_mass_fixer_dyn +! Tracer mass fixer for moldiff_MPF + + write(*,*) "call tracer mass fixer for moldiff_MPF?" + call getin_p("call_mass_fixer_moldiff_MPF", + & call_mass_fixer_moldiff_MPF) + write(*,*) "call_mass_fixer_moldiff_MPF = ", + & call_mass_fixer_moldiff_MPF + ! Test of incompatibility: ! if photochem is used, then water should be used too