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

Timestamp:

Sep 23, 2013, 9:56:47 AM (11 years ago)

Author:

emillour

Message:

Mars GCM:

• IMPORTANT CHANGE: Removed all reference/use of ngridmx (dimphys.h) in routines (necessary prerequisite to using parallel dynamics); in most cases this just means adding 'ngrid' as routine argument, and making local saved variables allocatable (and allocated at first call). In the process, had to convert many *.h files to equivalent modules: yomaer.h => yomaer_h.F90 , surfdat.h => surfdat_h.F90 , comsaison.h => comsaison_h.F90 , yomlw.h => yomlw_h.F90 , comdiurn.h => comdiurn_h.F90 , dimradmars.h => dimradmars_mod.F90 , comgeomfi.h => comgeomfi_h.F90, comsoil.h => comsoil_h.F90 , slope.h => slope_mod.F90
• Also updated EOF routines, everything is now in eofdump_mod.F90
• Removed unused routine lectfux.F (in dyn3d)

EM

File:

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

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

@@ -1 @@
-subroutine moldiff_red(pplay,pplev,pt,pdt,pq,pdq,ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
-
-use tracer_mod, only: nqmx, noms, mmol
+subroutine moldiff_red(ngrid,nlayer,nq,pplay,pplev,pt,pdt,pq,pdq,&
+                       ptimestep,zzlay,pdteuv,pdtconduc,pdqdiff)
+
+use tracer_mod, only: noms, mmol
 
 implicit none
 
-#include "dimensions.h"
-#include "dimphys.h"
+!#include "dimensions.h"
+!#include "dimphys.h"
 #include "comcstfi.h"
-#include "callkeys.h"
-#include "comdiurn.h"
-#include "chimiedata.h"
+!#include "callkeys.h"
+!#include "comdiurn.h"
+!#include "chimiedata.h"
 !#include "tracer.h"
-#include "conc.h"
+!#include "conc.h"
 #include "diffusion.h"
 
@@ -19 +20 @@
 ! Input/Output
 !
+      integer,intent(in) :: ngrid ! number of atmospheric columns
+      integer,intent(in) :: nlayer ! number of atmospheric layers
+      integer,intent(in) :: nq ! number of advected tracers
       real ptimestep
-      real pplay(ngridmx,nlayermx)
-      real zzlay(ngridmx,nlayermx)
-      real pplev(ngridmx,nlayermx+1)
-      real pq(ngridmx,nlayermx,nqmx)
-      real pdq(ngridmx,nlayermx,nqmx)
-      real pt(ngridmx,nlayermx)
-      real pdt(ngridmx,nlayermx)
-      real pdteuv(ngridmx,nlayermx)
-      real pdtconduc(ngridmx,nlayermx)
-      real pdqdiff(ngridmx,nlayermx,nqmx)
+      real pplay(ngrid,nlayer)
+      real zzlay(ngrid,nlayer)
+      real pplev(ngrid,nlayer+1)
+      real pq(ngrid,nlayer,nq)
+      real pdq(ngrid,nlayer,nq)
+      real pt(ngrid,nlayer)
+      real pdt(ngrid,nlayer)
+      real pdteuv(ngrid,nlayer)
+      real pdtconduc(ngrid,nlayer)
+      real pdqdiff(ngrid,nlayer,nq)
 !
 ! Local
@@ -37 +41 @@
 !      real hco2(ncompdiff),ho
 
-      integer,dimension(nqmx) :: indic_diff
+      integer,dimension(nq) :: indic_diff
       integer ig,iq,nz,l,k,n,nn,p,ij0
       integer istep,il,gcn,ntime,nlraf
@@ -44 +48 @@
       real*8 rho0,D0,T0,H0,time0,dZ,time,dZraf,tdiff,Zmin,Zmax
       real*8 FacEsc,invsgmu
-      real*8 hp(nlayermx)
-      real*8 pp(nlayermx)
-      real*8 pint(nlayermx)
-      real*8 tt(nlayermx),tnew(nlayermx),tint(nlayermx)
-      real*8 zz(nlayermx)
+      real*8 hp(nlayer)
+      real*8 pp(nlayer)
+      real*8 pint(nlayer)
+      real*8 tt(nlayer),tnew(nlayer),tint(nlayer)
+      real*8 zz(nlayer)
       real*8,dimension(:,:),allocatable,save :: qq,qnew,qint,FacMass
       real*8,dimension(:,:),allocatable,save :: rhoK,rhokinit
-      real*8 rhoT(nlayermx)
-      real*8 dmmeandz(nlayermx)
-      real*8 massemoy(nlayermx)
+      real*8 rhoT(nlayer)
+      real*8 dmmeandz(nlayer)
+      real*8 massemoy(nlayer)
       real*8,dimension(:),allocatable :: Praf,Traf,Rraf,Mraf,Nraf,Draf,Hraf,Wraf
       real*8,dimension(:),allocatable :: Zraf,Tdiffraf
@@ -130 +134 @@
 
         ncompdiff=0
-        indic_diff(1:nqmx)=0
+        indic_diff(1:nq)=0
         
-        do nn=1,nqmx
+        do nn=1,nq
         do n=1,14
         if (ListeDiff(n) .eq. noms(nn)) then
@@ -151 +155 @@
 ! Store gcm indexes in gcmind
         n=0
-        do nn=1,nqmx
+        do nn=1,nq
         if (indic_diff(nn) .eq. 1) then
         n=n+1
@@ -162 +166 @@
 ! find vertical index above which diffusion is computed
 
-        do l=1,nlayermx
+        do l=1,nlayer
         if (pplay(1,l) .gt. Pdiff) then
         il0=l
@@ -176 +180 @@
 
 ! allocatation des tableaux dependants du nombre d especes diffusees
-        allocate(qq(nlayermx,ncompdiff))
-        allocate(qnew(nlayermx,ncompdiff))
-        allocate(qint(nlayermx,ncompdiff))
-        allocate(FacMass(nlayermx,ncompdiff))
-        allocate(rhok(nlayermx,ncompdiff))
-        allocate(rhokinit(nlayermx,ncompdiff))
+        allocate(qq(nlayer,ncompdiff))
+        allocate(qnew(nlayer,ncompdiff))
+        allocate(qint(nlayer,ncompdiff))
+        allocate(FacMass(nlayer,ncompdiff))
+        allocate(rhok(nlayer,ncompdiff))
+        allocate(rhokinit(nlayer,ncompdiff))
 
         allocate(wi(ncompdiff))
@@ -212 +216 @@
         
 !       print*,'moldiff',i_h2,i_h,ncompdiff
-      do ig=1,ngridmx
+      do ig=1,ngrid
         pp=dble(pplay(ig,:))
 
@@ -218 +222 @@
 
 !       CALL TMNEW(pt(ig,:),pdt(ig,:),pdtconduc(ig,:),pdteuv(ig,:)      & 
-!     &  ,tt,ptimestep,nlayermx,ig)
-        do l=1,nlayermx
+!     &  ,tt,ptimestep,nlayer,ig)
+        do l=1,nlayer
           tt(l)=pt(ig,l)*1D0+(pdt(ig,l)*dble(ptimestep)+                &
                               pdtconduc(ig,l)*dble(ptimestep)+          &
@@ -228 +232 @@
               pdt(ig,l),pdtconduc(ig,l),pdteuv(ig,l),dble(ptimestep)
           endif
-        enddo ! of do l=1,nlayermx
+        enddo ! of do l=1,nlayer
 
 ! Update the mass mixing ratios modified by other processes
 
-!       CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayermx,        &
+!       CALL QMNEW(pq(ig,:,:),pdq(ig,:,:),qq,ptimestep,nlayer,        &
 !     &  ncompdiff,gcmind,ig)
         do iq=1,ncompdiff
-         do l=1,nlayermx
+         do l=1,nlayer
           qq(l,iq)=pq(ig,l,gcmind(iq))*1D0+(                            &
                            pdq(ig,l,gcmind(iq))*dble(ptimestep))
           qq(l,iq)=max(qq(l,iq),1d-30)
-         enddo ! of do l=1,nlayermx
+         enddo ! of do l=1,nlayer
         enddo ! of do iq=1,ncompdiff
 
 ! Compute the Pressure scale height
 
-        CALL HSCALE(pp,hp,nlayermx)
+        CALL HSCALE(pp,hp,nlayer)
 
 ! Compute the atmospheric mass (in Dalton)
 
-        CALL MMOY(massemoy,mmol,qq,gcmind,nlayermx,ncompdiff)
+        CALL MMOY(massemoy,mmol,qq,gcmind,nlayer,ncompdiff)
 
 ! Compute the vertical gradient of atmospheric mass
 
-        CALL DMMOY(massemoy,hp,dmmeandz,nlayermx)
+        CALL DMMOY(massemoy,hp,dmmeandz,nlayer)
 
 ! Compute the altitude of each layer
 
-        CALL ZVERT(pp,tt,massemoy,zz,nlayermx,ig)
+        CALL ZVERT(pp,tt,massemoy,zz,nlayer,ig)
 
 ! Compute the total mass density (kg/m3)
         
-        CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayermx,ncompdiff)
+        CALL RHOTOT(pp,tt,massemoy,qq,RHOT,RHOK,nlayer,ncompdiff)
         RHOKINIT=RHOK
 
