source: LMDZ.3.3/branches/rel-1-0-patch/libf/dyn3d/gcm.F @ 1284

Last change on this file since 1284 was 253, checked in by (none), 23 years ago

This commit was manufactured by cvs2svn to create branch
'rel-1-0-patch'.

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