Changeset 710 for trunk/LMDZ.MARS/libf

Timestamp:

Jun 28, 2012, 10:02:21 AM (13 years ago)

Author:

emillour

Message:

Mars GCM:

• Decreased time step for molecular diffusion (diffusion.h) for better stability
• Added test in moldiff_red.F90 to enforce that the system can't be stuck in an infinite loop.
• The creation and output of coefficients in file 'coeffs.dat' in moldiff_red.F is now controled by an internal flag (by default set to false).

JYC + EM

Location:

trunk/LMDZ.MARS/libf/aeronomars

Files:

• diffusion.h (modified) (1 diff)
• moldiff_red.F90 (modified) (10 diffs)
• moldiffcoeff_red.F (modified) (3 diffs)

trunk/LMDZ.MARS/libf/aeronomars/diffusion.h

@@ -10 +10 @@
 
       parameter (Pdiff=15.)      ! pressure below which diffusion is computed
-      parameter (tdiffmin=10d0)
+      parameter (tdiffmin=5d0)
       parameter (dzres=2d0)      ! grid resolution (km) for diffusion
       

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

@@ -97 +97 @@
       integer,save :: ncompdiff
       integer,dimension(:),allocatable,save :: gcmind ! array of GCM indexes
-      integer ierr
+      integer ierr,compteur
 
       logical,save :: firstcall=.true.
@@ -312 +312 @@
 
         CALL GCMGRID_P(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qint,tt,tint  &
-     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayermx)
+     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,ig)
 
 ! We compute the mass correction factor of each specie at each pressure level
@@ -554 +554 @@
 
         CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
-     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,FacMass)
+     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,FacMass,ig)
 
         CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayermx,ncompdiff)
@@ -561 +561 @@
         do l=il0,nlayermx
         write(*,'(i2,1x,19(e12.4,1x))') l,zz(l),tt(l),RHOK(l,1)/sum(RHOK(l,:)),RHOKINIT(l,1)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,2)/sum(RHOK(l,:)),RHOKINIT(l,2)/sum(RHOKINIT(l,:)),&
+     &  RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
      &  RHOK(l,5)/sum(RHOK(l,:)),RHOKINIT(l,5)/sum(RHOKINIT(l,:)),&
-     &  RHOK(l,3)/sum(RHOK(l,:)),RHOKINIT(l,3)/sum(RHOKINIT(l,:)),&
-     &  RHOK(l,6)/sum(RHOK(l,:)),RHOKINIT(l,6)/sum(RHOKINIT(l,:)),&
      &  RHOK(l,7)/sum(RHOK(l,:)),RHOKINIT(l,7)/sum(RHOKINIT(l,:))
         enddo
@@ -1352 +1352 @@
 
         SUBROUTINE GCMGRID_P(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew,             &
-     &    pp,M,gc,nl,nq,nlx)
+     &    pp,M,gc,nl,nq,nlx,ig)
         IMPLICIT NONE
 #include "dimensions.h"
-        INTEGER :: nl,nq,nlx,il,nn,iP
+        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
         INTEGER,DIMENSION(1) :: indP
         INTEGER,DIMENSION(nq) :: gc 
@@ -1412 +1412 @@
         Pnew2=Pnew
 
+        compteur=0
         do while (pnew2 .ge. pp(il))    
+        compteur=compteur+1
         do nn=1,nq
         Hi=Konst*T(nl)/dble(M(gc(nn)))/g
@@ -1421 +1423 @@
         Znew=Znew2
         Znew2=Znew2+Dz
+        if (compteur .ge. 100000) then
+        print*,'error moldiff_red infinite loop'
+        print*,ig,il,pp(il),tt(nl),pnew2,qnew(il,:),Znew2
+        stop
+        endif
 !       print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
         enddo
@@ -1446 +1453 @@
 
         SUBROUTINE GCMGRID_P2(Z,P,Q,T,Nk,Rk,qq,qnew,tt,tnew             &
-     &   ,pp,M,gc,nl,nq,nlx,facM)
+     &   ,pp,M,gc,nl,nq,nlx,facM,ig)
         IMPLICIT NONE
 #include "dimensions.h"
-        INTEGER :: nl,nq,nlx,il,nn,iP
+        INTEGER :: nl,nq,nlx,il,nn,iP,ig,compteur
         INTEGER,DIMENSION(1) :: indP
         INTEGER,DIMENSION(nq) :: gc
@@ -1505 +1512 @@
         Pnew2=Pnew
 
+        compteur=0
         do while (pnew2 .ge. pp(il))    
+        compteur=compteur+1
         do nn=1,nq
         Hi=Konst*T(nl)/dble(M(gc(nn)))/g
@@ -1514 +1523 @@
         Znew=Znew2
         Znew2=Znew2+Dz
+        if (compteur .ge. 100000) then
+        print*,'pb moldiff_red infinite loop'
+        print*,ig,nl,T(nl),pnew2,qnew(il,:),Znew2
+        stop
+        endif
+
 !       print*,'test',Pnew2,Znew2,Nnew(nq),pp(il)
         enddo

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

@@ -87 +87 @@
       real dnh
       logical,save :: firstcall=.true.
+      logical,parameter :: outputcoeffs=.false. ! to output 'coeffs.dat' file,
+                                                ! set outputcoeffs=.true.
 
 ! Initializations at first call (and some sanity checks)
@@ -253 +255 @@
         if (noms(n) .eq. 'o') g_o=n
         enddo
-        print*,'gh2',g_h2,g_h,g_o
-
-       print*,'moldiffcoef: COEFF CALC'
-       open(56,file='coeffs.dat',status='unknown')
+        print*,'moldiffcoeff_red: gh2',g_h2,g_h,g_o
+
+       print*,'moldiffcoeff_red: COEFF CALC'
+
       do n=1,ncompdiff2
         dijref=0.
@@ -283 +285 @@
       enddo 
 
-      do n=1,ncompdiff2
+      if (outputcoeffs) then
+       ! output coefficients in 'coeffs.dat' file
+       open(56,file='coeffs.dat',status='unknown')
+       do n=1,ncompdiff2
         do nn=n,ncompdiff2
           write(56,*) n,nn,dij(n,nn)    !*1.e5/1.381e-23/(273**1.75)
         enddo
-      enddo
-      close(56)
+       enddo
+       close(56)
+      endif