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

Last change on this file since 193 was 190, checked in by lmdzadmin, 24 years ago

Rajout de ok_nudge pour les runs nudges
LF

  • 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
RevLine 
[2]1
2      PROGRAM gcm
3
[57]4      USE IOIPSL
5
[2]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=======================================================================
[78]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.
[2]28
[78]29c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
30c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
[2]31c
[78]32c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
33c
[2]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
[57]121      INTEGER fluxid, fluxvid,fluxdid
[2]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"
[57]130      PARAMETER (offline=.true.)
[2]131
132      character*80 dynhist_file, dynhistave_file
133      character*20 modname
134      character*80 abort_message
135
[57]136C Calendrier
137      LOGICAL true_calendar
138      PARAMETER (true_calendar = .false.)
[190]139C Run nudge
140      LOGICAL ok_nudge
141      PARAMETER (ok_nudge = .false.)
[57]142
[2]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'
[57]152
153      if (true_calendar) then
154        call ioconf_calendar('gregorian')
155      else
156        call ioconf_calendar('360d')
157      endif
[2]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
[78]166c     iadv = 4    schema  Van-leer + humidite specifique
167c                        Modif F.Codron
[2]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
[78]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
[2]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
[78]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
[2]210       PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
211     * . '
[78]212         CALL ABORT
[2]213       ENDIF
214      ENDDO
215c
216      first = .TRUE.
217      numvanle = nqmx + 1
218      DO  iq = 1, nqmx
[78]219        IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
[2]220          numvanle = iq
221          first    = .FALSE.
222        ENDIF
223      ENDDO
224c
225      DO  iq = 1, nqmx
[78]226      IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
[2]227          PRINT *,' Il y a discontinuite dans le choix du shema de ',
228     *    'Van-leer pour les traceurs . Corriger et repasser . '
[78]229           CALL ABORT
230      ENDIF
231        IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 )    THEN
[2]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
[190]327
328      if (ok_nudge) then
329        call nudge(itau,ucov,vcov,teta,masse,ps)
330      endif
[2]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,
[78]412     *                      p, masse, dq,  iadv(1), teta, pk     )
[2]413c
[78]414c                   ...  Modif  F.Codron  ....
415c
[2]416         ENDIF
417c
418 15     CONTINUE
419C
420         IF (offline) THEN
421Cmaf stokage du flux de masse pour traceurs OFF-LINE
[57]422
423            CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
424     .   time_step, itau,fluxid, fluxvid,fluxdid )
425
[2]426         ENDIF
427c
428      ENDIF
429
430
431c-----------------------------------------------------------------------
432c   integrations dynamique et traceurs:
433c   ----------------------------------
434
435
436       CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
437     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
438     $              finvmaold                                    )
439
440c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
441c
442c-----------------------------------------------------------------------
443c   calcul des tendances physiques:
444c   -------------------------------
445c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
446c
447       IF( purmats )  THEN
448          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
449       ELSE
450          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
451       ENDIF
452c
453c
454       IF( apphys )  THEN
455c
456c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
457c
458
459         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
460         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
461
462           rdaym_ini  = itau * dtvr / daysec
463           rdayvrai   = rdaym_ini  + day_ini
464
465           IF ( ecritphy.LT.1. )  THEN
466             rday_ecri = rdaym_ini
467           ELSE
468             rday_ecri = NINT( rdayvrai )
469           ENDIF
470c
471        CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
472     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
473     $     du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi  )
474
475c      ajout des tendances physiques:
476c      ------------------------------
477          CALL addfi( nqmx, dtphys, leapf, forward   ,
478     $                  ucov, vcov, teta , q   ,ps ,
479     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
480c
481       ENDIF
482
483        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
484        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
485
486c-----------------------------------------------------------------------
487c
488c
489c   dissipation horizontale et verticale  des petites echelles:
490c   ----------------------------------------------------------
491
492      IF(apdiss) THEN
493
494         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
495
496         CALL addit( ijp1llm,ucov ,dudis,ucov )
497         CALL addit( ijmllm ,vcov ,dvdis,vcov )
498         CALL addit( ijp1llm,teta ,dhdis,teta )
499
500
501c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
502c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
503c
504
505        DO l  =  1, llm
506          DO ij =  1,iim
507           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
508           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
509          ENDDO
510           tpn  = SSUM(iim,tppn,1)/apoln
511           tps  = SSUM(iim,tpps,1)/apols
512
513          DO ij = 1, iip1
514           teta(  ij    ,l) = tpn
515           teta(ij+ip1jm,l) = tps
516          ENDDO
517        ENDDO
518
[49]519        DO ij =  1,iim
520          tppn(ij)  = aire(  ij    ) * ps (  ij    )
521          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
522        ENDDO
523          tpn  = SSUM(iim,tppn,1)/apoln
524          tps  = SSUM(iim,tpps,1)/apols
[2]525
[49]526        DO ij = 1, iip1
527          ps(  ij    ) = tpn
528          ps(ij+ip1jm) = tps
529        ENDDO
530
531
[2]532      END IF
533       
534c   ********************************************************************
535c   ********************************************************************
536c   .... fin de l'integration dynamique  et physique pour le pas itau ..
537c   ********************************************************************
538c   ********************************************************************
539
540c   preparation du pas d'integration suivant  ......
541
542      IF ( .NOT.purmats ) THEN
543c       ........................................................
544c       ..............  schema matsuno + leapfrog  ..............
545c       ........................................................
546
547            IF(forward. OR. leapf) THEN
548              itau= itau + 1
549              iday= day_ini+itau/day_step
550              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
551                IF(time.GT.1.) THEN
552                  time = time-1.
553                  iday = iday+1
554                ENDIF
555            ENDIF
556
557
558            IF( itau. EQ. itaufinp1 ) then 
559              abort_message = 'Simulation finished'
560              call abort_gcm(modname,abort_message,0)
561            ENDIF
562c-----------------------------------------------------------------------
563c   ecriture du fichier histoire moyenne:
564c   -------------------------------------
565
566            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
567               IF(itau.EQ.itaufin) THEN
568                  iav=1
569               ELSE
570                  iav=0
571               ENDIF
572               CALL writedynav(histaveid, nqmx, itau,vcov ,
573     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
574            ENDIF
575
576c-----------------------------------------------------------------------
577c   ecriture de la bande histoire:
578c   ------------------------------
579
580            IF( MOD(itau,iecri*day_step).EQ.0) THEN
581
582               nbetat = nbetatdem
583       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
584       CALL writehist( histid, histvid, nqmx, itau,vcov ,
585     ,                          ucov,teta,phi,q,masse,ps,phis)
586
587
588            ENDIF
589
590            IF(itau.EQ.itaufin) THEN
591
592
593       PRINT *,' Appel test_period avant redem ', itau
594       CALL  test_period ( ucov,vcov,teta,q,p,phis )
595       CALL dynredem1("restart.nc",0.0,
596     .                     vcov,ucov,teta,q,nqmx,masse,ps)
597
598              CLOSE(99)
599            ENDIF
600
601c-----------------------------------------------------------------------
602c   gestion de l'integration temporelle:
603c   ------------------------------------
604
605            IF( MOD(itau,iperiod).EQ.0 )    THEN
606                    GO TO 1
607            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
608
609                   IF( forward )  THEN
610c      fin du pas forward et debut du pas backward
611
612                      forward = .FALSE.
613                        leapf = .FALSE.
614                           GO TO 2
615
616                   ELSE
617c      fin du pas backward et debut du premier pas leapfrog
618
619                        leapf =  .TRUE.
620                        dt  =  2.*dtvr
621                        GO TO 2
622                   END IF
623            ELSE
624
625c      ......   pas leapfrog  .....
626
627                 leapf = .TRUE.
628                 dt  = 2.*dtvr
629                 GO TO 2
630            END IF
631
632      ELSE
633
634c       ........................................................
635c       ..............       schema  matsuno        ...............
636c       ........................................................
637            IF( forward )  THEN
638
639             itau =  itau + 1
640             iday = day_ini+itau/day_step
641             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
642
643                  IF(time.GT.1.) THEN
644                   time = time-1.
645                   iday = iday+1
646                  ENDIF
647
648               forward =  .FALSE.
649               IF( itau. EQ. itaufinp1 ) then 
650                 abort_message = 'Simulation finished'
651                 call abort_gcm(modname,abort_message,0)
652               ENDIF
653               GO TO 2
654
655            ELSE
656
657            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
658               IF(itau.EQ.itaufin) THEN
659                  iav=1
660               ELSE
661                  iav=0
662               ENDIF
663               CALL writedynav(histaveid, nqmx, itau,vcov ,
664     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
665
666            ENDIF
667
668               IF(MOD(itau,iecri*day_step).EQ.0) THEN
669                  nbetat = nbetatdem
670       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
671       CALL writehist( histid, histvid, nqmx, itau,vcov ,
672     ,                          ucov,teta,phi,q,masse,ps,phis)
673               ENDIF
674
675                 IF(itau.EQ.itaufin)
676     . CALL dynredem1("restart.nc",0.0,
677     .                     vcov,ucov,teta,q,nqmx,masse,ps)
678
679                 forward = .TRUE.
680                 GO TO  1
681
682            ENDIF
683
684      END IF
685
686 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
687     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
688
689      STOP
690      END
Note: See TracBrowser for help on using the repository browser.