@@ -271 +275 @@
         Mtot1(1:ncompdiff)=0d0
 
-        do l=il0,nlayermx
+        do l=il0,nlayer
         do nn=1,ncompdiff
         Mtot1(nn)=Mtot1(nn)+1d0/g*qq(l,nn)*                             &
@@ -279 +283 @@
 
         Zmin=zz(il0)
-        Zmax=zz(nlayermx)
+        Zmax=zz(nlayer)
 
 
@@ -286 +290 @@
         if (Zmax .gt. 4000000.) then
         Print*,'Zmax too high',ig,zmax,zmin
-        do l=1,nlayermx
+        do l=1,nlayer
         print*,'old',zz(l),pt(ig,l),pdteuv(ig,l),pdq(ig,l,:)
         print*,'l',l,rhot(l),tt(l),pp(l),massemoy(l),qq(l,:)
@@ -321 +325 @@
         CALL UPPER_RESOL(pp,tt,zz,massemoy,RHOT,RHOK,                   &
      &  qq,mmol,gcmind,Praf,Traf,Qraf,Mraf,Zraf,                        &
-     &  Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayermx,ig)
+     &  Nraf,Nrafk,Rraf,Rrafk,il0,nlraf,ncompdiff,nlayer,ig)
 
         Prafold=Praf
