Index: LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcmtest.F
===================================================================
--- LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcmtest.F	(revision 341)
+++ 	(revision )
@@ -1,656 +1,0 @@
-
-      PROGRAM gcm
-
-      USE IOIPSL
-
-      IMPLICIT NONE
-
-c      ......   Version  du 10/01/98    ..........
-
-c             avec  coordonnees  verticales hybrides 
-c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
-
-c=======================================================================
-c
-c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
-c   -------
-c
-c   Objet:
-c   ------
-c
-c   GCM LMD nouvelle grille
-c
-c=======================================================================
-
-c   ...  Dans inigeom , nouveaux calculs pour les elongations  cu , cv 
-c        et possibilite d'appeler une fonction f(y)  a derivee tangente 
-c        hyperbolique a la  place de la fonction a derivee sinusoidale.         
-
-c   ...  Possibilite de choisir le shema de Van-leer pour l'advection de
-c         q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
-c
-c-----------------------------------------------------------------------
-c   Declarations:
-c   -------------
-
-#include "dimensions.h"
-#include "paramet.h"
-#include "comconst.h"
-#include "comdissnew.h"
-#include "comvert.h"
-#include "comgeom.h"
-#include "logic.h"
-#include "temps.h"
-#include "control.h"
-#include "ener.h"
-#include "netcdf.inc"
-#include "description.h"
-#include "serre.h"
-#include "tracstoke.h"
-
-      INTEGER nid
-      INTEGER idim_rlonu, idim_rlonv, idim_rlatu, idim_latv
-      INTEGER idim_s, idim_sig, idim_tim
-
-      INTEGER         longcles
-      PARAMETER     ( longcles = 20 )
-      REAL  clesphy0( longcles )
-      SAVE  clesphy0
-
-      INTEGER*4  iday ! jour julien
-      REAL       time ! Heure de la journee en fraction d'1 jour
-      REAL zdtvr
-      INTEGER nbetatmoy, nbetatdem,nbetat
-
-c   variables dynamiques
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-      REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
-      REAL ps(ip1jmp1)                       ! pression  au sol
-      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
-      REAL pks(ip1jmp1)                      ! exner au  sol
-      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
-      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
-      REAL masse(ip1jmp1,llm)                ! masse d'air
-      REAL phis(ip1jmp1)                     ! geopotentiel au sol
-      REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL w(ip1jmp1,llm)                    ! vitesse verticale
-
-c variables dynamiques intermediaire pour le transport
-      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
-
-c   variables dynamiques au pas -1
-      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
-      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
-      REAL massem1(ip1jmp1,llm)
-
-c   tendances dynamiques
-      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
-      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
-
-c   tendances de la dissipation
-      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
-      REAL dhdis(ip1jmp1,llm)
-
-c   tendances physiques
-      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
-      REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
-
-c   variables pour le fichier histoire
-      REAL dtav      ! intervalle de temps elementaire
-
-      REAL tppn(iim),tpps(iim),tpn,tps
-c
-      INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
-
-      INTEGER itau,itaufinp1,iav
-
-      EXTERNAL caldyn, traceur
-      EXTERNAL dissip,geopot,iniconst,inifilr
-      EXTERNAL integrd,SCOPY
-      EXTERNAL inigeom
-      EXTERNAL exner_hyb,addit
-      EXTERNAL defrun_new, test_period
-      REAL  SSUM
-      REAL time_0 , finvmaold(ip1jmp1,llm)
-
-      LOGICAL lafin
-      INTEGER ij,iq,j,i,l,numvanle,iapp_tracvl
-
-
-      integer histid, histvid, histaveid
-      real time_step, t_wrt, t_ops
-      integer itime_step
-      character*80 dynhist_file, dynhistave_file
-
-
-      REAL rdayvrai,rdaym_ini,rday_ecri
-      LOGICAL first
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
-CMAF pour rajout injection
-      real lon_stat,lat_stat
-      integer i_stat,j_stat
-      LOGICAL injecte
-      DATA injecte/.true./
-      integer ijstat
-CMAF
-cAA
-      character*10 file
-cAA
-
-c-----------------------------------------------------------------------
-c   Initialisations:
-c   ----------------
-
-      descript = 'Run GCM LMDZ'
-      lafin    = .FALSE.
-
-c-----------------------------------------------------------------------
-c
-c  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
-c  ...................................................................
-c
-c     iadv = 1    shema  transport type "humidite specifique LMD"  
-c     iadv = 2    shema   amont
-c     iadv = 3    shema  Van-leer
-c 
-c     dans le tableau q(ij,l,iq) , iq = 1  pour l'eau vapeur
-c                                , iq = 2  pour l'eau liquide
-c         et eventuellement      , iq = 3, nqmx pour les autres traceurs
-c
-c     iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
-c
-      DO iq = 1, nqmx
-       iadv( iq ) = 3
-      ENDDO
-c
-      DO  iq = 1, nqmx
-       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
-     * ,' pour le traceur no ', iq
-       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'
-     * ,' traceur no ', iq
-       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
-     * ,'le traceur no ', iq
-       IF( iadv(iq).LE.0.OR.iadv(iq).GT.3 )   THEN
-       PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
-     * . '
-         STOP 
-       ENDIF
-      ENDDO
-c
-      first = .TRUE.
-      numvanle = nqmx + 1
-      DO  iq = 1, nqmx
-        IF( iadv(iq).EQ.3.AND.first ) THEN
-          numvanle = iq
-          first    = .FALSE.
-        ENDIF
-      ENDDO
-c
-      DO  iq = 1, nqmx
-        IF( iadv(iq).NE.3.AND.iq.GT.numvanle )   THEN
-          PRINT *,' Il y a discontinuite dans le choix du shema de ',
-     *    'Van-leer pour les traceurs . Corriger et repasser . '
-           STOP
-        ENDIF
-        IF( iadv(iq).LT.1.OR.iadv(iq).GT.3 )    THEN
-          PRINT *,' Le choix de  iadv  est errone pour le traceur ',
-     *    iq
-         STOP
-        ENDIF
-      ENDDO
-c
-
-      CALL dynetat0("start.nc",nqmx,vcov,ucov,
-     .              teta,q,masse,ps,phis, time_0)
-
-      OPEN( 99,file ='run.def',status='old',form='formatted')
-      CALL defrun_new( 99, .TRUE. , clesphy0 )
-      CLOSE(99)
-
-c  on recalcule eventuellement le pas de temps
-
-      IF(MOD(day_step,iperiod).NE.0)
-     * STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
-
-      IF(MOD(day_step,iphysiq).NE.0)
-     * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
-
-      zdtvr    = daysec/FLOAT(day_step)
-        IF(dtvr.NE.zdtvr) THEN
-         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
-        ENDIF
-
-c  nombre d'etats dans les fichiers demarrage et histoire
-      nbetatdem = nday / iecri
-      nbetatmoy = nday / periodav + 1
-
-      dtvr = zdtvr
-      CALL iniconst
-      CALL inigeom
-
-      CALL inifilr
-c
-c   ......    P.Le Van    ( modif  le 29/04/97 )   .........  
-c
-      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
-     *                tetagdiv, tetagrot , tetatemp              )
-c
-
-      CALL pression ( ip1jmp1, ap, bp, ps, p       )
-      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-c
-
-c  numero de stockage pour les fichiers de redemarrage:
-
-c-----------------------------------------------------------------------
-c   temps de depart et de fin:
-c   --------------------------
-
-      itau = 0
-      iday = day_ini+itau/day_step
-      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
-         IF(time.GT.1.) THEN
-          time = time-1.
-          iday = iday+1
-         ENDIF
-      itaufin   = nday*day_step
-      itaufinp1 = itaufin +1
-
-      day_end = day_ini + nday
-      PRINT 300, itau,itaufin,day_ini,day_end
-
-      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
-      CALL dynhist0("histoire.nc",day_ini,anne_ini,phis,nqmx)
-
-      ecripar = .TRUE.
-
-      dtav = iperiod*dtvr/daysec
-      dynhistave_file = 'histmoy.nc'
-      t_ops = iperiod*dtvr
-      t_wrt = periodav * daysec
-      CALL initdynav(dynhistave_file,day_ini,anne_ini,dtvr,
-     .              t_ops, t_wrt, nqmx, histaveid)
-
-
-c   Quelques initialisations pour les traceurs
-      call initial0(ijp1llm*nqmx,dq)
-      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
-      istphy=istdyn/iphysiq
-
-c-----------------------------------------------------------------------
-c   Debut de l'integration temporelle:
-c   ----------------------------------
-
-   1  CONTINUE
-c
-      IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
-         CALL  test_period ( ucov,vcov,teta,q,p,phis )
-         PRINT *,' ----   Test_period   OK ! -----', itau
-      ENDIF 
-c
-      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
-      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
-      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
-      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
-      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
-
-      forward = .TRUE.
-      leapf   = .FALSE.
-      dt      =  dtvr
-
-c   ...    P.Le Van .26/04/94  ....
-
-      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
-      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
-
-
-   2  CONTINUE
-
-c-----------------------------------------------------------------------
-
-c   date:
-c   -----
-
-
-c   gestion des appels de la physique et des dissipations:
-c   ------------------------------------------------------
-c
-c   ...    P.Le Van  ( 6/02/95 )  ....
-
-      apphys = .FALSE.
-      statcl = .FALSE.
-      conser = .FALSE.
-      apdiss = .FALSE.
-
-      IF( purmats ) THEN
-         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
-         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
-         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward 
-     $                              .AND.   physic    ) apphys = .TRUE.
-      ELSE
-         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
-         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
-         IF( MOD(itau+1,iphysiq).EQ.0. AND. physic    ) apphys = .TRUE.
-      END IF
-
-c-----------------------------------------------------------------------
-c   calcul des tendances dynamiques:
-c   --------------------------------
-
-      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
-c
-c
-      CALL caldyn 
-     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
-     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
-
-c-----------------------------------------------------------------------
-c   calcul des tendances advection des traceurs (dont l'humidite)
-c   -------------------------------------------------------------
-
-      IF( forward. OR . leapf )  THEN
-c
-        if (injecte) then
-            injecte=.false.
-c On injecte tous les traceurs de la meme facon
-C injection au debut de la simulation
-C !! MAF dans code pi en radians et lat lon en degrees => conversion
-C injection au point -2, 48 en lat lon
-            lon_stat= -2.
-            lat_stat= 48.
-C
-            pi=2.*asin(1.)
-            lon_stat=lon_stat*pi/180.
-            lat_stat=lat_stat*pi/180.
-            CALL coordij(lon_stat,lat_stat,i_stat,j_stat)
-
-            print*,'INJECTION AU POINT ',lon_stat,lat_stat
-     s      ,'     I,J=',i_stat,j_stat,jjp1-j_stat+1
-            ijstat=(j_stat-1)*iip1+i_stat
-C cas avec rnpb on commence a injecter pour le traceur 5
-            do iq=5,nqmx
-               q(ijstat,1,iq)=
-     s         +1.e12/masse(ijstat,1)
-            enddo
-        endif
-
-        DO 15  iq = 1, nqmx
-c
-         IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
-           CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
-
-         ELSE IF( iq.EQ. nqmx )   THEN
-c
-            iapp_tracvl = 5
-c
-cccc     iapp_tracvl est la frequence en pas du groupement des flux
-cccc      de masse pour  Van-Leer dans la routine  tracvl  .
-c
-            CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
-     *                             p,  masse, dq                  )
-c
-         ENDIF
-c
- 15     CONTINUE
-C
-Cmaf stokage du flux de masse pour traceurs OFF-LINE
-          CALL fluxstoke(pbaru,pbarv,masse,teta,phi,phis)
-c
-      ENDIF
-
-
-c-----------------------------------------------------------------------
-c   integrations dynamique et traceurs:
-c   ----------------------------------
-
-       CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
-     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
-     $              finvmaold                                    )
-
-c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
-c
-c-----------------------------------------------------------------------
-c   calcul des tendances physiques:
-c   -------------------------------
-c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
-c
-       IF( purmats )  THEN
-          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
-       ELSE
-          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
-       ENDIF
-c
-c
-       IF( apphys )  THEN
-c
-c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
-c
-
-         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
-         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
-
-           rdaym_ini  = itau * dtvr / daysec
-           rdayvrai   = rdaym_ini  + day_ini
-
-           IF ( ecritphy.LT.1. )  THEN
-             rday_ecri = rdaym_ini
-           ELSE
-             rday_ecri = NINT( rdayvrai )
-           ENDIF
-c
-C MAF le dernier traceur ne passe pas dans la physique 
-C nqmx => nqmx -1 dans calfis et addfi
-C
-        CALL calfis( nqmx-1, lafin ,rdayvrai,rday_ecri,time  ,
-     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
-     $     du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi  )
-
-c      ajout des tendances physiques:
-c      ------------------------------
-          CALL addfi( nqmx-1, dtphys, leapf, forward   ,
-     $                  ucov, vcov, teta , q   ,ps ,
-     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
-c
-       ENDIF
-
-        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-
-c-----------------------------------------------------------------------
-c
-c
-c   dissipation horizontale et verticale  des petites echelles:
-c   ----------------------------------------------------------
-
-      IF(apdiss) THEN
-
-         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
-
-         CALL addit( ijp1llm,ucov ,dudis,ucov )
-         CALL addit( ijmllm ,vcov ,dvdis,vcov )
-         CALL addit( ijp1llm,teta ,dhdis,teta )
-
-
-c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
-c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
-c
-      IF( alphax.NE. 0. ) THEN
-
-        DO l  =  1, llm
-          DO ij =  1,iim
-           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
-           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
-          ENDDO
-           tpn  = SSUM(iim,tppn,1)/apoln
-           tps  = SSUM(iim,tpps,1)/apols
-
-          DO ij = 1, iip1
-           teta(  ij    ,l) = tpn
-           teta(ij+ip1jm,l) = tps
-          ENDDO
-        ENDDO
-
-      ENDIF
-
-      END IF
-        
-c   ********************************************************************
-c   ********************************************************************
-c   .... fin de l'integration dynamique  et physique pour le pas itau ..
-c   ********************************************************************
-c   ********************************************************************
-
-c   preparation du pas d'integration suivant  ......
-
-      IF ( .NOT.purmats ) THEN
-c       ........................................................
-c       ..............  schema matsuno + leapfrog  ..............
-c       ........................................................
-
-            IF(forward. OR. leapf) THEN
-              itau= itau + 1
-              iday= day_ini+itau/day_step
-              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
-                IF(time.GT.1.) THEN
-                  time = time-1.
-                  iday = iday+1
-                ENDIF
-            ENDIF
-
-
-C            IF( itau. EQ. itaufinp1 )    STOP
-            IF( itau. EQ. itaufinp1 )    goto 200
-
-
-c-----------------------------------------------------------------------
-c   ecriture du fichier histoire moyenne:
-c   -------------------------------------
-
-            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
-               IF(itau.EQ.itaufin) THEN
-                  iav=1
-               ELSE
-                  iav=0
-               ENDIF
-               CALL writedynav(histaveid, nqmx, itau,vcov ,
-     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
-
-            ENDIF
-
-c-----------------------------------------------------------------------
-c   ecriture de la bande histoire:
-c   ------------------------------
-
-            IF( MOD(itau,iecri*day_step).EQ.0) THEN
-
-               nbetat = nbetatdem
-       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
-       CALL dynhist1("histoire.nc",time+iday-day_ini,
-     .                     vcov,ucov,teta,q,nqmx,masse,ps)
-
-
-            ENDIF
-
-            IF(itau.EQ.itaufin) THEN
-
-       CALL dynredem1("restart.nc",0.0,
-     .                     vcov,ucov,teta,q,nqmx,masse,ps)
-
-              CLOSE(99)
-            ENDIF
-
-c-----------------------------------------------------------------------
-c   gestion de l'integration temporelle:
-c   ------------------------------------
-
-            IF( MOD(itau,iperiod).EQ.0 )    THEN
-                    GO TO 1
-            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
-
-                   IF( forward )  THEN
-c      fin du pas forward et debut du pas backward
-
-                      forward = .FALSE.
-                        leapf = .FALSE.
-                           GO TO 2
-
-                   ELSE
-c      fin du pas backward et debut du premier pas leapfrog
-
-                        leapf =  .TRUE.
-                        dt  =  2.*dtvr
-                        GO TO 2
-                   END IF
-            ELSE
-
-c      ......   pas leapfrog  .....
-
-                 leapf = .TRUE.
-                 dt  = 2.*dtvr
-                 GO TO 2
-            END IF
-
-      ELSE
-
-c       ........................................................
-c       ..............       schema  matsuno        ...............
-c       ........................................................
-            IF( forward )  THEN
-
-             itau =  itau + 1
-             iday = day_ini+itau/day_step
-             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
-
-                  IF(time.GT.1.) THEN
-                   time = time-1.
-                   iday = iday+1
-                  ENDIF
-
-               forward =  .FALSE.
-C            IF( itau. EQ. itaufinp1 )    STOP
-            IF( itau. EQ. itaufinp1 )    goto 200
-
-               GO TO 2
-
-            ELSE
-
-            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
-               IF(itau.EQ.itaufin) THEN
-                  iav=1
-               ELSE
-                  iav=0
-               ENDIF
-               CALL writedynav(histaveid, nqmx, itau,vcov ,
-     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
-
-            ENDIF
-
-               IF(MOD(itau,iecri*day_step).EQ.0) THEN
-                  nbetat = nbetatdem
-       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
-      CALL dynhist1("histoire.nc",time+iday-day_ini,
-     .                     vcov,ucov,teta,q,nqmx,masse,ps)
-
-               ENDIF
-
-                 IF(itau.EQ.itaufin)
-     . CALL dynredem1("restart.nc",0.0,
-     .                     vcov,ucov,teta,q,nqmx,masse,ps)
-
-                 forward = .TRUE.
-                 GO TO  1
-
-            ENDIF
-
-      END IF
-
- 200  continue
-      call histclo
-
-
- 300  FORMAT('1'/15x'run du pas'i7,2x,'au pas'i7,2x,  'c"est a dire du
-     * jour'i7,3x'au jour'i7//)
-
-      STOP
-      END
