source: LMDZ4/branches/V3_test/libf/dyn3d/leapfrog.F @ 4291

Last change on this file since 4291 was 726, checked in by Laurent Fairhead, 18 years ago

Modifications pour rendre INCA plus independant de LMDZ ACo
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.7 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      itaufin   = nday*day_step
188      itaufinp1 = itaufin +1
189
190
191      itau = 0
192      iday = day_ini+itau/day_step
193      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
194         IF(time.GT.1.) THEN
195          time = time-1.
196          iday = iday+1
197         ENDIF
198
199
200c-----------------------------------------------------------------------
201c   On initialise la pression et la fonction d'Exner :
202c   --------------------------------------------------
203
204      dq=0.
205      CALL pression ( ip1jmp1, ap, bp, ps, p       )
206      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
207
208c-----------------------------------------------------------------------
209c   Debut de l'integration temporelle:
210c   ----------------------------------
211
212   1  CONTINUE
213
214
215#ifdef CPP_IOIPSL
216      if (ok_guide.and.(itaufin-itau-1)*dtvr.gt.21600) then
217        call guide(itau,ucov,vcov,teta,q,masse,ps)
218      else
219        IF(prt_level>9)WRITE(*,*)'attention on ne guide pas les ',
220     .    '6 dernieres heures'
221      endif
222#endif
223c
224c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
225c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
226c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
227c     ENDIF
228c
229      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
230      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
231      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
232      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
233      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
234
235      forward = .TRUE.
236      leapf   = .FALSE.
237      dt      =  dtvr
238
239c   ...    P.Le Van .26/04/94  ....
240
241      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
242      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
243
244      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
245
246   2  CONTINUE
247
248c-----------------------------------------------------------------------
249
250c   date:
251c   -----
252
253
254c   gestion des appels de la physique et des dissipations:
255c   ------------------------------------------------------
256c
257c   ...    P.Le Van  ( 6/02/95 )  ....
258
259      apphys = .FALSE.
260      statcl = .FALSE.
261      conser = .FALSE.
262      apdiss = .FALSE.
263
264      IF( purmats ) THEN
265         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
266         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
267         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
268     s          .and. iflag_phys.NE.0                 ) apphys = .TRUE.
269      ELSE
270         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
271         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
272         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.NE.0) apphys=.TRUE.
273      END IF
274
275c-----------------------------------------------------------------------
276c   calcul des tendances dynamiques:
277c   --------------------------------
278
279      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
280
281      CALL caldyn
282     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
283     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
284
285c-----------------------------------------------------------------------
286c   calcul des tendances advection des traceurs (dont l'humidite)
287c   -------------------------------------------------------------
288
289      IF( forward. OR . leapf )  THEN
290
291c
292#ifdef INCA
293             CALL caladvtrac(q,pbaru,pbarv,
294     *                      p, masse, dq,  teta,
295     .             flxw,
296     .             pk)
297#else
298             CALL caladvtrac(q,pbaru,pbarv,
299     *                      p, masse, dq,  teta,
300     .             pk)
301#endif
302
303         IF (offline) THEN
304Cmaf stokage du flux de masse pour traceurs OFF-LINE
305
306#ifdef CPP_IOIPSL
307           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
308     .   dtvr, itau)
309#endif
310
311
312         ENDIF
313c
314      ENDIF
315
316
317c-----------------------------------------------------------------------
318c   integrations dynamique et traceurs:
319c   ----------------------------------
320
321
322       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
323     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
324     $              finvmaold                                    )
325
326
327c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
328c
329c-----------------------------------------------------------------------
330c   calcul des tendances physiques:
331c   -------------------------------
332c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
333c
334       IF( purmats )  THEN
335          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
336       ELSE
337          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
338       ENDIF
339c
340c
341       IF( apphys )  THEN
342c
343c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
344c
345
346         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
347         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
348
349           rdaym_ini  = itau * dtvr / daysec
350           rdayvrai   = rdaym_ini  + day_ini
351
352
353c rajout debug
354c       lafin = .true.
355
356
357c   Inbterface avec les routines de phylmd (phymars ... )
358c   -----------------------------------------------------
359
360#ifdef CPP_PHYS
361c+jld
362
363c  Diagnostique de conservation de l'énergie : initialisation
364      IF (ip_ebil_dyn.ge.1 ) THEN
365          ztit='bil dyn'
366          CALL diagedyn(ztit,2,1,1,dtphys
367     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
368      ENDIF
369c-jld
370cIM : pour sortir les param. du modele dans un fis. netcdf 110106
371      IF (first) THEN
372       first=.false.
373#include "ini_paramLMDZ_dyn.h"
374      ENDIF
375c
376#include "write_paramLMDZ_dyn.h"
377c
378
379        CALL calfis( nq, lafin ,rdayvrai,time  ,
380     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
381     $               du,dv,dteta,dq,w,
382#ifdef INCA
383     $               flxw,
384#endif
385     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
386
387c      ajout des tendances physiques:
388c      ------------------------------
389          CALL addfi( nqmx, dtphys, leapf, forward   ,
390     $                  ucov, vcov, teta , q   ,ps ,
391     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
392c
393c  Diagnostique de conservation de l'énergie : difference
394      IF (ip_ebil_dyn.ge.1 ) THEN
395          ztit='bil phys'
396          CALL diagedyn(ztit,2,1,1,dtphys
397     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
398      ENDIF
399#else
400
401c   Calcul academique de la physique = Rappel Newtonien + fritcion
402c   --------------------------------------------------------------
403       teta(:,:)=teta(:,:)
404     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
405       call friction(ucov,vcov,iphysiq*dtvr)
406
407#endif
408
409c-jld
410       ENDIF
411
412        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
413        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
414
415
416c-----------------------------------------------------------------------
417c   dissipation horizontale et verticale  des petites echelles:
418c   ----------------------------------------------------------
419
420      IF(apdiss) THEN
421
422
423c   calcul de l'energie cinetique avant dissipation
424        call covcont(llm,ucov,vcov,ucont,vcont)
425        call enercin(vcov,ucov,vcont,ucont,ecin0)
426
427c   dissipation
428        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
429        ucov=ucov+dudis
430        vcov=vcov+dvdis
431c       teta=teta+dtetadis
432
433
434c------------------------------------------------------------------------
435        if (dissip_conservative) then
436C       On rajoute la tendance due a la transform. Ec -> E therm. cree
437C       lors de la dissipation
438            call covcont(llm,ucov,vcov,ucont,vcont)
439            call enercin(vcov,ucov,vcont,ucont,ecin)
440            dtetaecdt= (ecin0-ecin)/ pk
441c           teta=teta+dtetaecdt
442            dtetadis=dtetadis+dtetaecdt
443        endif
444        teta=teta+dtetadis
445c------------------------------------------------------------------------
446
447
448c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
449c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
450c
451
452        DO l  =  1, llm
453          DO ij =  1,iim
454           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
455           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
456          ENDDO
457           tpn  = SSUM(iim,tppn,1)/apoln
458           tps  = SSUM(iim,tpps,1)/apols
459
460          DO ij = 1, iip1
461           teta(  ij    ,l) = tpn
462           teta(ij+ip1jm,l) = tps
463          ENDDO
464        ENDDO
465
466        DO ij =  1,iim
467          tppn(ij)  = aire(  ij    ) * ps (  ij    )
468          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
469        ENDDO
470          tpn  = SSUM(iim,tppn,1)/apoln
471          tps  = SSUM(iim,tpps,1)/apols
472
473        DO ij = 1, iip1
474          ps(  ij    ) = tpn
475          ps(ij+ip1jm) = tps
476        ENDDO
477
478
479      END IF
480
481c ajout debug
482c              IF( lafin ) then 
483c                abort_message = 'Simulation finished'
484c                call abort_gcm(modname,abort_message,0)
485c              ENDIF
486       
487c   ********************************************************************
488c   ********************************************************************
489c   .... fin de l'integration dynamique  et physique pour le pas itau ..
490c   ********************************************************************
491c   ********************************************************************
492
493c   preparation du pas d'integration suivant  ......
494
495      IF ( .NOT.purmats ) THEN
496c       ........................................................
497c       ..............  schema matsuno + leapfrog  ..............
498c       ........................................................
499
500            IF(forward. OR. leapf) THEN
501              itau= itau + 1
502              iday= day_ini+itau/day_step
503              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
504                IF(time.GT.1.) THEN
505                  time = time-1.
506                  iday = iday+1
507                ENDIF
508            ENDIF
509
510
511            IF( itau. EQ. itaufinp1 ) then 
512c$$$       write(79,*) 'ucov',ucov
513c$$$       write(80,*) 'vcov',vcov
514c$$$       write(81,*) 'teta',teta
515c$$$       write(82,*) 'ps',ps
516c$$$       write(83,*) 'q',q
517c$$$       WRITE(85,*) 'q1 = ',q(:,:,1)
518c$$$       WRITE(86,*) 'q3 = ',q(:,:,3)
519
520              abort_message = 'Simulation finished'
521
522              call abort_gcm(modname,abort_message,0)
523            ENDIF
524c-----------------------------------------------------------------------
525c   ecriture du fichier histoire moyenne:
526c   -------------------------------------
527
528            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
529               IF(itau.EQ.itaufin) THEN
530                  iav=1
531               ELSE
532                  iav=0
533               ENDIF
534#ifdef CPP_IOIPSL
535              CALL writedynav(histaveid, nqmx, itau,vcov ,
536     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
537               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
538     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
539#endif
540
541            ENDIF
542
543c-----------------------------------------------------------------------
544c   ecriture de la bande histoire:
545c   ------------------------------
546
547            IF( MOD(itau,iecri         ).EQ.0) THEN
548c           IF( MOD(itau,iecri*day_step).EQ.0) THEN
549
550               nbetat = nbetatdem
551       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
552        unat=0.
553        do l=1,llm
554           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
555           vnat(:,l)=vcov(:,l)/cv(:)
556        enddo
557#ifdef CPP_IOIPSL
558c        CALL writehist(histid,histvid, nqmx,itau,vcov,
559c     s                       ucov,teta,phi,q,masse,ps,phis)
560#else
561#include "write_grads_dyn.h"
562#endif
563
564
565            ENDIF
566
567            IF(itau.EQ.itaufin) THEN
568
569
570#ifdef CPP_IOIPSL
571       CALL dynredem1("restart.nc",0.0,
572     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
573#endif
574
575              CLOSE(99)
576            ENDIF
577
578c-----------------------------------------------------------------------
579c   gestion de l'integration temporelle:
580c   ------------------------------------
581
582            IF( MOD(itau,iperiod).EQ.0 )    THEN
583                    GO TO 1
584            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
585
586                   IF( forward )  THEN
587c      fin du pas forward et debut du pas backward
588
589                      forward = .FALSE.
590                        leapf = .FALSE.
591                           GO TO 2
592
593                   ELSE
594c      fin du pas backward et debut du premier pas leapfrog
595
596                        leapf =  .TRUE.
597                        dt  =  2.*dtvr
598                        GO TO 2
599                   END IF
600            ELSE
601
602c      ......   pas leapfrog  .....
603
604                 leapf = .TRUE.
605                 dt  = 2.*dtvr
606                 GO TO 2
607            END IF
608
609      ELSE
610
611c       ........................................................
612c       ..............       schema  matsuno        ...............
613c       ........................................................
614            IF( forward )  THEN
615
616             itau =  itau + 1
617             iday = day_ini+itau/day_step
618             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
619
620                  IF(time.GT.1.) THEN
621                   time = time-1.
622                   iday = iday+1
623                  ENDIF
624
625               forward =  .FALSE.
626               IF( itau. EQ. itaufinp1 ) then 
627                 abort_message = 'Simulation finished'
628                 call abort_gcm(modname,abort_message,0)
629               ENDIF
630               GO TO 2
631
632            ELSE
633
634            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
635               IF(itau.EQ.itaufin) THEN
636                  iav=1
637               ELSE
638                  iav=0
639               ENDIF
640#ifdef CPP_IOIPSL
641              CALL writedynav(histaveid, nqmx, itau,vcov ,
642     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
643               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
644     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
645#endif
646
647            ENDIF
648
649               IF(MOD(itau,iecri         ).EQ.0) THEN
650c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
651                  nbetat = nbetatdem
652       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
653        unat=0.
654        do l=1,llm
655           unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
656           vnat(:,l)=vcov(:,l)/cv(:)
657        enddo
658#ifdef CPP_IOIPSL
659c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
660c    ,                           ucov,teta,phi,q,masse,ps,phis)
661#else
662#include "write_grads_dyn.h"
663#endif
664
665
666               ENDIF
667
668#ifdef CPP_IOIPSL
669                 IF(itau.EQ.itaufin)
670     . CALL dynredem1("restart.nc",0.0,
671     .                     vcov,ucov,teta,q,nqmx,masse,ps)
672#endif
673
674                 forward = .TRUE.
675                 GO TO  1
676
677            ENDIF
678
679      END IF
680
681      STOP
682      END
Note: See TracBrowser for help on using the repository browser.