@@ -330 +334 @@
 
         CALL GCMGRID_P(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qint,tt,tint  &
-     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,ig)
+     &    ,pp,mmol,gcmind,nlraf,ncompdiff,nlayer,ig)
 
 ! We compute the mass correction factor of each specie at each pressure level
 
-        CALL CORRMASS(qq,qint,FacMass,nlayermx,ncompdiff)
+        CALL CORRMASS(qq,qint,FacMass,nlayer,ncompdiff)
 
 ! Altitude step
@@ -367 +371 @@
 !       enddo
 
-!       do l=1,nlayermx
+!       do l=1,nlayer
 !       print*,'l',l,zz(l),pp(l),tt(l),sum(qq(l,:)),massemoy(l)
 !       enddo
@@ -374 +378 @@
 ! No change below il0
         
-        do l=1,nlayermx
+        do l=1,nlayer
          qnew(l,:)=qq(l,:) ! No effet below il0
        enddo
@@ -527 +531 @@
         print*,'Mraf',Mraf
         stop
-!       pdqdiff(1:ngridmx,1:nlayermx,1:nqmx)=0.
+!       pdqdiff(1:ngrid,1:nlayer,1:nq)=0.
 !       return
 !       Rrafk(l,nn)=1D-30*Rraf(l)
@@ -572 +576 @@
 
         CALL GCMGRID_P2(Zraf,Praf,Qraf,Traf,Nrafk,Rrafk,qq,qnew,tt,tnew,&
-     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayermx,FacMass,ig)
-
-        CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayermx,ncompdiff)
+     &   pp,mmol,gcmind,nlraf,ncompdiff,nlayer,FacMass,ig)
+
+        CALL RHOTOT(pp,tt,massemoy,qnew,RHOT,RHOK,nlayer,ncompdiff)
 
         if (ig .eq. ij0) then 
-        do l=il0,nlayermx
+        do l=il0,nlayer
         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,:)),&
@@ -590 +594 @@
         Mtot2(1:ncompdiff)=0d0
 
-        do l=il0,nlayermx
+        do l=il0,nlayer
         do nn=1,ncompdiff
         Mtot2(nn)=Mtot2(nn)+1d0/g*qnew(l,nn)*                           &
@@ -600 +604 @@
 
 !       do nn=1,ncompdiff
-!       CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayermx,nn,ncompdiff)
+!       CALL CheckMass2(qq,qnew,pplev(ig,:),il0,nlayer,nn,ncompdiff)
 !       enddo
 
 ! Compute the diffusion trends du to diffusion
 
-        do l=1,nlayermx
+        do l=1,nlayer
         do nn=1,ncompdiff
         pdqdiff(ig,l,gcmind(nn))=(qnew(l,nn)-qq(l,nn))/ptimestep