Changeset 2461 for trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90

Timestamp:

Feb 16, 2021, 12:31:33 PM (4 years ago)

Author:

emillour

Message:

Mars GCM:

• Adding the deuterium chemistry now that the HDO cycle is included.
• Chemistry still works as before if deuterium tracers are not present.
• Added handling of hdo in molecular diffusion (moldiff_red).

FGG+JYC+EM

File:

• trunk/LMDZ.MARS/libf/aeronomars/moldiff_red.F90 (modified) (10 diffs)

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

@@ -43 +43 @@
       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
+      real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2,PhiEscD
       real*8 hp(nlayer)
       real*8 pp(nlayer)
@@ -64 +64 @@
       real*8,dimension(:),allocatable,save ::  wi,Wad,Uthermal,Lambdaexo,Hspecie
       real*8,dimension(:),allocatable,save :: Mtot1,Mtot2,Mraf1,Mraf2
-      integer,parameter :: ListeDiffNb=15
+      integer,parameter :: ListeDiffNb=16
       character(len=20),dimension(ListeDiffNb) :: ListeDiff
       real*8,parameter :: pi=3.141592653
@@ -75 +75 @@
       integer,save :: i_h2  
       integer,save :: i_h
+      integer,save :: i_d
 ! vertical index limit for the molecular diffusion
       integer,save :: il0  
@@ -130 +131 @@
         ListeDiff(14)='n'
         ListeDiff(15)='he'
+        ListeDiff(16)='hdo_vap'
         i_h=1000
         i_h2=1000
+        i_d=1000
 ! On regarde quelles especes diffusables sont presentes
 
@@ -162 +165 @@
         if (noms(nn) .eq. 'h') i_h=n
         if (noms(nn) .eq. 'h2') i_h2=n
+        if (noms(nn) .eq. 'd') i_d=n
         endif
         enddo
@@ -222 +226 @@
         PhiEscH=0D0
         PhiEscH2=0D0
+        PhiEscD=0D0
 
       do ig=1,ngrid
@@ -403 +408 @@
      &  (Rmars+Zraf(nlraf))/kbolt/Traf(nlraf)
         wi(nn)=0D0
-        if (nn .eq. i_h .or. nn .eq. i_h2) then 
+        if (nn .eq. i_h .or. nn .eq. i_h2 .or. nn .eq. i_d) then 
         wi(nn)=Uthermal(nn)/2d0/sqrt(PI)*exp(-Lambdaexo(nn))*           &
      &  (Lambdaexo(nn)+1d0)
@@ -409 +414 @@
         enddo
 
-!       print*,'wi',wi(i_h),wi(i_h2),Uthermal,Lambdaexo,mmol(gcmind(:))
+!       print*,'wi',wi(i_h),wi(i_h2),wi(i_d),Uthermal,Lambdaexo,mmol(gcmind(:))
 !       print*,'wi',wi
-!       stop
 
 ! Compute time step for diffusion
@@ -594 +598 @@
 ! the trend only at the end
 
-        PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*cell_area(ig) ! in s-1
-        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)
-!       print*,'test',ig,wi(i_h),Nrafk(nlraf,i_h),wi(i_h2),Nrafk(nlraf,i_h2),cell_area(ig),PhiEscH,PhiEscH2,i_h,i_h2
+        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)
+!       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
 
@@ -651 +656 @@
 
        enddo  ! ig loop         
-!       print*,'Escape flux H, H2 (s-1)',PhiEscH,PhiEscH2
+!       print*,'Escape flux H, H2,D (s-1)',PhiEscH,PhiEscH2,PhiEscD
 
       return