source: LMDZ4/trunk/libf/dyn3d/leapfrog.F @ 945

Last change on this file since 945 was 825, checked in by Laurent Fairhead, 17 years ago

Probleme de calcul des pas de temps d'ecriture dans le fichier paramLMDZ_dyn.nc
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.8 KB
Line 
1!
2! $Header$
3!
4c
5c
6      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
7     &                    time_0)
8
9
10cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11      USE IOIPSL
12
13      IMPLICIT NONE
14
15c      ......   Version  du 10/01/98    ..........
16
17c             avec  coordonnees  verticales hybrides
18c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
19
20c=======================================================================
21c
22c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
23c   -------
24c
25c   Objet:
26c   ------
27c
28c   GCM LMD nouvelle grille
29c
30c=======================================================================
31c
32c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
33c      et possibilite d'appeler une fonction f(y)  a derivee tangente
34c      hyperbolique a la  place de la fonction a derivee sinusoidale.
35
36c  ... Possibilite de choisir le shema pour l'advection de
37c        q  , en modifiant iadv dans traceur.def  (10/02) .
38c
39c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
40c      Pour Van-Leer iadv=10
41c
42c-----------------------------------------------------------------------
43c   Declarations:
44c   -------------
45
46#include "dimensions.h"
47#include "paramet.h"
48#include "comconst.h"
49#include "comdissnew.h"
50#include "comvert.h"
51#include "comgeom.h"
52#include "logic.h"
53#include "temps.h"
54#include "control.h"
55#include "ener.h"
56#include "description.h"
57#include "serre.h"
58#include "com_io_dyn.h"
59#include "iniprint.h"
60#include "advtrac.h"
61c#include "tracstoke.h"
62
63#include "academic.h"
64
65#include "clesphys.h"
66
67      integer nq
68
69      INTEGER         longcles
70      PARAMETER     ( longcles = 20 )
71      REAL  clesphy0( longcles )
72
73      real zqmin,zqmax
74      INTEGER nbetatmoy, nbetatdem,nbetat
75
76c   variables dynamiques
77      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
78      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
79      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
80      REAL ps(ip1jmp1)                       ! pression  au sol
81      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
82      REAL pks(ip1jmp1)                      ! exner au  sol
83      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
84      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
85      REAL masse(ip1jmp1,llm)                ! masse d'air
86      REAL phis(ip1jmp1)                     ! geopotentiel au sol
87      REAL phi(ip1jmp1,llm)                  ! geopotentiel
88      REAL w(ip1jmp1,llm)                    ! vitesse verticale
89
90c variables dynamiques intermediaire pour le transport
91      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
92
93c   variables dynamiques au pas -1
94      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
95      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
96      REAL massem1(ip1jmp1,llm)
97
98c   tendances dynamiques
99      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
100      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqmx),dp(ip1jmp1)
101
102c   tendances de la dissipation
103      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
104      REAL dtetadis(ip1jmp1,llm)
105
106c   tendances physiques
107      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
108      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqmx),dpfi(ip1jmp1)
109
110c   variables pour le fichier histoire
111      REAL dtav      ! intervalle de temps elementaire
112
113      REAL tppn(iim),tpps(iim),tpn,tps
114c
115      INTEGER itau,itaufinp1,iav
116      INTEGER*4  iday ! jour julien
117      REAL       time ! Heure de la journee en fraction d'1 jour
118
119      REAL  SSUM
120      REAL time_0 , finvmaold(ip1jmp1,llm)
121
122cym      LOGICAL  lafin
123      LOGICAL :: lafin=.false.
124      INTEGER ij,iq,l
125      INTEGER ik
126
127      real time_step, t_wrt, t_ops
128
129      REAL rdayvrai,rdaym_ini
130      LOGICAL first,callinigrads
131cIM : pour sortir les param. du modele dans un fis. netcdf 110106
132      save first
133      data first/.true./
134      real dt_cum, zjulian
135      character*10 infile
136      integer zan, tau0, thoriid
137      integer nid_ctesGCM
138      save nid_ctesGCM
139      real degres
140      real rlong(iip1), rlatg(jjp1)
141      real zx_tmp_2d(iip1,jjp1)
142      integer ndex2d(iip1*jjp1)
143      logical ok_sync
144      parameter (ok_sync = .true.)
145
146      data callinigrads/.true./
147      character*10 string10
148
149      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
150#ifdef INCA
151      REAL :: flxw(ip1jmp1,llm)
152#endif
153
154c+jld variables test conservation energie
155      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
156C     Tendance de la temp. potentiel d (theta)/ d t due a la
157C     tansformation d'energie cinetique en energie thermique
158C     cree par la dissipation
159      REAL dtetaecdt(ip1jmp1,llm)
160      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
161      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
162      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
163      CHARACTER*15 ztit
164!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
165!IM   SAVE      ip_ebil_dyn
166!IM   DATA      ip_ebil_dyn/0/
167c-jld
168
169      character*80 dynhist_file, dynhistave_file
170      character*20 modname
171      character*80 abort_message
172
173C Calendrier
174      LOGICAL true_calendar
175      PARAMETER (true_calendar = .false.)
176
177      logical dissip_conservative
178      save dissip_conservative
179      data dissip_conservative/.true./
180
181      LOGICAL prem
182      save prem
183      DATA prem/.true./
184      INTEGER testita
185      PARAMETER (testita = 9)
186
187      integer itau_w   ! pas de temps ecriture = itap + itau_phy
188
189
190      itaufin   = nday*day_step
191      itaufinp1 = itaufin +1
192
193
194      itau = 0
195      iday = day_ini+itau/day_step
196      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
197         IF(time.GT.1.) THEN
198          time = time-1.
199          iday = iday+1
200         ENDIF
201
202
203c-----------------------------------------------------------------------
204c   On initialise la pression et la fonction d'Exner :
205c   --------------------------------------------------
206
207      dq=0.
208      CALL pression ( ip1jmp1, ap, bp, ps, p       )
209      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
210
211c-----------------------------------------------------------------------
212c   Debut de l'integration temporelle:
213c   ----------------------------------
214
215   1  CONTINUE
216
217
218#ifdef CPP_IOIPSL
219      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
220        call guide(itau,ucov,vcov,teta,q,masse,ps)
221      else
222        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
223     .    '6 dernieres heures'
224      endif
225#endif
226c
227c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
228c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
229c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
230c     ENDIF
231c
232      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
233      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
234      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
235      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
236      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
237
238      forward = .TRUE.
239      leapf   = .FALSE.
240      dt      =  dtvr
241
242c   ...    P.Le Van .26/04/94  ....
243
244      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
245      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
246
247      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
248
249   2  CONTINUE
250
251c-----------------------------------------------------------------------
252
253c   date:
254c   -----
255
256
257c   gestion des appels de la physique et des dissipations:
258c   ------------------------------------------------------
259c
260c   ...    P.Le Van  ( 6/02/95 )  ....
261
262      apphys = .FALSE.
263      statcl = .FALSE.
264      conser = .FALSE.
265      apdiss = .FALSE.
266
267      IF( purmats ) THEN
268         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
269         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
270         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
271     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
272      ELSE
273         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
274         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
275         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
276      END IF
277
278c-----------------------------------------------------------------------
279c   calcul des tendances dynamiques:
280c   --------------------------------
281
282      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
283
284      CALL caldyn
285     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
286     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
287
288c-----------------------------------------------------------------------
289c   calcul des tendances advection des traceurs (dont l'humidite)
290c   -------------------------------------------------------------
291
292      IF( forward. OR . leapf )  THEN
293
294c
295#ifdef INCA
296             CALL caladvtrac(q,pbaru,pbarv,
297     *                      p, masse, dq,  teta,
298     .             flxw,
299     .             pk)
300#else
301             CALL caladvtrac(q,pbaru,pbarv,
302     *                      p, masse, dq,  teta,
303     .             pk)
304#endif
305
306         IF (offline) THEN
307Cmaf stokage du flux de masse pour traceurs OFF-LINE
308
309#ifdef CPP_IOIPSL
310           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
311     .   dtvr, itau)
312#endif
313
314
315         ENDIF
316c
317      ENDIF
318
319
320c-----------------------------------------------------------------------
321c   integrations dynamique et traceurs:
322c   ----------------------------------
323
324
325       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
326     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
327     $              finvmaold                                    )
328
329
330c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
331c
332c-----------------------------------------------------------------------
333c   calcul des tendances physiques:
334c   -------------------------------
335c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
336c
337       IF( purmats )  THEN
338          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
339       ELSE
340          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
341       ENDIF
342c
343c
344       IF( apphys )  THEN
345c
346c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
347c
348
349         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
350         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
351
352           rdaym_ini  = itau * dtvr / daysec
353           rdayvrai   = rdaym_ini  + day_ini
354
355
356c rajout debug
357c       lafin = .true.
358
359
360c   Inbterface avec les routines de phylmd (phymars ... )
361c   -----------------------------------------------------
362
363#ifdef CPP_PHYS
364c+jld
365
366c  Diagnostique de conservation de l'énergie : initialisation
367      IF (ip_ebil_dyn.ge.1 ) THEN
368          ztit='bil dyn'
369          CALL diagedyn(ztit,2,1,1,dtphys
370     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
371      ENDIF
372c-jld
373cIM : pour sortir les param. du modele dans un fis. netcdf 110106
374      IF (first) THEN
375       first=.false.
376#include "ini_paramLMDZ_dyn.h"
377      ENDIF
378c
379#include "write_paramLMDZ_dyn.h"
380c
381
382        CALL calfis( nq, lafin ,rdayvrai,time  ,
383     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
384     $               du,dv,dteta,dq,w,
385#ifdef INCA
386     $               flxw,
387#endif
388     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
389
390c      ajout des tendances physiques:
391c      ------------------------------
392          CALL addfi( nqmx, dtphys, leapf, forward   ,
393     $                  ucov, vcov, teta , q   ,ps ,
394     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
395c
396c  Diagnostique de conservation de l'énergie : difference
397      IF (ip_ebil_dyn.ge.1 ) THEN
398          ztit='bil phys'
399          CALL diagedyn(ztit,2,1,1,dtphys
400     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
401      ENDIF
402#else
403
404c   Calcul academique de la physique = Rappel Newtonien + fritcion
405c   --------------------------------------------------------------
406       teta(:,:)=teta(:,:)
407     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
408       call friction(ucov,vcov,iphysiq*dtvr)
409
410#endif
411
412c-jld
413       ENDIF
414
415        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
416        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
417
418
419c-----------------------------------------------------------------------
420c   dissipation horizontale et verticale  des petites echelles:
421c   ----------------------------------------------------------
422
423      IF(apdiss) THEN
424
425
426c   calcul de l'energie cinetique avant dissipation
427        call covcont(llm,ucov,vcov,ucont,vcont)
428        call enercin(vcov,ucov,vcont,ucont,ecin0)
429
430c   dissipation
431        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
432        ucov=ucov+dudis
433        vcov=vcov+dvdis
434c       teta=teta+dtetadis
435
436
437c------------------------------------------------------------------------
438        if (dissip_conservative) then
439C       On rajoute la tendance due a la transform. Ec -> E therm. cree
440C       lors de la dissipation
441            call covcont(llm,ucov,vcov,ucont,vcont)
442            call enercin(vcov,ucov,vcont,ucont,ecin)
443            dtetaecdt= (ecin0-ecin)/ pk
444c           teta=teta+dtetaecdt
445            dtetadis=dtetadis+dtetaecdt
446        endif
447        teta=teta+dtetadis
448c------------------------------------------------------------------------
449
450
451c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
452c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
453c
454
455        DO l  =  1, llm
456          DO ij =  1,iim
457           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
458           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
459          ENDDO
460           tpn  = SSUM(iim,tppn,1)/apoln
461           tps  = SSUM(iim,tpps,1)/apols
462
463          DO ij = 1, iip1
464           teta(  ij    ,l) = tpn
465           teta(ij+ip1jm,l) = tps
466          ENDDO
467        ENDDO
468
469        DO ij =  1,iim
470          tppn(ij)  = aire(  ij    ) * ps (  ij    )
471          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
472        ENDDO
473          tpn  = SSUM(iim,tppn,1)/apoln
474          tps  = SSUM(iim,tpps,1)/apols
475
476        DO ij = 1, iip1
477          ps(  ij    ) = tpn
478          ps(ij+ip1jm) = tps
479        ENDDO
480
481
482      END IF
483
484c ajout debug
485c              IF( lafin ) then 
486c                abort_message = 'Simulation finished'
487c                call abort_gcm(modname,abort_message,0)
488c              ENDIF
489       
490c   ********************************************************************
491c   ********************************************************************
492c   .... fin de l'integration dynamique  et physique pour le pas itau ..
493c   ********************************************************************
494c   ********************************************************************
495
496c   preparation du pas d'integration suivant  ......
497
498      IF ( .NOT.purmats ) THEN
499c       ........................................................
500c       ..............  schema matsuno + leapfrog  ..............
501c       ........................................................
502
503            IF(forward. OR. leapf) THEN
504              itau= itau + 1
505              iday= day_ini+itau/day_step
506              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
507                IF(time.GT.1.) THEN
508                  time = time-1.
509                  iday = iday+1
510                ENDIF
511            ENDIF
512
513
514            IF( itau. EQ. itaufinp1 ) then 
515c              write(79,*) 'ucov',ucov
516c              write(80,*) 'vcov',vcov
517c              write(81,*) 'teta',teta
518c              write(82,*) 'ps',ps
519c              write(83,*) 'q',q
520c              WRITE(85,*) 'q1 = ',q(:,:,1)
521c              WRITE(86,*) 'q3 = ',q(:,:,3)
522
523              abort_message = 'Simulation finished'
524
525              call abort_gcm(modname,abort_message,0)
526            ENDIF
527c-----------------------------------------------------------------------
528c   ecriture du fichier histoire moyenne:
529c   -------------------------------------
530
531            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
532               IF(itau.EQ.itaufin) THEN
533                  iav=1
534               ELSE
535                  iav=0
536               ENDIF
537#ifdef CPP_IOIPSL
538              CALL writedynav(histaveid, nqmx, itau,vcov ,
539     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
540               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
541     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
542#endif
543
544            ENDIF
545
546c-----------------------------------------------------------------------
547c   ecriture de la bande histoire:
548c   ------------------------------
549
550            IF( MOD(itau,iecri         ).EQ.0) THEN
551c           IF( MOD(itau,iecri*day_step).EQ.0) THEN
552
553               nbetat = nbetatdem
554       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
555        unat=0.
556        do l=1,llm
557           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
558           vnat(:,l)=vcov(:,l)/cv(:)
559        enddo
560#ifdef CPP_IOIPSL
561c        CALL writehist(histid,histvid, nqmx,itau,vcov,
562c     s                       ucov,teta,phi,q,masse,ps,phis)
563#else
564#include "write_grads_dyn.h"
565#endif
566
567
568            ENDIF
569
570            IF(itau.EQ.itaufin) THEN
571
572
573#ifdef CPP_IOIPSL
574       CALL dynredem1("restart.nc",0.0,
575     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
576#endif
577
578              CLOSE(99)
579            ENDIF
580
581c-----------------------------------------------------------------------
582c   gestion de l'integration temporelle:
583c   ------------------------------------
584
585            IF( MOD(itau,iperiod).EQ.0 )    THEN
586                    GO TO 1
587            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
588
589                   IF( forward )  THEN
590c      fin du pas forward et debut du pas backward
591
592                      forward = .FALSE.
593                        leapf = .FALSE.
594                           GO TO 2
595
596                   ELSE
597c      fin du pas backward et debut du premier pas leapfrog
598
599                        leapf =  .TRUE.
600                        dt  =  2.*dtvr
601                        GO TO 2
602                   END IF
603            ELSE
604
605c      ......   pas leapfrog  .....
606
607                 leapf = .TRUE.
608                 dt  = 2.*dtvr
609                 GO TO 2
610            END IF
611
612      ELSE
613
614c       ........................................................
615c       ..............       schema  matsuno        ...............
616c       ........................................................
617            IF( forward )  THEN
618
619             itau =  itau + 1
620             iday = day_ini+itau/day_step
621             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
622
623                  IF(time.GT.1.) THEN
624                   time = time-1.
625                   iday = iday+1
626                  ENDIF
627
628               forward =  .FALSE.
629               IF( itau. EQ. itaufinp1 ) then 
630                 abort_message = 'Simulation finished'
631                 call abort_gcm(modname,abort_message,0)
632               ENDIF
633               GO TO 2
634
635            ELSE
636
637            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
638               IF(itau.EQ.itaufin) THEN
639                  iav=1
640               ELSE
641                  iav=0
642               ENDIF
643#ifdef CPP_IOIPSL
644              CALL writedynav(histaveid, nqmx, itau,vcov ,
645     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
646               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
647     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
648#endif
649
650            ENDIF
651
652               IF(MOD(itau,iecri         ).EQ.0) THEN
653c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
654                  nbetat = nbetatdem
655       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
656        unat=0.
657        do l=1,llm
658           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
659           vnat(:,l)=vcov(:,l)/cv(:)
660        enddo
661#ifdef CPP_IOIPSL
662c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
663c    ,                           ucov,teta,phi,q,masse,ps,phis)
664#else
665#include "write_grads_dyn.h"
666#endif
667
668
669               ENDIF
670
671#ifdef CPP_IOIPSL
672                 IF(itau.EQ.itaufin)
673     . CALL dynredem1("restart.nc",0.0,
674     .                     vcov,ucov,teta,q,nqmx,masse,ps)
675#endif
676
677                 forward = .TRUE.
678                 GO TO  1
679
680            ENDIF
681
682      END IF
683
684      STOP
685      END
Note: See TracBrowser for help on using the repository browser.