Changeset 1357

Timestamp:

Oct 20, 2014, 8:18:06 AM (10 years ago)

Author:

emillour

Message:

Mars GCM:

• Update and bug correction on the escape fluxes of H and H2.

JYC

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/aeronomars/moldiff_red.F90 (modified) (10 diffs)

trunk/LMDZ.MARS/README

@@ -2147 +2147 @@
   measurements, along with dustiropacity='mcs'. Future evolutions of this option
   could include other instruments, or be used for water ice.
+
+== 20/10/2014 == JYC
+- Update and bug correction on the escape fluxes of H and H2.

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

@@ -3 +3 @@
 
 use tracer_mod, only: noms, mmol
-USE comcstfi_h
-
+use comgeomphy, only: airephy
 implicit none
 
+!#include "dimensions.h"
+!#include "dimphys.h"
+#include "comcstfi.h"
+!#include "callkeys.h"
+!#include "comdiurn.h"
+!#include "chimiedata.h"
+!#include "tracer.h"
+!#include "conc.h"
 #include "diffusion.h"
+
+! July 2014 JYC ADD BALISTIC Transport coupling to compute wup for H and H2
+
+
 
 !
@@ -39 +50 @@
       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
+      real*8 FacEsc,invsgmu,PhiEscH,PhiEscH2
       real*8 hp(nlayer)
       real*8 pp(nlayer)
@@ -132 +143 @@
         if (ListeDiff(n) .eq. noms(nn)) then
         indic_diff(nn)=1
-        if (noms(nn) .eq. 'h') i_h=n
-        if (noms(nn) .eq. 'h2') i_h2=n
+!       if (noms(nn) .eq. 'h') i_h=n
+!       if (noms(nn) .eq. 'h2') i_h2=n
         endif
 
@@ -151 +162 @@
         n=n+1
         gcmind(n)=nn
-        endif
-        enddo
-
-        print*,'gcmind',gcmind
+        if (noms(nn) .eq. 'h') i_h=n
+        if (noms(nn) .eq. 'h2') i_h2=n
+        endif
+        enddo
+
+!       print*,'gcmind',gcmind,i_h,i_h2
 
 ! find vertical index above which diffusion is computed
@@ -207 +220 @@
         invsgmu=1d0/g/masseU
         
-!       print*,'moldiff',i_h2,i_h,ncompdiff
+! Compute the wup(ig) for H and H2 using the balistic code from R Yelle
+
+        PhiEscH=0D0
+        PhiEscH2=0D0
+
       do ig=1,ngrid
         pp=dble(pplay(ig,:))
@@ -393 +410 @@
         endif
         enddo
+
+!       print*,'wi',wi(i_h),wi(i_h2),Uthermal,Lambdaexo,mmol(gcmind(:))
+!       print*,'wi',wi
+!       stop
 
 ! Compute time step for diffusion
@@ -572 +593 @@
         CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
 
+! Update total escape flux of H and H2 (if q was really qnew, but not forget we will output
+! the trend only at the end
+
+        PhiEscH=PhiEscH+wi(i_h)*Nrafk(nlraf,i_h)*airephy(ig) ! in s-1
+        PhiEscH2=PhiEscH2+wi(i_h2)*Nrafk(nlraf,i_h2)*airephy(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),airephy(ig),PhiEscH,PhiEscH2,i_h,i_h2
+!       stop
+
+
         if (ig .eq. ij0) then 
         do l=il0,nlayer
@@ -623 +653 @@
 
        enddo  ! ig loop         
-
+!       print*,'Escape flux H, H2 (s-1)',PhiEscH,PhiEscH2
 
       return
@@ -1336 +1366 @@
         Ytot=sum(Y)
 
-        if (abs((Xtot-Ytot)/Xtot) .gt. 1d-4) then
+        if (abs((Xtot-Ytot)/Xtot) .gt. 1d-3) then
         print*,'no conservation for mass',Xtot,Ytot,nn
         endif
@@ -1357 +1387 @@
 !       print*,'l',l,qold(l,nn),qnew(l,nn),Mold,Mnew,DM,P(l),P(l+1)
         ENDDO
-        IF (abs(DM/Mold) .gt. 1d-5) THEN
+        IF (abs(DM/Mold) .gt. 1d-2) THEN
         Print*,'We dont conserve mas',nn,DM,Mold,Mnew,DM/Mold
         ENDIF