Changeset 2467 for trunk/LMDZ.MARS/libf/aeronomars

Timestamp:

Mar 3, 2021, 1:48:21 PM (4 years ago)

Author:

emillour

Message:

Mars GCM:

• bug fix in the diagnostic of H, H2, and D total escape which was wrong when running in parallel. The escape rates can now also be included as outputs in diagfi.

FGG

Location:

trunk/LMDZ.MARS/libf/aeronomars

Files:

• moldiff_red.F90 (modified) (6 diffs)
• thermosphere.F (modified) (3 diffs)

trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90

@@ -1 +1 @@
 subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
-                       ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
+                       ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff,&
+                       PhiEscH,PhiEscH2,PhiEscD)
 
 use tracer_mod, only: noms, mmol
 use geometry_mod, only: cell_area
+use planetwide_mod, only: planetwide_sumval
 
 implicit none
@@ -30 +32 @@
       real pdtconduc(ngrid,nlayer)
       real pdqdiff(ngrid,nlayer,nq)
+      real*8 PhiEscH,PhiEscH2,PhiEscD
 !
 ! Local
@@ -43 +46 @@
       real*8 masseU,kBolt,RgazP,Rmars,Grav,Mmars
       real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
-      real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2,PhiEscD
+      real*8 FacEsc,invsgmu 
+      real*8 PhiauxH(ngrid),PhiauxH2(ngrid),PhiauxD(ngrid)
       real*8 hp(nlayer)
       real*8 pp(nlayer)
@@ -598 +602 @@
 ! the trend only at the end
 
-        if (i_h .ne. 1000) PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1
-        if (i_h2 .ne. 1000) PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3)
-        if (i_d .ne. 1000) PhiEscD=PhiEscD+wi(i_d)*Nrafk(nlraf,i_d)*cell_area(ig)
+!       if (i_h .ne. 1000) PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1
+        if (i_h .ne. 1000) PhiauxH(ig)=wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1
+!       if (i_h2 .ne. 1000) PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig) ! in s-1 (U in m/s, aire in m2, Nrafk in m-3)
+        if (i_h2 .ne. 1000) PhiauxH2(ig)=wi(i_h2)*Nrafk(nlraf,i_h2)*cell_area(ig)
+!        if (i_d .ne. 1000) PhiEscD=PhiEscD+wi(i_d)*Nrafk(nlraf,i_d)*cell_area(ig)
+        if (i_d .ne. 1000) PhiauxD(ig)=wi(i_d)*Nrafk(nlraf,i_d)*cell_area(ig)
 !       print*,'test',ig,wi(i_h),wi(i_h2),Nrafk(nlraf,i_h),Nrafk(nlraf,i_h2),Nrafk(nlraf,i_d),cell_area(ig),PhiEscH,PhiEscH2,i_h,i_h2,i_d,PhiEscD
 !       stop
@@ -657 +664 @@
        enddo  ! ig loop         
 !       print*,'Escape flux H, H2,D (s-1)',PhiEscH,PhiEscH2,PhiEscD
+       if (i_h.ne.1000) call planetwide_sumval(PhiauxH,PhiEscH)
+       if (i_h2.ne.1000) call planetwide_sumval(PhiauxH2,PhiEscH2)
+       if (i_d.ne.1000) call planetwide_sumval(PhiauxD,PhiEscD)
 
       return
@@ -1178 +1188 @@
         endif
         enddo
-!       D(l)=1d0/interm
+        !Temporary: eliminate modification to include Wilke's formulation
+        !back to the old scheme to check effect 
+        !D(l)=1d0/interm
         D(l)=(1D0-Nk(l,nn)/N(l))/interm
         enddo

trunk/LMDZ.MARS/libf/aeronomars/thermosphere.F

@@ -3 +3 @@
      $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
      &     pt,pq,pu,pv,pdt,pdq,
-     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
+     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff,
+     $     PhiEscH,PhiEscH2,PhiEscD)
 
       use conc_mod, only: rnew, cpnew
@@ -34 +35 @@
       REAL,INTENT(out) :: zdvmolvis(ngrid,nlayer)
       REAL,INTENT(out) :: zdqmoldiff(ngrid,nlayer,nq)
+      REAL*8,INTENT(out) :: PhiEscH,PhiEscH2,PhiEscD
 
       INTEGER :: l,ig
@@ -80 +82 @@
         call moldiff_red(ngrid,nlayer,nq,
      &                   pplay,pplev,pt,pdt,pq,pdq,ptimestep,
-     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff)
+     &                   zzlay,zdteuv,zdtconduc,zdqmoldiff,
+     &                   PhiEscH,PhiEscH2,PhiEscD)
       endif