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

Last change on this file since 959 was 956, checked in by lmdzadmin, 17 years ago

Nettoyage du controle des parametres physiques. FH

Les parametres cycle_diurne, soil_model, new_oliq, ok_orodr, ok_orolf, ok_limitvrai, nbapp_rad et iflag_con
sont maintenant geres par la physique uniquement.
ecritphy est elimine.
dimphy.F90 et clesphys.h ne sont plus utilises par le code dynamique.
Le test academique obtenu en compilant avec
makegcm -p nophys gcm
fonctionne. FH
IM

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