source: trunk/LMDZ.MARS/libf/dyn3d/gcm.F @ 1036

Last change on this file since 1036 was 1036, checked in by emillour, 11 years ago

Mars GCM: (a first step towards using parallel dynamics)

  • IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
    • advtrac.h (in dyn3d) removed and replaced by module infotrac.F
    • tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
  • Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
  • Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

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