source: LMDZ.3.3/trunk/libf/dyn3d/gcm.F @ 507

Last change on this file since 507 was 344, checked in by lmdz, 23 years ago

Inclusion des modifs de D. Hauglustaine pour la version 1 de INCA
LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 21.8 KB
Line 
1
2      PROGRAM gcm
3
4#ifdef INCA
5      USE transport_controls, ONLY : adv_flg, mmt_adj
6#endif 
7      USE IOIPSL
8
9      IMPLICIT NONE
10
11c      ......   Version  du 10/01/98    ..........
12
13c             avec  coordonnees  verticales hybrides
14c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
15
16c=======================================================================
17c
18c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
19c   -------
20c
21c   Objet:
22c   ------
23c
24c   GCM LMD nouvelle grille
25c
26c=======================================================================
27c
28c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
29c      et possibilite d'appeler une fonction f(y)  a derivee tangente
30c      hyperbolique a la  place de la fonction a derivee sinusoidale.
31
32c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
33c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
34c
35c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
36c
37c-----------------------------------------------------------------------
38c   Declarations:
39c   -------------
40
41#include "dimensions.h"
42#include "paramet.h"
43#include "comconst.h"
44#include "comdissnew.h"
45#include "comvert.h"
46#include "comgeom.h"
47#include "logic.h"
48#include "temps.h"
49#include "control.h"
50#include "ener.h"
51#include "netcdf.inc"
52#include "description.h"
53#include "serre.h"
54#include "tracstoke.h"
55
56
57      INTEGER         longcles
58      PARAMETER     ( longcles = 20 )
59      REAL  clesphy0( longcles )
60      SAVE  clesphy0
61
62      INTEGER*4  iday ! jour julien
63      REAL       time ! Heure de la journee en fraction d'1 jour
64      REAL zdtvr
65      INTEGER nbetatmoy, nbetatdem,nbetat
66
67c   variables dynamiques
68      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
69      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
70      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
71      REAL ps(ip1jmp1)                       ! pression  au sol
72      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
73      REAL pks(ip1jmp1)                      ! exner au  sol
74      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
75      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
76      REAL masse(ip1jmp1,llm)                ! masse d'air
77      REAL phis(ip1jmp1)                     ! geopotentiel au sol
78      REAL phi(ip1jmp1,llm)                  ! geopotentiel
79      REAL w(ip1jmp1,llm)                    ! vitesse verticale
80
81c variables dynamiques intermediaire pour le transport
82      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
83
84c   variables dynamiques au pas -1
85      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
86      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
87      REAL massem1(ip1jmp1,llm)
88
89c   tendances dynamiques
90      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
91      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
92
93c   tendances de la dissipation
94      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
95      REAL dhdis(ip1jmp1,llm)
96
97c   tendances physiques
98      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
99      REAL dhfi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
100
101c   variables pour le fichier histoire
102      REAL dtav      ! intervalle de temps elementaire
103
104      REAL tppn(iim),tpps(iim),tpn,tps
105c
106      INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
107
108      INTEGER itau,itaufinp1,iav
109
110      EXTERNAL caldyn, traceur
111      EXTERNAL dissip,geopot,iniconst,inifilr
112      EXTERNAL integrd,SCOPY
113      EXTERNAL iniav,writeav,writehist
114      EXTERNAL inigeom
115      EXTERNAL exner_hyb,addit
116      EXTERNAL defrun_new, test_period
117      REAL  SSUM
118      REAL time_0 , finvmaold(ip1jmp1,llm)
119
120      LOGICAL lafin
121      INTEGER ij,iq,l,numvanle,iapp_tracvl
122
123
124      INTEGER fluxid, fluxvid,fluxdid
125      integer histid, histvid, histaveid
126      real time_step, t_wrt, t_ops
127
128      REAL rdayvrai,rdaym_ini,rday_ecri
129      LOGICAL first
130      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
131
132#ifdef INCA_CH4
133      REAL :: flxw(ip1jmp1,llm)
134#endif
135
136      LOGICAL offline  ! Controle du stockage ds "fluxmass"
137      PARAMETER (offline=.true.)
138
139      character*80 dynhist_file, dynhistave_file
140      character*20 modname
141      character*80 abort_message
142
143C Calendrier
144      LOGICAL true_calendar
145      PARAMETER (true_calendar = .false.)
146C Run nudge
147      LOGICAL ok_nudge
148      PARAMETER (ok_nudge = .false.)
149
150c-----------------------------------------------------------------------
151c   Initialisations:
152c   ----------------
153
154      modname = 'gcm'
155      descript = 'Run GCM LMDZ'
156      lafin    = .FALSE.
157      dynhist_file = 'dyn_hist.nc'
158      dynhistave_file = 'dyn_hist_ave.nc'
159
160      if (true_calendar) then
161        call ioconf_calendar('gregorian')
162      else
163        call ioconf_calendar('360d')
164      endif
165c-----------------------------------------------------------------------
166c
167c  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
168c  ...................................................................
169c
170c     iadv = 1    shema  transport type "humidite specifique LMD" 
171c     iadv = 2    shema   amont
172c     iadv = 3    shema  Van-leer
173c     iadv = 4    schema  Van-leer + humidite specifique
174c                        Modif F.Codron
175c
176c     dans le tableau q(ij,l,iq) , iq = 1  pour l'eau vapeur
177c                                , iq = 2  pour l'eau liquide
178c         et eventuellement      , iq = 3, nqmx pour les autres traceurs
179c
180c     iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
181c
182      DO iq = 1, nqmx
183       iadv( iq ) = 3
184      ENDDO
185c
186       IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
187     * ,' pour la vapeur d''eau'
188       IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'
189     * ,' vapeur d''eau '
190       IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
191     * ,'la vapeur d''eau'
192       IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + humidite'
193     * ,' specifique pour la vapeur d''eau'
194c
195       IF( iadv(1).LE.0.OR.iadv(1).GT.4 )   THEN
196        PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
197     * ,  car  iadv(1) = ', iadv(1)
198         CALL ABORT
199       ENDIF
200
201      DO  iq = 2, nqmx
202       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
203     * ,' pour le traceur no ', iq
204       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'
205     * ,' traceur no ', iq
206       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
207     * ,'le traceur no ', iq
208
209       IF( iadv(iq).EQ.4 )  THEN
210         PRINT *,' Le shema  Van-Leer + humidite specifique ',
211     * ' est  uniquement pour la vapeur d eau .'
212         PRINT *,' Corriger iadv( ',iq, ')  et repasser ! '
213         CALL ABORT
214       ENDIF
215
216       IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 )   THEN
217       PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
218     * . '
219         CALL ABORT
220       ENDIF
221      ENDDO
222c
223      first = .TRUE.
224      numvanle = nqmx + 1
225      DO  iq = 1, nqmx
226        IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
227          numvanle = iq
228          first    = .FALSE.
229        ENDIF
230      ENDDO
231c
232      DO  iq = 1, nqmx
233      IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
234          PRINT *,' Il y a discontinuite dans le choix du shema de ',
235     *    'Van-leer pour les traceurs . Corriger et repasser . '
236           CALL ABORT
237      ENDIF
238        IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 )    THEN
239          PRINT *,' Le choix de  iadv  est errone pour le traceur ',
240     *    iq
241         STOP
242        ENDIF
243      ENDDO
244c
245
246      CALL dynetat0("start.nc",nqmx,vcov,ucov,
247     .              teta,q,masse,ps,phis, time_0)
248
249
250c     OPEN( 99,file ='run.def',status='old',form='formatted')
251      CALL defrun_new( 99, .TRUE. , clesphy0 )
252c     CLOSE(99)
253
254c  on recalcule eventuellement le pas de temps
255
256      IF(MOD(day_step,iperiod).NE.0)
257     * STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
258
259      IF(MOD(day_step,iphysiq).NE.0)
260     * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
261
262      zdtvr    = daysec/FLOAT(day_step)
263        IF(dtvr.NE.zdtvr) THEN
264         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
265        ENDIF
266
267c  nombre d'etats dans les fichiers demarrage et histoire
268      nbetatdem = nday / iecri
269      nbetatmoy = nday / periodav + 1
270
271      dtvr = zdtvr
272      CALL iniconst
273      CALL inigeom
274
275      CALL inifilr
276c
277c   ......    P.Le Van    ( modif  le 29/04/97 )   ......... 
278c
279      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
280     *                tetagdiv, tetagrot , tetatemp              )
281c
282
283      CALL pression ( ip1jmp1, ap, bp, ps, p       )
284      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
285c
286
287c  numero de stockage pour les fichiers de redemarrage:
288
289c-----------------------------------------------------------------------
290c   temps de depart et de fin:
291c   --------------------------
292
293      itau = 0
294      iday = day_ini+itau/day_step
295      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
296         IF(time.GT.1.) THEN
297          time = time-1.
298          iday = iday+1
299         ENDIF
300      itaufin   = nday*day_step
301      itaufinp1 = itaufin +1
302
303      day_end = day_ini + nday
304      PRINT 300, itau,itaufin,day_ini,day_end
305
306      CALL dynredem0("restart.nc",day_end,anne_ini,phis,nqmx)
307
308      ecripar = .TRUE.
309
310      time_step = zdtvr
311      t_ops = iecri * daysec
312      t_wrt = iecri * daysec
313      CALL inithist(dynhist_file,day_ini,anne_ini,time_step,
314     .              t_ops, t_wrt, nqmx, histid, histvid)
315
316      t_ops = iperiod * time_step
317      t_wrt = periodav * daysec
318      CALL initdynav(dynhistave_file,day_ini,anne_ini,time_step,
319     .              t_ops, t_wrt, nqmx, histaveid)
320
321      dtav = iperiod*dtvr/daysec
322
323
324c   Quelques initialisations pour les traceurs
325      call initial0(ijp1llm*nqmx,dq)
326      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
327      istphy=istdyn/iphysiq
328
329c-----------------------------------------------------------------------
330c   Debut de l'integration temporelle:
331c   ----------------------------------
332
333   1  CONTINUE
334
335      if (ok_nudge) then
336        call nudge(itau,ucov,vcov,teta,masse,ps)
337      endif
338c
339c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
340        CALL  test_period ( ucov,vcov,teta,q,p,phis )
341        PRINT *,' ----   Test_period apres continue   OK ! -----', itau
342c     ENDIF
343c
344      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
345      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
346      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
347      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
348      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
349
350      forward = .TRUE.
351      leapf   = .FALSE.
352      dt      =  dtvr
353
354c   ...    P.Le Van .26/04/94  ....
355
356      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
357      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
358
359
360   2  CONTINUE
361
362c-----------------------------------------------------------------------
363
364c   date:
365c   -----
366
367
368c   gestion des appels de la physique et des dissipations:
369c   ------------------------------------------------------
370c
371c   ...    P.Le Van  ( 6/02/95 )  ....
372
373      apphys = .FALSE.
374      statcl = .FALSE.
375      conser = .FALSE.
376      apdiss = .FALSE.
377
378      IF( purmats ) THEN
379         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
380         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
381         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
382     $                              .AND.   physic    ) apphys = .TRUE.
383      ELSE
384         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
385         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
386         IF( MOD(itau+1,iphysiq).EQ.0. AND. physic    ) apphys = .TRUE.
387      END IF
388
389c-----------------------------------------------------------------------
390c   calcul des tendances dynamiques:
391c   --------------------------------
392
393      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
394c
395c
396      CALL caldyn
397     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
398     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
399
400c-----------------------------------------------------------------------
401c   calcul des tendances advection des traceurs (dont l'humidite)
402c   -------------------------------------------------------------
403
404      IF( forward. OR . leapf )  THEN
405c
406        DO 15  iq = 1, nqmx
407c
408         IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
409           CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
410
411         ELSE IF( iq.EQ. nqmx )   THEN
412c
413            iapp_tracvl = 5
414c
415cccc     iapp_tracvl est la frequence en pas du groupement des flux
416cccc      de masse pour  Van-Leer dans la routine  tracvl  .
417c
418      CALL vanleer(numvanle,
419     .             iapp_tracvl,
420     .             nqmx,
421     .             q,
422     .             pbaru,
423     .             pbarv,
424     .             p,
425     .             masse,
426     .             dq, 
427     .             iadv(1),
428     .             teta,
429#ifdef INCA_CH4
430     .             flxw,
431     .             pk,
432     .             mmt_adj,
433     .             adv_flg)
434#else
435     .             pk)
436#endif
437c
438c                   ...  Modif  F.Codron  ....
439c
440         ENDIF
441c
442 15     CONTINUE
443C
444         IF (offline) THEN
445Cmaf stokage du flux de masse pour traceurs OFF-LINE
446
447            CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
448     .   time_step, itau)
449
450
451         ENDIF
452c
453      ENDIF
454
455
456c-----------------------------------------------------------------------
457c   integrations dynamique et traceurs:
458c   ----------------------------------
459
460
461       CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
462     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
463     $              finvmaold                                    )
464
465c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
466c
467c-----------------------------------------------------------------------
468c   calcul des tendances physiques:
469c   -------------------------------
470c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
471c
472       IF( purmats )  THEN
473          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
474       ELSE
475          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
476       ENDIF
477c
478c
479       IF( apphys )  THEN
480c
481c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
482c
483
484         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
485         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
486
487           rdaym_ini  = itau * dtvr / daysec
488           rdayvrai   = rdaym_ini  + day_ini
489
490           IF ( ecritphy.LT.1. )  THEN
491             rday_ecri = rdaym_ini
492           ELSE
493             rday_ecri = NINT( rdayvrai )
494           ENDIF
495c
496        CALL calfis( nqmx,
497     $               lafin,
498     $               rdayvrai,
499     $               rday_ecri,
500     $               time,
501     $               ucov,
502     $               vcov,
503     $               teta,
504     $               q,
505     $               masse,
506     $               ps,
507     $               p,
508     $               pk,
509     $               phis,
510     $               phi,
511     $               du,
512     $               dv,
513     $               dteta,
514     $               dq,
515     $               w,
516#ifdef INCA_CH4
517     $               flxw,
518#endif
519     $               clesphy0,
520     $               dufi,
521     $               dvfi,
522     $               dhfi,
523     $               dqfi,
524     $               dpfi)
525
526c      ajout des tendances physiques:
527c      ------------------------------
528          CALL addfi( nqmx, dtphys, leapf, forward   ,
529     $                  ucov, vcov, teta , q   ,ps ,
530     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
531c
532       ENDIF
533
534        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
535        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
536
537c-----------------------------------------------------------------------
538c
539c
540c   dissipation horizontale et verticale  des petites echelles:
541c   ----------------------------------------------------------
542
543      IF(apdiss) THEN
544
545         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
546
547         CALL addit( ijp1llm,ucov ,dudis,ucov )
548         CALL addit( ijmllm ,vcov ,dvdis,vcov )
549         CALL addit( ijp1llm,teta ,dhdis,teta )
550
551
552c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
553c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
554c
555
556        DO l  =  1, llm
557          DO ij =  1,iim
558           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
559           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
560          ENDDO
561           tpn  = SSUM(iim,tppn,1)/apoln
562           tps  = SSUM(iim,tpps,1)/apols
563
564          DO ij = 1, iip1
565           teta(  ij    ,l) = tpn
566           teta(ij+ip1jm,l) = tps
567          ENDDO
568        ENDDO
569
570        DO ij =  1,iim
571          tppn(ij)  = aire(  ij    ) * ps (  ij    )
572          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
573        ENDDO
574          tpn  = SSUM(iim,tppn,1)/apoln
575          tps  = SSUM(iim,tpps,1)/apols
576
577        DO ij = 1, iip1
578          ps(  ij    ) = tpn
579          ps(ij+ip1jm) = tps
580        ENDDO
581
582
583      END IF
584       
585c   ********************************************************************
586c   ********************************************************************
587c   .... fin de l'integration dynamique  et physique pour le pas itau ..
588c   ********************************************************************
589c   ********************************************************************
590
591c   preparation du pas d'integration suivant  ......
592
593      IF ( .NOT.purmats ) THEN
594c       ........................................................
595c       ..............  schema matsuno + leapfrog  ..............
596c       ........................................................
597
598            IF(forward. OR. leapf) THEN
599              itau= itau + 1
600              iday= day_ini+itau/day_step
601              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
602                IF(time.GT.1.) THEN
603                  time = time-1.
604                  iday = iday+1
605                ENDIF
606            ENDIF
607
608
609            IF( itau. EQ. itaufinp1 ) then 
610              abort_message = 'Simulation finished'
611              call abort_gcm(modname,abort_message,0)
612            ENDIF
613c-----------------------------------------------------------------------
614c   ecriture du fichier histoire moyenne:
615c   -------------------------------------
616
617            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
618               IF(itau.EQ.itaufin) THEN
619                  iav=1
620               ELSE
621                  iav=0
622               ENDIF
623               CALL writedynav(histaveid, nqmx, itau,vcov ,
624     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
625            ENDIF
626
627c-----------------------------------------------------------------------
628c   ecriture de la bande histoire:
629c   ------------------------------
630
631            IF( MOD(itau,iecri*day_step).EQ.0) THEN
632
633               nbetat = nbetatdem
634       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
635       CALL writehist( histid, histvid, nqmx, itau,vcov ,
636     ,                          ucov,teta,phi,q,masse,ps,phis)
637
638
639            ENDIF
640
641            IF(itau.EQ.itaufin) THEN
642
643
644       PRINT *,' Appel test_period avant redem ', itau
645       CALL  test_period ( ucov,vcov,teta,q,p,phis )
646       CALL dynredem1("restart.nc",0.0,
647     .                     vcov,ucov,teta,q,nqmx,masse,ps)
648
649              CLOSE(99)
650            ENDIF
651
652c-----------------------------------------------------------------------
653c   gestion de l'integration temporelle:
654c   ------------------------------------
655
656            IF( MOD(itau,iperiod).EQ.0 )    THEN
657                    GO TO 1
658            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
659
660                   IF( forward )  THEN
661c      fin du pas forward et debut du pas backward
662
663                      forward = .FALSE.
664                        leapf = .FALSE.
665                           GO TO 2
666
667                   ELSE
668c      fin du pas backward et debut du premier pas leapfrog
669
670                        leapf =  .TRUE.
671                        dt  =  2.*dtvr
672                        GO TO 2
673                   END IF
674            ELSE
675
676c      ......   pas leapfrog  .....
677
678                 leapf = .TRUE.
679                 dt  = 2.*dtvr
680                 GO TO 2
681            END IF
682
683      ELSE
684
685c       ........................................................
686c       ..............       schema  matsuno        ...............
687c       ........................................................
688            IF( forward )  THEN
689
690             itau =  itau + 1
691             iday = day_ini+itau/day_step
692             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
693
694                  IF(time.GT.1.) THEN
695                   time = time-1.
696                   iday = iday+1
697                  ENDIF
698
699               forward =  .FALSE.
700               IF( itau. EQ. itaufinp1 ) then 
701                 abort_message = 'Simulation finished'
702                 call abort_gcm(modname,abort_message,0)
703               ENDIF
704               GO TO 2
705
706            ELSE
707
708            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
709               IF(itau.EQ.itaufin) THEN
710                  iav=1
711               ELSE
712                  iav=0
713               ENDIF
714               CALL writedynav(histaveid, nqmx, itau,vcov ,
715     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
716
717            ENDIF
718
719               IF(MOD(itau,iecri*day_step).EQ.0) THEN
720                  nbetat = nbetatdem
721       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
722       CALL writehist( histid, histvid, nqmx, itau,vcov ,
723     ,                          ucov,teta,phi,q,masse,ps,phis)
724               ENDIF
725
726                 IF(itau.EQ.itaufin)
727     . CALL dynredem1("restart.nc",0.0,
728     .                     vcov,ucov,teta,q,nqmx,masse,ps)
729
730                 forward = .TRUE.
731                 GO TO  1
732
733            ENDIF
734
735      END IF
736
737 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
738     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
739
740      STOP
741      END
Note: See TracBrowser for help on using the repository browser.