Changeset 4159
- Timestamp:
- Mar 27, 2026, 12:04:04 PM (6 weeks ago)
- Location:
- trunk/LMDZ.MARS
- Files:
-
- 3 edited
-
changelog.txt (modified) (1 diff)
-
libf/aeronomars/moldiff_MPF.F90 (modified) (17 diffs)
-
libf/phymars/conf_phys.F (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/changelog.txt
r4156 r4159 5130 5130 == 26/03/2026 == YCL 5131 5131 Add a flag call_mass_fixer_dyn for correcting tracer mass non-conservation in dynamics to reference def files in deftank. 5132 5133 == 27/03/2026 == YCL 5134 Add a mass-conservation fixer for moldiff_MPF to correct tracer mass changes introduced by molecular diffusion. 5135 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. 5136 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. 5137 GCMGRID_P2 is revised to account for the fact that the mixing ratios of all diffused species do not sum to one. 5138 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. 5139 Add a flag call_mass_fixer_moldiff_MPF to enable/disable the above corrections (default: .false.). -
trunk/LMDZ.MARS/libf/aeronomars/moldiff_MPF.F90
r4035 r4159 4 4 real*8,parameter :: tdiffmin=5d0 5 5 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) 6 9 7 10 CONTAINS … … 51 54 52 55 integer,dimension(nq) :: indic_diff 53 integer ig,iq,nz,l,k,n,nn, p,ij056 integer ig,iq,nz,l,k,n,nn,nn_o,nn_o2,nn_h,nn_h2,nn_n,nn_n2,nn_co2,p,ij0 54 57 integer istep,il,gcn,ntime,nlraf 55 58 real*8 masse … … 80 83 real*8,dimension(:),allocatable,save :: wi,Wad,Uthermal,Lambdaexo,Hspecie,alphaT 81 84 !$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) 84 93 integer,parameter :: ListeDiffNb=16 85 94 character(len=20),dimension(ListeDiffNb) :: ListeDiff … … 97 106 ! vertical index limit for the molecular diffusion 98 107 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) 100 111 101 112 ! Tracer indexes in the GCM: … … 125 136 integer ierr,compteur 126 137 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 127 142 logical,save :: firstcall=.true. 128 143 ! real abfac(ncompdiff) … … 130 145 real,save :: step 131 146 132 !$OMP THREADPRIVATE(ncompdiff, gcmind,firstcall,dij,step)147 !$OMP THREADPRIVATE(ncompdiff, gcmind, masscorrfac, firstcall, dij, step) 133 148 134 149 ! Initializations at first call … … 178 193 print*,'number of diffused species:',ncompdiff 179 194 allocate(gcmind(ncompdiff)) 195 if (call_mass_fixer_moldiff_MPF) then 196 allocate(masscorrfac(ncompdiff)) 197 endif 180 198 181 199 ! Store gcm indexes in gcmind … … 202 220 enddo 203 221 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) 204 231 endif ! (is_master) 205 232 CALL bcast(il0) 233 if (call_mass_fixer_moldiff_MPF) CALL bcast(il1) 206 234 print*,'vertical index for diffusion',il0,pplay(1,il0) 207 235 … … 225 253 allocate(lambdaexo(ncompdiff)) 226 254 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 231 264 232 265 firstcall= .false. … … 317 350 ! e.g. OH, O(1D) 318 351 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) 327 370 328 371 Zmin=zz(il0) … … 398 441 ! Total mass computed on the altitude grid 399 442 400 401 do nn=1,ncompdiff402 do l=1,nlraf403 Mraf1(nn)=Mraf1(nn)+Rrafk(l,nn)*Dzraf404 enddo405 enddo443 ! 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 406 449 407 450 ! Reupdate values for mass conservation … … 637 680 ! Compute the total mass of each species to check mass conservation 638 681 639 640 do nn=1,ncompdiff641 do l=1,nlraf642 Mraf2(nn)=Mraf2(nn)+Rrafk(l,nn)*Dzraf643 enddo644 enddo682 ! 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 645 688 646 689 ! print*,'Mraf',Mraf2 … … 679 722 ! Compute total mass of each specie on the GCM vertical grid 680 723 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) 689 807 690 808 ! Check mass conservation of each specie on column … … 696 814 ! Compute the diffusion trends du to diffusion 697 815 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) 698 843 do l=1,nlayer 699 844 do nn=1,ncompdiff … … 701 846 enddo 702 847 enddo 848 endif ! if (call_mass_fixer_moldiff_MPF) 703 849 704 850 ! deallocation des tableaux … … 1569 1715 tnew(il)=T(iP)+(T(ip+1)-T(iP))*facP 1570 1716 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 1571 1723 do nn=1,nq 1572 1724 qnew(il,nn)=rhoknew(il,nn)/rhonew(il) 1573 1725 enddo 1726 endif ! if (call_mass_fixer_moldiff_MPF) 1574 1727 1575 1728 else ! pp < P(nl) need to extrapolate density of each specie … … 1607 1760 enddo 1608 1761 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 1609 1768 do nn=1,nq 1610 1769 qnew(il,nn)=rhoknew(il,nn)/rhonew(il) 1611 1770 enddo 1771 endif ! if (call_mass_fixer_moldiff_MPF) 1612 1772 tnew(il)=T(nl) 1613 1773 -
trunk/LMDZ.MARS/libf/phymars/conf_phys.F
r4154 r4159 51 51 use lwxn_mod, only: lwxn_linear, lwxn_alphan, lwxn_ncouche 52 52 use tracer_mass_fixer_dyn_mod, only: call_mass_fixer_dyn 53 use moldiff_MPF_mod, only: call_mass_fixer_moldiff_MPF 53 54 use callkeys_mod, only: startphy_file, activice, activeco2ice, 54 55 & calladj, callatke, callcond, … … 1255 1256 write(*,*) "call_mass_fixer_dyn = ", call_mass_fixer_dyn 1256 1257 1258 ! Tracer mass fixer for moldiff_MPF 1259 1260 write(*,*) "call tracer mass fixer for moldiff_MPF?" 1261 call getin_p("call_mass_fixer_moldiff_MPF", 1262 & call_mass_fixer_moldiff_MPF) 1263 write(*,*) "call_mass_fixer_moldiff_MPF = ", 1264 & call_mass_fixer_moldiff_MPF 1265 1257 1266 ! Test of incompatibility: 1258 1267 ! if photochem is used, then water should be used too
Note: See TracChangeset
for help on using the changeset viewer.
