Ignore:
Timestamp:
Oct 21, 2024, 2:58:45 PM (23 hours ago)
Author:
abarral
Message:

Convert fixed-form to free-form sources .F -> .{f,F}90
(WIP: some .F remain, will be handled in subsequent commits)

File:
1 moved

Legend:

Unmodified
Added
Removed

  • LMDZ6/trunk/libf/dyn3d/leapfrog.F90

    r5245 r5246  
    22! $Id$
    33!
    4 c
    5 c
    6       SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
    7 
    8 
    9 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
    10 #ifdef CPP_IOIPSL
    11       use IOIPSL
    12 #endif
    13       USE infotrac, ONLY: nqtot, isoCheck
    14       USE guide_mod, ONLY : guide_main
    15       USE write_field, ONLY: writefield
    16       USE control_mod, ONLY: nday, day_step, planet_type, offline,
    17      &                       iconser, iphysiq, iperiod, dissip_period,
    18      &                       iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
    19      &                       periodav, ok_dyn_ave, output_grads_dyn
    20       use exner_hyb_m, only: exner_hyb
    21       use exner_milieu_m, only: exner_milieu
    22       USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs
    23       USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
    24       USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys,
    25      &                     statcl,conser,apdiss,purmats,ok_strato
    26       USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref,
    27      &                        start_time,dt
    28       USE strings_mod, ONLY: msg
    29 
    30       IMPLICIT NONE
    31 
    32 c      ......   Version  du 10/01/98    ..........
    33 
    34 c             avec  coordonnees  verticales hybrides
    35 c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
    36 
    37 c=======================================================================
    38 c
    39 c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
    40 c   -------
    41 c
    42 c   Objet:
    43 c   ------
    44 c
    45 c   GCM LMD nouvelle grille
    46 c
    47 c=======================================================================
    48 c
    49 c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
    50 c      et possibilite d'appeler une fonction f(y)  a derivee tangente
    51 c      hyperbolique a la  place de la fonction a derivee sinusoidale.
    52 
    53 c  ... Possibilite de choisir le shema pour l'advection de
    54 c        q  , en modifiant iadv dans traceur.def  (10/02) .
    55 c
    56 c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
    57 c      Pour Van-Leer iadv=10
    58 c
    59 c-----------------------------------------------------------------------
    60 c   Declarations:
    61 c   -------------
    62 
    63       include "dimensions.h"
    64       include "paramet.h"
    65       include "comdissnew.h"
    66       include "comgeom.h"
    67       include "description.h"
    68       include "iniprint.h"
    69       include "academic.h"
    70 
    71       REAL,INTENT(IN) :: time_0 ! not used
    72 
    73 c   dynamical variables:
    74       REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm)    ! zonal covariant wind
    75       REAL,INTENT(INOUT) :: vcov(ip1jm,llm)      ! meridional covariant wind
    76       REAL,INTENT(INOUT) :: teta(ip1jmp1,llm)    ! potential temperature
    77       REAL,INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
    78       REAL,INTENT(INOUT) :: masse(ip1jmp1,llm)   ! air mass
    79       REAL,INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
    80       REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
    81 
    82       REAL p (ip1jmp1,llmp1  )               ! interlayer pressure
    83       REAL pks(ip1jmp1)                      ! exner at the surface
    84       REAL pk(ip1jmp1,llm)                   ! exner at mid-layer
    85       REAL pkf(ip1jmp1,llm)                  ! filtered exner at mid-layer
    86       REAL phi(ip1jmp1,llm)                  ! geopotential
    87       REAL w(ip1jmp1,llm)                    ! vertical velocity
    88 
    89       real zqmin,zqmax
    90 
    91 c variables dynamiques intermediaire pour le transport
    92       REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
    93 
    94 c   variables dynamiques au pas -1
    95       REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
    96       REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
    97       REAL massem1(ip1jmp1,llm)
    98 
    99 c   tendances dynamiques
    100       REAL dv(ip1jm,llm),du(ip1jmp1,llm)
    101       REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
    102 
    103 c   tendances de la dissipation
    104       REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
    105       REAL dtetadis(ip1jmp1,llm)
    106 
    107 c   tendances physiques
    108       REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
    109       REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
    110 
    111 c   variables pour le fichier histoire
    112       REAL dtav      ! intervalle de temps elementaire
    113 
    114       REAL tppn(iim),tpps(iim),tpn,tps
    115 c
    116       INTEGER itau,itaufinp1,iav
    117 !      INTEGER  iday ! jour julien
    118       REAL       time
    119 
    120       REAL  SSUM
    121 !     REAL finvmaold(ip1jmp1,llm)
    122 
    123 cym      LOGICAL  lafin
    124       LOGICAL :: lafin=.false.
    125       INTEGER ij,iq,l
    126       INTEGER ik
    127 
    128       real time_step, t_wrt, t_ops
    129 
    130 !      REAL rdayvrai,rdaym_ini
    131 ! jD_cur: jour julien courant
    132 ! jH_cur: heure julienne courante
    133       REAL :: jD_cur, jH_cur
    134       INTEGER :: an, mois, jour
    135       REAL :: secondes
    136 
    137       LOGICAL first,callinigrads
    138 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
    139       save first
    140       data first/.true./
    141       real dt_cum
    142       character*10 infile
    143       integer zan, tau0, thoriid
    144       integer nid_ctesGCM
    145       save nid_ctesGCM
    146       real degres
    147       real rlong(iip1), rlatg(jjp1)
    148       real zx_tmp_2d(iip1,jjp1)
    149       integer ndex2d(iip1*jjp1)
    150       logical ok_sync
    151       parameter (ok_sync = .true.)
    152       logical physic
    153 
    154       data callinigrads/.true./
    155       character*10 string10
    156 
    157       REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
    158 
    159 c+jld variables test conservation energie
    160       REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
    161 C     Tendance de la temp. potentiel d (theta)/ d t due a la
    162 C     tansformation d'energie cinetique en energie thermique
    163 C     cree par la dissipation
    164       REAL dtetaecdt(ip1jmp1,llm)
    165       REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
    166       REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
    167       REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
    168       CHARACTER*15 ztit
    169 !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
    170 !IM   SAVE      ip_ebil_dyn
    171 !IM   DATA      ip_ebil_dyn/0/
    172 c-jld
    173 
    174       character*80 dynhist_file, dynhistave_file
    175       character(len=*),parameter :: modname="leapfrog"
    176       character*80 abort_message
    177 
    178       logical dissip_conservative
    179       save dissip_conservative
    180       data dissip_conservative/.true./
    181 
    182       LOGICAL prem
    183       save prem
    184       DATA prem/.true./
    185       INTEGER testita
    186       PARAMETER (testita = 9)
    187 
    188       logical , parameter :: flag_verif = .false.
    189      
    190 
    191       integer itau_w   ! pas de temps ecriture = itap + itau_phy
    192 
    193 
    194       if (nday>=0) then
    195          itaufin   = nday*day_step
     4!
     5!
     6SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
     7
     8
     9  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
     10#ifdef CPP_IOIPSL
     11  use IOIPSL
     12#endif
     13  USE infotrac, ONLY: nqtot, isoCheck
     14  USE guide_mod, ONLY : guide_main
     15  USE write_field, ONLY: writefield
     16  USE control_mod, ONLY: nday, day_step, planet_type, offline, &
     17        iconser, iphysiq, iperiod, dissip_period, &
     18        iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins, &
     19        periodav, ok_dyn_ave, output_grads_dyn
     20  use exner_hyb_m, only: exner_hyb
     21  use exner_milieu_m, only: exner_milieu
     22  USE comvert_mod, ONLY: ap,bp,pressure_exner,presnivs
     23  USE comconst_mod, ONLY: cpp, dtphys, dtvr, pi, ihf
     24  USE logic_mod, ONLY: iflag_phys,ok_guide,forward,leapf,apphys, &
     25        statcl,conser,apdiss,purmats,ok_strato
     26  USE temps_mod, ONLY: jD_ref,jH_ref,itaufin,day_ini,day_ref, &
     27        start_time,dt
     28  USE strings_mod, ONLY: msg
     29
     30  IMPLICIT NONE
     31
     32   ! ......   Version  du 10/01/98    ..........
     33
     34   !        avec  coordonnees  verticales hybrides
     35  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
     36
     37  !=======================================================================
     38  !
     39  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
     40  !   -------
     41  !
     42  !   Objet:
     43  !   ------
     44  !
     45  !   GCM LMD nouvelle grille
     46  !
     47  !=======================================================================
     48  !
     49  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
     50  !  et possibilite d'appeler une fonction f(y)  a derivee tangente
     51  !  hyperbolique a la  place de la fonction a derivee sinusoidale.
     52
     53  !  ... Possibilite de choisir le shema pour l'advection de
     54  !    q  , en modifiant iadv dans traceur.def  (10/02) .
     55  !
     56  !  Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
     57  !  Pour Van-Leer iadv=10
     58  !
     59  !-----------------------------------------------------------------------
     60  !   Declarations:
     61  !   -------------
     62
     63  include "dimensions.h"
     64  include "paramet.h"
     65  include "comdissnew.h"
     66  include "comgeom.h"
     67  include "description.h"
     68  include "iniprint.h"
     69  include "academic.h"
     70
     71  REAL,INTENT(IN) :: time_0 ! not used
     72
     73  !   dynamical variables:
     74  REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm)    ! zonal covariant wind
     75  REAL,INTENT(INOUT) :: vcov(ip1jm,llm)      ! meridional covariant wind
     76  REAL,INTENT(INOUT) :: teta(ip1jmp1,llm)    ! potential temperature
     77  REAL,INTENT(INOUT) :: ps(ip1jmp1)          ! surface pressure (Pa)
     78  REAL,INTENT(INOUT) :: masse(ip1jmp1,llm)   ! air mass
     79  REAL,INTENT(INOUT) :: phis(ip1jmp1)        ! geopotentiat at the surface
     80  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
     81
     82  REAL :: p (ip1jmp1,llmp1  )               ! interlayer pressure
     83  REAL :: pks(ip1jmp1)                      ! exner at the surface
     84  REAL :: pk(ip1jmp1,llm)                   ! exner at mid-layer
     85  REAL :: pkf(ip1jmp1,llm)                  ! filtered exner at mid-layer
     86  REAL :: phi(ip1jmp1,llm)                  ! geopotential
     87  REAL :: w(ip1jmp1,llm)                    ! vertical velocity
     88
     89  real :: zqmin,zqmax
     90
     91  ! variables dynamiques intermediaire pour le transport
     92  REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
     93
     94  !   variables dynamiques au pas -1
     95  REAL :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
     96  REAL :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
     97  REAL :: massem1(ip1jmp1,llm)
     98
     99  !   tendances dynamiques
     100  REAL :: dv(ip1jm,llm),du(ip1jmp1,llm)
     101  REAL :: dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
     102
     103  !   tendances de la dissipation
     104  REAL :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
     105  REAL :: dtetadis(ip1jmp1,llm)
     106
     107  !   tendances physiques
     108  REAL :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
     109  REAL :: dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
     110
     111  !   variables pour le fichier histoire
     112  REAL :: dtav      ! intervalle de temps elementaire
     113
     114  REAL :: tppn(iim),tpps(iim),tpn,tps
     115  !
     116  INTEGER :: itau,itaufinp1,iav
     117   ! INTEGER  iday ! jour julien
     118  REAL :: time
     119
     120  REAL :: SSUM
     121  ! REAL finvmaold(ip1jmp1,llm)
     122
     123  !ym      LOGICAL  lafin
     124  LOGICAL :: lafin=.false.
     125  INTEGER :: ij,iq,l
     126  INTEGER :: ik
     127
     128  real :: time_step, t_wrt, t_ops
     129
     130   ! REAL rdayvrai,rdaym_ini
     131  ! jD_cur: jour julien courant
     132  ! jH_cur: heure julienne courante
     133  REAL :: jD_cur, jH_cur
     134  INTEGER :: an, mois, jour
     135  REAL :: secondes
     136
     137  LOGICAL :: first,callinigrads
     138  !IM : pour sortir les param. du modele dans un fis. netcdf 110106
     139  save first
     140  data first/.true./
     141  real :: dt_cum
     142  character(len=10) :: infile
     143  integer :: zan, tau0, thoriid
     144  integer :: nid_ctesGCM
     145  save nid_ctesGCM
     146  real :: degres
     147  real :: rlong(iip1), rlatg(jjp1)
     148  real :: zx_tmp_2d(iip1,jjp1)
     149  integer :: ndex2d(iip1*jjp1)
     150  logical :: ok_sync
     151  parameter (ok_sync = .true.)
     152  logical :: physic
     153
     154  data callinigrads/.true./
     155  character(len=10) :: string10
     156
     157  REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
     158
     159  !+jld variables test conservation energie
     160  REAL :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
     161  ! Tendance de la temp. potentiel d (theta)/ d t due a la
     162  ! tansformation d'energie cinetique en energie thermique
     163  ! cree par la dissipation
     164  REAL :: dtetaecdt(ip1jmp1,llm)
     165  REAL :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
     166  REAL :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
     167  REAL :: d_h_vcol, d_qt, d_qw, d_ql, d_ec
     168  CHARACTER(len=15) :: ztit
     169  !IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
     170  !IM   SAVE      ip_ebil_dyn
     171  !IM   DATA      ip_ebil_dyn/0/
     172  !-jld
     173
     174  character(len=80) :: dynhist_file, dynhistave_file
     175  character(len=*),parameter :: modname="leapfrog"
     176  character(len=80) :: abort_message
     177
     178  logical :: dissip_conservative
     179  save dissip_conservative
     180  data dissip_conservative/.true./
     181
     182  LOGICAL :: prem
     183  save prem
     184  DATA prem/.true./
     185  INTEGER :: testita
     186  PARAMETER (testita = 9)
     187
     188  logical , parameter :: flag_verif = .false.
     189
     190
     191  integer :: itau_w   ! pas de temps ecriture = itap + itau_phy
     192
     193
     194  if (nday>=0) then
     195     itaufin   = nday*day_step
     196  else
     197     itaufin   = -nday
     198  endif
     199  itaufinp1 = itaufin +1
     200  itau = 0
     201  physic=.true.
     202  if (iflag_phys==0.or.iflag_phys==2) physic=.false.
     203
     204   ! iday = day_ini+itau/day_step
     205   ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
     206   !    IF(time.GT.1.) THEN
     207   !     time = time-1.
     208   !     iday = iday+1
     209   !    ENDIF
     210
     211
     212  !-----------------------------------------------------------------------
     213  !   On initialise la pression et la fonction d'Exner :
     214  !   --------------------------------------------------
     215
     216  dq(:,:,:)=0.
     217  CALL pression ( ip1jmp1, ap, bp, ps, p       )
     218  if (pressure_exner) then
     219    CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
     220  else
     221    CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
     222  endif
     223
     224  !-----------------------------------------------------------------------
     225  !   Debut de l'integration temporelle:
     226  !   ----------------------------------
     227
     228   1   CONTINUE ! Matsuno Forward step begins here
     229
     230  !   date: (NB: date remains unchanged for Backward step)
     231  !   -----
     232
     233  jD_cur = jD_ref + day_ini - day_ref +                             &
     234        (itau+1)/day_step
     235  jH_cur = jH_ref + start_time +                                    &
     236        mod(itau+1,day_step)/float(day_step)
     237  jD_cur = jD_cur + int(jH_cur)
     238  jH_cur = jH_cur - int(jH_cur)
     239
     240  call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
     241
     242#ifdef CPP_IOIPSL
     243  if (ok_guide) then
     244    call guide_main(itau,ucov,vcov,teta,q,masse,ps)
     245  endif
     246#endif
     247
     248
     249  !
     250  ! IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
     251  !   CALL  test_period ( ucov,vcov,teta,q,p,phis )
     252  !   PRINT *,' ----   Test_period apres continue   OK ! -----', itau
     253  ! ENDIF
     254  !
     255
     256  ! Save fields obtained at previous time step as '...m1'
     257  CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
     258  CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
     259  CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
     260  CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
     261  CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
     262
     263  forward = .TRUE.
     264  leapf   = .FALSE.
     265  dt      =  dtvr
     266
     267  !   ...    P.Le Van .26/04/94  ....
     268  ! Ehouarn: finvmaold is actually not used
     269   ! CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
     270   ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
     271
     272  call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
     273
     274   2   CONTINUE ! Matsuno backward or leapfrog step begins here
     275
     276  !-----------------------------------------------------------------------
     277
     278  !   date: (NB: only leapfrog step requires recomputing date)
     279  !   -----
     280
     281  IF (leapf) THEN
     282    jD_cur = jD_ref + day_ini - day_ref + &
     283          (itau+1)/day_step
     284    jH_cur = jH_ref + start_time + &
     285          mod(itau+1,day_step)/float(day_step)
     286    jD_cur = jD_cur + int(jH_cur)
     287    jH_cur = jH_cur - int(jH_cur)
     288  ENDIF
     289
     290
     291  !   gestion des appels de la physique et des dissipations:
     292  !   ------------------------------------------------------
     293  !
     294  !   ...    P.Le Van  ( 6/02/95 )  ....
     295
     296  apphys = .FALSE.
     297  statcl = .FALSE.
     298  conser = .FALSE.
     299  apdiss = .FALSE.
     300
     301  IF( purmats ) THEN
     302  ! ! Purely Matsuno time stepping
     303     IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
     304     IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward ) &
     305           apdiss = .TRUE.
     306     IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward &
     307           .and. physic                        ) apphys = .TRUE.
     308  ELSE
     309  ! ! Leapfrog/Matsuno time stepping
     310     IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
     311     IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward ) &
     312           apdiss = .TRUE.
     313     IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
     314  END IF
     315
     316  ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
     317       ! supress dissipation step
     318  if (llm.eq.1) then
     319    apdiss=.false.
     320  endif
     321
     322
     323  call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
     324
     325  !-----------------------------------------------------------------------
     326  !   calcul des tendances dynamiques:
     327  !   --------------------------------
     328
     329  ! ! compute geopotential phi()
     330  CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
     331
     332  time = jD_cur + jH_cur
     333  CALL caldyn &
     334        ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis , &
     335        phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
     336
     337
     338  !-----------------------------------------------------------------------
     339  !   calcul des tendances advection des traceurs (dont l'humidite)
     340  !   -------------------------------------------------------------
     341
     342  call check_isotopes_seq(q,ip1jmp1, &
     343        'leapfrog 686: avant caladvtrac')
     344
     345  IF( forward.OR. leapf )  THEN
     346  ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
     347     CALL caladvtrac(q,pbaru,pbarv, &
     348           p, masse, dq,  teta, &
     349           flxw, pk)
     350      ! !write(*,*) 'caladvtrac 346'
     351
     352
     353     IF (offline) THEN
     354  !maf stokage du flux de masse pour traceurs OFF-LINE
     355
     356#ifdef CPP_IOIPSL
     357       CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis, &
     358             dtvr, itau)
     359#endif
     360
     361
     362     ENDIF ! of IF (offline)
     363  !
     364  ENDIF ! of IF( forward.OR. leapf )
     365
     366
     367  !-----------------------------------------------------------------------
     368  !   integrations dynamique et traceurs:
     369  !   ----------------------------------
     370
     371   CALL msg('720', modname, isoCheck)
     372   call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
     373
     374   CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 , &
     375         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
     376  ! $              finvmaold                                    )
     377
     378   CALL msg('724', modname, isoCheck)
     379   call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
     380
     381  ! .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
     382  !
     383  !-----------------------------------------------------------------------
     384  !   calcul des tendances physiques:
     385  !   -------------------------------
     386  !    ########   P.Le Van ( Modif le  6/02/95 )   ###########
     387  !
     388   IF( purmats )  THEN
     389      IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
     390   ELSE
     391      IF( itau+1.EQ. itaufin )              lafin = .TRUE.
     392   ENDIF
     393  !
     394  !
     395   IF( apphys )  THEN
     396  !
     397  ! .......   Ajout   P.Le Van ( 17/04/96 )   ...........
     398  !
     399
     400     CALL pression (  ip1jmp1, ap, bp, ps,  p      )
     401     if (pressure_exner) then
     402       CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
     403     else
     404       CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
     405     endif
     406
     407  ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
     408  ! avec dyn3dmem
     409     CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
     410
     411        ! rdaym_ini  = itau * dtvr / daysec
     412        ! rdayvrai   = rdaym_ini  + day_ini
     413        ! jD_cur = jD_ref + day_ini - day_ref
     414  ! $        + int (itau * dtvr / daysec)
     415  !       jH_cur = jH_ref +                                            &
     416  ! &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
     417       jD_cur = jD_ref + day_ini - day_ref +                        &
     418             (itau+1)/day_step
     419
     420       IF (planet_type .eq."generic") THEN
     421          ! ! AS: we make jD_cur to be pday
     422          jD_cur = int(day_ini + itau/day_step)
     423       ENDIF
     424
     425       jH_cur = jH_ref + start_time +                               &
     426             mod(itau+1,day_step)/float(day_step)
     427       jD_cur = jD_cur + int(jH_cur)
     428       jH_cur = jH_cur - int(jH_cur)
     429      ! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
     430      ! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
     431      ! write(lunout,*)'current date = ',an, mois, jour, secondes
     432
     433  ! rajout debug
     434    ! lafin = .true.
     435
     436
     437  !   Inbterface avec les routines de phylmd (phymars ... )
     438  !   -----------------------------------------------------
     439
     440  !+jld
     441
     442  !  Diagnostique de conservation de l'energie : initialisation
     443     IF (ip_ebil_dyn.ge.1 ) THEN
     444      ztit='bil dyn'
     445  ! Ehouarn: be careful, diagedyn is Earth-specific!
     446       IF (planet_type.eq."earth") THEN
     447        CALL diagedyn(ztit,2,1,1,dtphys &
     448              , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
     449       ENDIF
     450     ENDIF ! of IF (ip_ebil_dyn.ge.1 )
     451  !-jld
     452#ifdef CPP_IOIPSL
     453  !IM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
     454  !IM uncomment next 6 lines to get some parameters for LMDZ dynamics
     455     ! IF (first) THEN
     456     !  first=.false.
     457  !#include "ini_paramLMDZ_dyn.h"
     458     ! ENDIF
     459  !
     460  !#include "write_paramLMDZ_dyn.h"
     461  !
     462#endif
     463  ! #endif of #ifdef CPP_IOIPSL
     464#ifdef CPP_PHYS
     465     CALL calfis( lafin , jD_cur, jH_cur, &
     466           ucov,vcov,teta,q,masse,ps,p,pk,phis,phi , &
     467           du,dv,dteta,dq, &
     468           flxw,dufi,dvfi,dtetafi,dqfi,dpfi  )
     469#endif
     470   ! ajout des tendances physiques:
     471   ! ------------------------------
     472      CALL addfi( dtphys, leapf, forward   , &
     473            ucov, vcov, teta , q   ,ps , &
     474            dufi, dvfi, dtetafi , dqfi ,dpfi  )
     475      ! ! since addfi updates ps(), also update p(), masse() and pk()
     476      CALL pression (ip1jmp1,ap,bp,ps,p)
     477      CALL massdair(p,masse)
     478      if (pressure_exner) then
     479        CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
    196480      else
    197          itaufin   = -nday
     481        CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
    198482      endif
    199       itaufinp1 = itaufin +1
    200       itau = 0
    201       physic=.true.
    202       if (iflag_phys==0.or.iflag_phys==2) physic=.false.
    203 
    204 c      iday = day_ini+itau/day_step
    205 c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    206 c         IF(time.GT.1.) THEN
    207 c          time = time-1.
    208 c          iday = iday+1
    209 c         ENDIF
    210 
    211 
    212 c-----------------------------------------------------------------------
    213 c   On initialise la pression et la fonction d'Exner :
    214 c   --------------------------------------------------
    215 
    216       dq(:,:,:)=0.
    217       CALL pression ( ip1jmp1, ap, bp, ps, p       )
    218       if (pressure_exner) then
    219         CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    220       else
    221         CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    222       endif
    223 
    224 c-----------------------------------------------------------------------
    225 c   Debut de l'integration temporelle:
    226 c   ----------------------------------
    227 
    228    1  CONTINUE ! Matsuno Forward step begins here
    229 
    230 c   date: (NB: date remains unchanged for Backward step)
    231 c   -----
    232 
    233       jD_cur = jD_ref + day_ini - day_ref +                             &
    234      &          (itau+1)/day_step
    235       jH_cur = jH_ref + start_time +                                    &
    236      &          mod(itau+1,day_step)/float(day_step)
    237       jD_cur = jD_cur + int(jH_cur)
    238       jH_cur = jH_cur - int(jH_cur)
    239 
    240       call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
    241 
    242 #ifdef CPP_IOIPSL
    243       if (ok_guide) then
    244         call guide_main(itau,ucov,vcov,teta,q,masse,ps)
    245       endif
    246 #endif
    247 
    248 
    249 c
    250 c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
    251 c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
    252 c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
    253 c     ENDIF
    254 c
    255 
    256 ! Save fields obtained at previous time step as '...m1'
    257       CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
    258       CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
    259       CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
    260       CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
    261       CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
    262 
    263       forward = .TRUE.
    264       leapf   = .FALSE.
    265       dt      =  dtvr
    266 
    267 c   ...    P.Le Van .26/04/94  ....
    268 ! Ehouarn: finvmaold is actually not used
    269 !      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
    270 !      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
    271 
    272       call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
    273 
    274    2  CONTINUE ! Matsuno backward or leapfrog step begins here
    275 
    276 c-----------------------------------------------------------------------
    277 
    278 c   date: (NB: only leapfrog step requires recomputing date)
    279 c   -----
    280 
    281       IF (leapf) THEN
    282         jD_cur = jD_ref + day_ini - day_ref +
    283      &            (itau+1)/day_step
    284         jH_cur = jH_ref + start_time +
    285      &            mod(itau+1,day_step)/float(day_step)
    286         jD_cur = jD_cur + int(jH_cur)
    287         jH_cur = jH_cur - int(jH_cur)
     483
     484     IF (ok_strato) THEN
     485       CALL top_bound( vcov,ucov,teta,masse,dtphys)
     486     ENDIF
     487
     488  !
     489  !  Diagnostique de conservation de l'energie : difference
     490     IF (ip_ebil_dyn.ge.1 ) THEN
     491      ztit='bil phys'
     492      IF (planet_type.eq."earth") THEN
     493       CALL diagedyn(ztit,2,1,1,dtphys &
     494             , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
    288495      ENDIF
    289 
    290 
    291 c   gestion des appels de la physique et des dissipations:
    292 c   ------------------------------------------------------
    293 c
    294 c   ...    P.Le Van  ( 6/02/95 )  ....
    295 
    296       apphys = .FALSE.
    297       statcl = .FALSE.
    298       conser = .FALSE.
    299       apdiss = .FALSE.
    300 
    301       IF( purmats ) THEN
    302       ! Purely Matsuno time stepping
    303          IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
    304          IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
    305      s        apdiss = .TRUE.
    306          IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
    307      s          .and. physic                        ) apphys = .TRUE.
    308       ELSE
    309       ! Leapfrog/Matsuno time stepping
    310          IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
    311          IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
    312      s        apdiss = .TRUE.
    313          IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
    314       END IF
    315 
    316 ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
    317 !          supress dissipation step
    318       if (llm.eq.1) then
    319         apdiss=.false.
    320       endif
    321 
    322 
    323       call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
    324 
    325 c-----------------------------------------------------------------------
    326 c   calcul des tendances dynamiques:
    327 c   --------------------------------
    328 
    329       ! compute geopotential phi()
    330       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    331 
    332       time = jD_cur + jH_cur
    333       CALL caldyn
    334      $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
    335      $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
    336 
    337 
    338 c-----------------------------------------------------------------------
    339 c   calcul des tendances advection des traceurs (dont l'humidite)
    340 c   -------------------------------------------------------------
    341 
    342       call check_isotopes_seq(q,ip1jmp1,
    343      &           'leapfrog 686: avant caladvtrac')
    344 
    345       IF( forward. OR . leapf )  THEN
    346 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
    347          CALL caladvtrac(q,pbaru,pbarv,
    348      *        p, masse, dq,  teta,
    349      .        flxw, pk)
    350           !write(*,*) 'caladvtrac 346'
    351 
    352          
    353          IF (offline) THEN
    354 Cmaf stokage du flux de masse pour traceurs OFF-LINE
    355 
    356 #ifdef CPP_IOIPSL
    357            CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
    358      .   dtvr, itau)
    359 #endif
    360 
    361 
    362          ENDIF ! of IF (offline)
    363 c
    364       ENDIF ! of IF( forward. OR . leapf )
    365 
    366 
    367 c-----------------------------------------------------------------------
    368 c   integrations dynamique et traceurs:
    369 c   ----------------------------------
    370 
    371        CALL msg('720', modname, isoCheck)
    372        call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
    373        
    374        CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
    375      $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
    376 !     $              finvmaold                                    )
    377 
    378        CALL msg('724', modname, isoCheck)
    379        call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
    380 
    381 c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
    382 c
    383 c-----------------------------------------------------------------------
    384 c   calcul des tendances physiques:
    385 c   -------------------------------
    386 c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
    387 c
    388        IF( purmats )  THEN
    389           IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
    390        ELSE
    391           IF( itau+1. EQ. itaufin )              lafin = .TRUE.
    392        ENDIF
    393 c
    394 c
    395        IF( apphys )  THEN
    396 c
    397 c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
    398 c
    399 
    400          CALL pression (  ip1jmp1, ap, bp, ps,  p      )
    401          if (pressure_exner) then
    402            CALL exner_hyb(  ip1jmp1, ps, p,pks, pk, pkf )
    403          else
    404            CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    405          endif
    406 
    407 ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
    408 ! avec dyn3dmem
    409          CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
    410 
    411 !           rdaym_ini  = itau * dtvr / daysec
    412 !           rdayvrai   = rdaym_ini  + day_ini
    413 !           jD_cur = jD_ref + day_ini - day_ref
    414 !     $        + int (itau * dtvr / daysec)
    415 !           jH_cur = jH_ref +                                            &
    416 !     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
    417            jD_cur = jD_ref + day_ini - day_ref +                        &
    418      &          (itau+1)/day_step
    419 
    420            IF (planet_type .eq."generic") THEN
    421               ! AS: we make jD_cur to be pday
    422               jD_cur = int(day_ini + itau/day_step)
     496     ENDIF ! of IF (ip_ebil_dyn.ge.1 )
     497
     498   ENDIF ! of IF( apphys )
     499
     500  IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
     501  !   Academic case : Simple friction and Newtonan relaxation
     502  !   -------------------------------------------------------
     503    DO l=1,llm
     504      DO ij=1,ip1jmp1
     505       teta(ij,l)=teta(ij,l)-dtvr* &
     506             (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
     507      ENDDO
     508    ENDDO ! of DO l=1,llm
     509
     510    if (planet_type.eq."giant") then
     511      ! ! add an intrinsic heat flux at the base of the atmosphere
     512      teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
     513    endif
     514
     515    call friction(ucov,vcov,dtvr)
     516
     517    ! ! Sponge layer (if any)
     518    IF (ok_strato) THEN
     519       ! dufi(:,:)=0.
     520       ! dvfi(:,:)=0.
     521       ! dtetafi(:,:)=0.
     522       ! dqfi(:,:,:)=0.
     523  !          dpfi(:)=0.
     524       ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     525       CALL top_bound( vcov,ucov,teta,masse,dtvr)
     526       ! CALL addfi( dtvr, leapf, forward   ,
     527  ! $                  ucov, vcov, teta , q   ,ps ,
     528  ! $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     529    ENDIF ! of IF (ok_strato)
     530  ENDIF ! of IF (iflag_phys.EQ.2)
     531
     532
     533  !-jld
     534
     535    CALL pression ( ip1jmp1, ap, bp, ps, p                  )
     536    if (pressure_exner) then
     537      CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
     538    else
     539      CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
     540    endif
     541    CALL massdair(p,masse)
     542
     543    call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
     544
     545  !-----------------------------------------------------------------------
     546  !   dissipation horizontale et verticale  des petites echelles:
     547  !   ----------------------------------------------------------
     548
     549  IF(apdiss) THEN
     550
     551
     552  !   calcul de l'energie cinetique avant dissipation
     553    call covcont(llm,ucov,vcov,ucont,vcont)
     554    call enercin(vcov,ucov,vcont,ucont,ecin0)
     555
     556  !   dissipation
     557    CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
     558    ucov=ucov+dudis
     559    vcov=vcov+dvdis
     560    ! teta=teta+dtetadis
     561
     562
     563  !------------------------------------------------------------------------
     564    if (dissip_conservative) then
     565    ! On rajoute la tendance due a la transform. Ec -> E therm. cree
     566    ! lors de la dissipation
     567        call covcont(llm,ucov,vcov,ucont,vcont)
     568        call enercin(vcov,ucov,vcont,ucont,ecin)
     569        dtetaecdt= (ecin0-ecin)/ pk
     570        ! teta=teta+dtetaecdt
     571        dtetadis=dtetadis+dtetaecdt
     572    endif
     573    teta=teta+dtetadis
     574  !------------------------------------------------------------------------
     575
     576
     577  !    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
     578  !   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
     579  !
     580
     581    DO l  =  1, llm
     582      DO ij =  1,iim
     583       tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
     584       tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
     585      ENDDO
     586       tpn  = SSUM(iim,tppn,1)/apoln
     587       tps  = SSUM(iim,tpps,1)/apols
     588
     589      DO ij = 1, iip1
     590       teta(  ij    ,l) = tpn
     591       teta(ij+ip1jm,l) = tps
     592      ENDDO
     593    ENDDO
     594
     595    if (1 == 0) then
     596  !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
     597  !!!                     2) should probably not be here anyway
     598  !!! but are kept for those who would want to revert to previous behaviour
     599       DO ij =  1,iim
     600         tppn(ij)  = aire(  ij    ) * ps (  ij    )
     601         tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
     602       ENDDO
     603         tpn  = SSUM(iim,tppn,1)/apoln
     604         tps  = SSUM(iim,tpps,1)/apols
     605
     606       DO ij = 1, iip1
     607         ps(  ij    ) = tpn
     608         ps(ij+ip1jm) = tps
     609       ENDDO
     610    endif ! of if (1 == 0)
     611
     612  END IF ! of IF(apdiss)
     613
     614  ! ajout debug
     615           ! IF( lafin ) then
     616           !   abort_message = 'Simulation finished'
     617           !   call abort_gcm(modname,abort_message,0)
     618           ! ENDIF
     619
     620  !   ********************************************************************
     621  !   ********************************************************************
     622  !   .... fin de l'integration dynamique  et physique pour le pas itau ..
     623  !   ********************************************************************
     624  !   ********************************************************************
     625
     626  !   preparation du pas d'integration suivant  ......
     627
     628  call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
     629
     630  IF ( .NOT.purmats ) THEN
     631    ! ........................................................
     632    ! ..............  schema matsuno + leapfrog  ..............
     633    ! ........................................................
     634
     635        IF(forward.OR. leapf) THEN
     636          itau= itau + 1
     637           ! iday= day_ini+itau/day_step
     638           ! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
     639           !   IF(time.GT.1.) THEN
     640           !     time = time-1.
     641           !     iday = iday+1
     642           !   ENDIF
     643        ENDIF
     644
     645
     646        IF( itau.EQ. itaufinp1 ) then
     647          if (flag_verif) then
     648            write(79,*) 'ucov',ucov
     649            write(80,*) 'vcov',vcov
     650            write(81,*) 'teta',teta
     651            write(82,*) 'ps',ps
     652            write(83,*) 'q',q
     653            WRITE(85,*) 'q1 = ',q(:,:,1)
     654            WRITE(86,*) 'q3 = ',q(:,:,3)
     655          endif
     656
     657          abort_message = 'Simulation finished'
     658
     659          call abort_gcm(modname,abort_message,0)
     660        ENDIF
     661  !-----------------------------------------------------------------------
     662  !   ecriture du fichier histoire moyenne:
     663  !   -------------------------------------
     664
     665        IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
     666           IF(itau.EQ.itaufin) THEN
     667              iav=1
     668           ELSE
     669              iav=0
    423670           ENDIF
    424671
    425            jH_cur = jH_ref + start_time +                               &
    426      &              mod(itau+1,day_step)/float(day_step)
    427            jD_cur = jD_cur + int(jH_cur)
    428            jH_cur = jH_cur - int(jH_cur)
    429 !         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
    430 !         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
    431 !         write(lunout,*)'current date = ',an, mois, jour, secondes
    432 
    433 c rajout debug
    434 c       lafin = .true.
    435 
    436 
    437 c   Inbterface avec les routines de phylmd (phymars ... )
    438 c   -----------------------------------------------------
    439 
    440 c+jld
    441 
    442 c  Diagnostique de conservation de l'energie : initialisation
    443          IF (ip_ebil_dyn.ge.1 ) THEN
    444           ztit='bil dyn'
    445 ! Ehouarn: be careful, diagedyn is Earth-specific!
    446            IF (planet_type.eq."earth") THEN
    447             CALL diagedyn(ztit,2,1,1,dtphys
    448      &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
     672           ! ! Ehouarn: re-compute geopotential for outputs
     673           CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     674
     675           IF (ok_dynzon) THEN
     676#ifdef CPP_IOIPSL
     677             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
     678                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
     679#endif
     680           END IF
     681           IF (ok_dyn_ave) THEN
     682#ifdef CPP_IOIPSL
     683             CALL writedynav(itau,vcov, &
     684                   ucov,teta,pk,phi,q,masse,ps,phis)
     685#endif
    449686           ENDIF
    450          ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    451 c-jld
    452 #ifdef CPP_IOIPSL
    453 cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
    454 cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
    455 c        IF (first) THEN
    456 c         first=.false.
    457 c#include "ini_paramLMDZ_dyn.h"
    458 c        ENDIF
    459 c
    460 c#include "write_paramLMDZ_dyn.h"
    461 c
    462 #endif
    463 ! #endif of #ifdef CPP_IOIPSL
    464 #ifdef CPP_PHYS
    465          CALL calfis( lafin , jD_cur, jH_cur,
    466      $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
    467      $               du,dv,dteta,dq,
    468      $               flxw,dufi,dvfi,dtetafi,dqfi,dpfi  )
    469 #endif
    470 c      ajout des tendances physiques:
    471 c      ------------------------------
    472           CALL addfi( dtphys, leapf, forward   ,
    473      $                  ucov, vcov, teta , q   ,ps ,
    474      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    475           ! since addfi updates ps(), also update p(), masse() and pk()
    476           CALL pression (ip1jmp1,ap,bp,ps,p)
    477           CALL massdair(p,masse)
    478           if (pressure_exner) then
    479             CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
    480           else
    481             CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
     687
     688        ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
     689
     690        call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
     691
     692  !-----------------------------------------------------------------------
     693  !   ecriture de la bande histoire:
     694  !   ------------------------------
     695
     696        IF( MOD(itau,iecri).EQ.0) THEN
     697         ! ! Ehouarn: output only during LF or Backward Matsuno
     698         if (leapf.or.(.not.leapf.and.(.not.forward))) then
     699          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     700          unat=0.
     701          do l=1,llm
     702            unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
     703            vnat(:,l)=vcov(:,l)/cv(:)
     704          enddo
     705#ifdef CPP_IOIPSL
     706          if (ok_dyn_ins) then
     707            ! write(lunout,*) "leapfrog: call writehist, itau=",itau
     708           CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
     709            ! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
     710            ! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
     711           ! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
     712           !  call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
     713           !  call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
     714          endif ! of if (ok_dyn_ins)
     715#endif
     716  ! For some Grads outputs of fields
     717          if (output_grads_dyn) then
     718#include "write_grads_dyn.h"
    482719          endif
    483 
    484          IF (ok_strato) THEN
    485            CALL top_bound( vcov,ucov,teta,masse,dtphys)
    486          ENDIF
    487        
    488 c
    489 c  Diagnostique de conservation de l'energie : difference
    490          IF (ip_ebil_dyn.ge.1 ) THEN
    491           ztit='bil phys'
    492           IF (planet_type.eq."earth") THEN
    493            CALL diagedyn(ztit,2,1,1,dtphys
    494      &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
    495           ENDIF
    496          ENDIF ! of IF (ip_ebil_dyn.ge.1 )
    497 
    498        ENDIF ! of IF( apphys )
    499 
    500       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
    501 !   Academic case : Simple friction and Newtonan relaxation
    502 !   -------------------------------------------------------
    503         DO l=1,llm   
    504           DO ij=1,ip1jmp1
    505            teta(ij,l)=teta(ij,l)-dtvr*
    506      &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
    507           ENDDO
    508         ENDDO ! of DO l=1,llm
    509        
    510         if (planet_type.eq."giant") then
    511           ! add an intrinsic heat flux at the base of the atmosphere
    512           teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
    513         endif
    514 
    515         call friction(ucov,vcov,dtvr)
    516        
    517         ! Sponge layer (if any)
    518         IF (ok_strato) THEN
    519 !          dufi(:,:)=0.
    520 !          dvfi(:,:)=0.
    521 !          dtetafi(:,:)=0.
    522 !          dqfi(:,:,:)=0.
    523 !          dpfi(:)=0.
    524 !          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    525            CALL top_bound( vcov,ucov,teta,masse,dtvr)
    526 !          CALL addfi( dtvr, leapf, forward   ,
    527 !     $                  ucov, vcov, teta , q   ,ps ,
    528 !     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    529         ENDIF ! of IF (ok_strato)
    530       ENDIF ! of IF (iflag_phys.EQ.2)
    531 
    532 
    533 c-jld
    534 
    535         CALL pression ( ip1jmp1, ap, bp, ps, p                  )
    536         if (pressure_exner) then
    537           CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
    538         else
    539           CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
    540         endif
    541         CALL massdair(p,masse)
    542 
    543         call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
    544 
    545 c-----------------------------------------------------------------------
    546 c   dissipation horizontale et verticale  des petites echelles:
    547 c   ----------------------------------------------------------
    548 
    549       IF(apdiss) THEN
    550 
    551 
    552 c   calcul de l'energie cinetique avant dissipation
    553         call covcont(llm,ucov,vcov,ucont,vcont)
    554         call enercin(vcov,ucov,vcont,ucont,ecin0)
    555 
    556 c   dissipation
    557         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
    558         ucov=ucov+dudis
    559         vcov=vcov+dvdis
    560 c       teta=teta+dtetadis
    561 
    562 
    563 c------------------------------------------------------------------------
    564         if (dissip_conservative) then
    565 C       On rajoute la tendance due a la transform. Ec -> E therm. cree
    566 C       lors de la dissipation
    567             call covcont(llm,ucov,vcov,ucont,vcont)
    568             call enercin(vcov,ucov,vcont,ucont,ecin)
    569             dtetaecdt= (ecin0-ecin)/ pk
    570 c           teta=teta+dtetaecdt
    571             dtetadis=dtetadis+dtetaecdt
    572         endif
    573         teta=teta+dtetadis
    574 c------------------------------------------------------------------------
    575 
    576 
    577 c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
    578 c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
    579 c
    580 
    581         DO l  =  1, llm
    582           DO ij =  1,iim
    583            tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
    584            tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
    585           ENDDO
    586            tpn  = SSUM(iim,tppn,1)/apoln
    587            tps  = SSUM(iim,tpps,1)/apols
    588 
    589           DO ij = 1, iip1
    590            teta(  ij    ,l) = tpn
    591            teta(ij+ip1jm,l) = tps
    592           ENDDO
    593         ENDDO
    594 
    595         if (1 == 0) then
    596 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
    597 !!!                     2) should probably not be here anyway
    598 !!! but are kept for those who would want to revert to previous behaviour
    599            DO ij =  1,iim
    600              tppn(ij)  = aire(  ij    ) * ps (  ij    )
    601              tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
    602            ENDDO
    603              tpn  = SSUM(iim,tppn,1)/apoln
    604              tps  = SSUM(iim,tpps,1)/apols
    605 
    606            DO ij = 1, iip1
    607              ps(  ij    ) = tpn
    608              ps(ij+ip1jm) = tps
    609            ENDDO
    610         endif ! of if (1 == 0)
    611 
    612       END IF ! of IF(apdiss)
    613 
    614 c ajout debug
    615 c              IF( lafin ) then 
    616 c                abort_message = 'Simulation finished'
    617 c                call abort_gcm(modname,abort_message,0)
    618 c              ENDIF
    619        
    620 c   ********************************************************************
    621 c   ********************************************************************
    622 c   .... fin de l'integration dynamique  et physique pour le pas itau ..
    623 c   ********************************************************************
    624 c   ********************************************************************
    625 
    626 c   preparation du pas d'integration suivant  ......
    627 
    628       call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
    629 
    630       IF ( .NOT.purmats ) THEN
    631 c       ........................................................
    632 c       ..............  schema matsuno + leapfrog  ..............
    633 c       ........................................................
    634 
    635             IF(forward. OR. leapf) THEN
    636               itau= itau + 1
    637 c              iday= day_ini+itau/day_step
    638 c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    639 c                IF(time.GT.1.) THEN
    640 c                  time = time-1.
    641 c                  iday = iday+1
    642 c                ENDIF
    643             ENDIF
    644 
    645 
    646             IF( itau. EQ. itaufinp1 ) then 
    647               if (flag_verif) then
    648                 write(79,*) 'ucov',ucov
    649                 write(80,*) 'vcov',vcov
    650                 write(81,*) 'teta',teta
    651                 write(82,*) 'ps',ps
    652                 write(83,*) 'q',q
    653                 WRITE(85,*) 'q1 = ',q(:,:,1)
    654                 WRITE(86,*) 'q3 = ',q(:,:,3)
    655               endif
    656 
    657               abort_message = 'Simulation finished'
    658 
    659               call abort_gcm(modname,abort_message,0)
    660             ENDIF
    661 c-----------------------------------------------------------------------
    662 c   ecriture du fichier histoire moyenne:
    663 c   -------------------------------------
    664 
    665             IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
    666                IF(itau.EQ.itaufin) THEN
    667                   iav=1
     720         endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
     721        ENDIF ! of IF(MOD(itau,iecri).EQ.0)
     722
     723        IF(itau.EQ.itaufin) THEN
     724
     725
     726           ! if (planet_type.eq."earth") then
     727  ! Write an Earth-format restart file
     728            CALL dynredem1("restart.nc",start_time, &
     729                  vcov,ucov,teta,q,masse,ps)
     730           ! endif ! of if (planet_type.eq."earth")
     731
     732          CLOSE(99)
     733          if (ok_guide) then
     734            ! ! set ok_guide to false to avoid extra output
     735            ! ! in following forward step
     736            ok_guide=.false.
     737          endif
     738          ! !!! Ehouarn: Why not stop here and now?
     739        ENDIF ! of IF (itau.EQ.itaufin)
     740
     741  !-----------------------------------------------------------------------
     742  !   gestion de l'integration temporelle:
     743  !   ------------------------------------
     744
     745        IF( MOD(itau,iperiod).EQ.0 )    THEN
     746                GO TO 1
     747        ELSE IF ( MOD(itau-1,iperiod).EQ. 0 ) THEN
     748
     749               IF( forward )  THEN
     750   ! fin du pas forward et debut du pas backward
     751
     752                  forward = .FALSE.
     753                    leapf = .FALSE.
     754                       GO TO 2
     755
    668756               ELSE
    669                   iav=0
    670                ENDIF
    671                
    672 !              ! Ehouarn: re-compute geopotential for outputs
    673                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    674 
    675                IF (ok_dynzon) THEN
    676 #ifdef CPP_IOIPSL
    677                  CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
    678      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
    679 #endif
    680                END IF
    681                IF (ok_dyn_ave) THEN
    682 #ifdef CPP_IOIPSL
    683                  CALL writedynav(itau,vcov,
    684      &                 ucov,teta,pk,phi,q,masse,ps,phis)
    685 #endif
    686                ENDIF
    687 
    688             ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
    689 
    690             call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
    691 
    692 c-----------------------------------------------------------------------
    693 c   ecriture de la bande histoire:
    694 c   ------------------------------
    695 
    696             IF( MOD(itau,iecri).EQ.0) THEN
    697              ! Ehouarn: output only during LF or Backward Matsuno
    698              if (leapf.or.(.not.leapf.and.(.not.forward))) then
    699               CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    700               unat=0.
    701               do l=1,llm
    702                 unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
    703                 vnat(:,l)=vcov(:,l)/cv(:)
    704               enddo
    705 #ifdef CPP_IOIPSL
    706               if (ok_dyn_ins) then
    707 !               write(lunout,*) "leapfrog: call writehist, itau=",itau
    708                CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
    709 !               call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
    710 !               call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
    711 !              call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
    712 !               call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
    713 !               call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
    714               endif ! of if (ok_dyn_ins)
    715 #endif
    716 ! For some Grads outputs of fields
    717               if (output_grads_dyn) then
     757   ! fin du pas backward et debut du premier pas leapfrog
     758
     759                    leapf =  .TRUE.
     760                    dt  =  2.*dtvr
     761                    GO TO 2
     762               END IF ! of IF (forward)
     763        ELSE
     764
     765   ! ......   pas leapfrog  .....
     766
     767             leapf = .TRUE.
     768             dt  = 2.*dtvr
     769             GO TO 2
     770        END IF ! of IF (MOD(itau,iperiod).EQ.0)
     771               ! !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
     772
     773  ELSE ! of IF (.not.purmats)
     774
     775        call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
     776
     777    ! ........................................................
     778    ! ..............       schema  matsuno        ...............
     779    ! ........................................................
     780        IF( forward )  THEN
     781
     782         itau =  itau + 1
     783          ! iday = day_ini+itau/day_step
     784          ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
     785  !
     786  !              IF(time.GT.1.) THEN
     787  !               time = time-1.
     788  !               iday = iday+1
     789  !              ENDIF
     790
     791           forward =  .FALSE.
     792           IF( itau.EQ. itaufinp1 ) then
     793             abort_message = 'Simulation finished'
     794             call abort_gcm(modname,abort_message,0)
     795           ENDIF
     796           GO TO 2
     797
     798        ELSE ! of IF(forward) i.e. backward step
     799
     800          call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
     801
     802          IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
     803           IF(itau.EQ.itaufin) THEN
     804              iav=1
     805           ELSE
     806              iav=0
     807           ENDIF
     808
     809           ! ! Ehouarn: re-compute geopotential for outputs
     810           CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     811
     812           IF (ok_dynzon) THEN
     813#ifdef CPP_IOIPSL
     814             CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav, &
     815                   ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
     816#endif
     817           ENDIF
     818           IF (ok_dyn_ave) THEN
     819#ifdef CPP_IOIPSL
     820             CALL writedynav(itau,vcov, &
     821                   ucov,teta,pk,phi,q,masse,ps,phis)
     822#endif
     823           ENDIF
     824
     825          ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
     826
     827          IF(MOD(itau,iecri         ).EQ.0) THEN
     828           ! IF(MOD(itau,iecri*day_step).EQ.0) THEN
     829            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
     830            unat=0.
     831            do l=1,llm
     832              unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
     833              vnat(:,l)=vcov(:,l)/cv(:)
     834            enddo
     835#ifdef CPP_IOIPSL
     836          if (ok_dyn_ins) then
     837             ! write(lunout,*) "leapfrog: call writehist (b)",
     838  ! &                        itau,iecri
     839            CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
     840          endif ! of if (ok_dyn_ins)
     841#endif
     842  ! For some Grads outputs
     843            if (output_grads_dyn) then
    718844#include "write_grads_dyn.h"
    719               endif
    720              endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
    721             ENDIF ! of IF(MOD(itau,iecri).EQ.0)
    722 
    723             IF(itau.EQ.itaufin) THEN
    724 
    725 
    726 !              if (planet_type.eq."earth") then
    727 ! Write an Earth-format restart file
    728                 CALL dynredem1("restart.nc",start_time,
    729      &                         vcov,ucov,teta,q,masse,ps)
    730 !              endif ! of if (planet_type.eq."earth")
    731 
    732               CLOSE(99)
    733               if (ok_guide) then
    734                 ! set ok_guide to false to avoid extra output
    735                 ! in following forward step
    736                 ok_guide=.false.
    737               endif
    738               !!! Ehouarn: Why not stop here and now?
    739             ENDIF ! of IF (itau.EQ.itaufin)
    740 
    741 c-----------------------------------------------------------------------
    742 c   gestion de l'integration temporelle:
    743 c   ------------------------------------
    744 
    745             IF( MOD(itau,iperiod).EQ.0 )    THEN
    746                     GO TO 1
    747             ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
    748 
    749                    IF( forward )  THEN
    750 c      fin du pas forward et debut du pas backward
    751 
    752                       forward = .FALSE.
    753                         leapf = .FALSE.
    754                            GO TO 2
    755 
    756                    ELSE
    757 c      fin du pas backward et debut du premier pas leapfrog
    758 
    759                         leapf =  .TRUE.
    760                         dt  =  2.*dtvr
    761                         GO TO 2
    762                    END IF ! of IF (forward)
    763             ELSE
    764 
    765 c      ......   pas leapfrog  .....
    766 
    767                  leapf = .TRUE.
    768                  dt  = 2.*dtvr
    769                  GO TO 2
    770             END IF ! of IF (MOD(itau,iperiod).EQ.0)
    771                    !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
    772 
    773       ELSE ! of IF (.not.purmats)
    774 
    775             call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
    776 
    777 c       ........................................................
    778 c       ..............       schema  matsuno        ...............
    779 c       ........................................................
    780             IF( forward )  THEN
    781 
    782              itau =  itau + 1
    783 c             iday = day_ini+itau/day_step
    784 c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
    785 c
    786 c                  IF(time.GT.1.) THEN
    787 c                   time = time-1.
    788 c                   iday = iday+1
    789 c                  ENDIF
    790 
    791                forward =  .FALSE.
    792                IF( itau. EQ. itaufinp1 ) then 
    793                  abort_message = 'Simulation finished'
    794                  call abort_gcm(modname,abort_message,0)
    795                ENDIF
    796                GO TO 2
    797 
    798             ELSE ! of IF(forward) i.e. backward step
    799  
    800               call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
    801 
    802               IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
    803                IF(itau.EQ.itaufin) THEN
    804                   iav=1
    805                ELSE
    806                   iav=0
    807                ENDIF
    808 
    809 !              ! Ehouarn: re-compute geopotential for outputs
    810                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    811 
    812                IF (ok_dynzon) THEN
    813 #ifdef CPP_IOIPSL
    814                  CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
    815      &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
    816 #endif
    817                ENDIF
    818                IF (ok_dyn_ave) THEN
    819 #ifdef CPP_IOIPSL
    820                  CALL writedynav(itau,vcov,
    821      &                 ucov,teta,pk,phi,q,masse,ps,phis)
    822 #endif
    823                ENDIF
    824 
    825               ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
    826 
    827               IF(MOD(itau,iecri         ).EQ.0) THEN
    828 c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
    829                 CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    830                 unat=0.
    831                 do l=1,llm
    832                   unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
    833                   vnat(:,l)=vcov(:,l)/cv(:)
    834                 enddo
    835 #ifdef CPP_IOIPSL
    836               if (ok_dyn_ins) then
    837 !                write(lunout,*) "leapfrog: call writehist (b)",
    838 !     &                        itau,iecri
    839                 CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
    840               endif ! of if (ok_dyn_ins)
    841 #endif
    842 ! For some Grads outputs
    843                 if (output_grads_dyn) then
    844 #include "write_grads_dyn.h"
    845                 endif
    846 
    847               ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
    848 
    849               IF(itau.EQ.itaufin) THEN
    850 !                if (planet_type.eq."earth") then
    851                   CALL dynredem1("restart.nc",start_time,
    852      &                           vcov,ucov,teta,q,masse,ps)
    853 !                endif ! of if (planet_type.eq."earth")
    854                 if (ok_guide) then
    855                   ! set ok_guide to false to avoid extra output
    856                   ! in following forward step
    857                   ok_guide=.false.
    858                 endif
    859               ENDIF ! of IF(itau.EQ.itaufin)
    860 
    861               forward = .TRUE.
    862               GO TO  1
    863 
    864             ENDIF ! of IF (forward)
    865 
    866       END IF ! of IF(.not.purmats)
    867 
    868       END
     845            endif
     846
     847          ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
     848
     849          IF(itau.EQ.itaufin) THEN
     850             ! if (planet_type.eq."earth") then
     851              CALL dynredem1("restart.nc",start_time, &
     852                    vcov,ucov,teta,q,masse,ps)
     853             ! endif ! of if (planet_type.eq."earth")
     854            if (ok_guide) then
     855              ! ! set ok_guide to false to avoid extra output
     856              ! ! in following forward step
     857              ok_guide=.false.
     858            endif
     859          ENDIF ! of IF(itau.EQ.itaufin)
     860
     861          forward = .TRUE.
     862          GO TO  1
     863
     864        ENDIF ! of IF (forward)
     865
     866  END IF ! of IF(.not.purmats)
     867
     868END SUBROUTINE leapfrog
Note: See TracChangeset for help on using the changeset viewer.