source: LMDZ.3.3/branches/rel-LF/libf/dyn3d/gcm.F @ 420

Last change on this file since 420 was 403, checked in by lmdzadmin, 22 years ago

Conservation de l'energie dans la dynamique (JLD) + diagnostics transports (FH)
IM/LF

  • Property svn:eol-style set to native
  • Property svn:executable set to *
  • Property svn:keywords set to Author Date Id Revision
File size: 23.6 KB
Line 
1
2      PROGRAM gcm
3
4      USE IOIPSL
5
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=======================================================================
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.
28
29c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
30c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
31c
32c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
33c
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
121      INTEGER fluxid, fluxvid,fluxdid
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)
128c+jld variables test conservation energie
129      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
130C     Tendance de la temp. potentiel d (theta)/ d t due a la
131C     tansformation d'energie cinetique en energie thermique
132C     cree par la dissipation
133      REAL dhecdt(ip1jmp1,llm)
134      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
135      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
136      CHARACTER*15 ztit
137      INTEGER   ip_ebil  ! PRINT level for energy conserv. diag.
138      SAVE      ip_ebil
139      DATA      ip_ebil/1/
140c-jld
141
142      LOGICAL offline  ! Controle du stockage ds "fluxmass"
143      PARAMETER (offline=.false.)
144
145      character*80 dynhist_file, dynhistave_file
146      character*20 modname
147      character*80 abort_message
148
149C Calendrier
150      LOGICAL true_calendar
151      PARAMETER (true_calendar = .false.)
152C Run nudge
153      LOGICAL ok_nudge
154      PARAMETER (ok_nudge = .false.)
155
156c-----------------------------------------------------------------------
157c   Initialisations:
158c   ----------------
159
160      abort_message = 'last timestep reached'
161      modname = 'gcm'
162      descript = 'Run GCM LMDZ'
163      lafin    = .FALSE.
164      dynhist_file = 'dyn_hist.nc'
165      dynhistave_file = 'dyn_hist_ave.nc'
166
167      if (true_calendar) then
168        call ioconf_calendar('gregorian')
169      else
170        call ioconf_calendar('360d')
171      endif
172c-----------------------------------------------------------------------
173c
174c  ....  Choix  des shemas d'advection pour l'eau et les traceurs  ...
175c  ...................................................................
176c
177c     iadv = 1    shema  transport type "humidite specifique LMD" 
178c     iadv = 2    shema   amont
179c     iadv = 3    shema  Van-leer
180c     iadv = 4    schema  Van-leer + humidite specifique
181c                        Modif F.Codron
182c
183c     dans le tableau q(ij,l,iq) , iq = 1  pour l'eau vapeur
184c                                , iq = 2  pour l'eau liquide
185c         et eventuellement      , iq = 3, nqmx pour les autres traceurs
186c
187c     iadv(1): choix pour l'eau vap. et  iadv(2) : choix pour l'eau liq.
188c
189      DO iq = 1, nqmx
190       iadv( iq ) = 3
191      ENDDO
192c
193       IF( iadv(1).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
194     * ,' pour la vapeur d''eau'
195       IF( iadv(1).EQ.2 ) PRINT *,' Choix du shema  amont',' pour la'
196     * ,' vapeur d''eau '
197       IF( iadv(1).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
198     * ,'la vapeur d''eau'
199       IF( iadv(1).EQ.4 ) PRINT *,' Choix du shema  Van-Leer + humidite'
200     * ,' specifique pour la vapeur d''eau'
201c
202       IF( iadv(1).LE.0.OR.iadv(1).GT.4 )   THEN
203        PRINT *,' Erreur dans le choix de iadv (1).Corriger et repasser
204     * ,  car  iadv(1) = ', iadv(1)
205         CALL ABORT
206       ENDIF
207
208      DO  iq = 2, nqmx
209       IF( iadv(iq).EQ.1 ) PRINT *,' Choix du shema humidite specifique'
210     * ,' pour le traceur no ', iq
211       IF( iadv(iq).EQ.2 ) PRINT *,' Choix du shema  amont',' pour le'
212     * ,' traceur no ', iq
213       IF( iadv(iq).EQ.3 ) PRINT *,' Choix du shema  Van-Leer ',' pour'
214     * ,'le traceur no ', iq
215
216       IF( iadv(iq).EQ.4 )  THEN
217         PRINT *,' Le shema  Van-Leer + humidite specifique ',
218     * ' est  uniquement pour la vapeur d eau .'
219         PRINT *,' Corriger iadv( ',iq, ')  et repasser ! '
220         CALL ABORT
221       ENDIF
222
223       IF( iadv(iq).LE.0.OR.iadv(iq).GT.4 )   THEN
224       PRINT *,' Erreur dans le  choix de  iadv . Corriger et repasser
225     * . '
226         CALL ABORT
227       ENDIF
228      ENDDO
229c
230      first = .TRUE.
231      numvanle = nqmx + 1
232      DO  iq = 1, nqmx
233        IF((iadv(iq).EQ.3.OR.iadv(iq).EQ.4).AND.first ) THEN
234          numvanle = iq
235          first    = .FALSE.
236        ENDIF
237      ENDDO
238c
239      DO  iq = 1, nqmx
240      IF( (iadv(iq).NE.3.AND.iadv(iq).NE.4).AND.iq.GT.numvanle )  THEN
241          PRINT *,' Il y a discontinuite dans le choix du shema de ',
242     *    'Van-leer pour les traceurs . Corriger et repasser . '
243           CALL ABORT
244      ENDIF
245        IF( iadv(iq).LT.1.OR.iadv(iq).GT.4 )    THEN
246          PRINT *,' Le choix de  iadv  est errone pour le traceur ',
247     *    iq
248         STOP
249        ENDIF
250      ENDDO
251c
252
253      CALL dynetat0("start.nc",nqmx,vcov,ucov,
254     .              teta,q,masse,ps,phis, time_0)
255
256
257c     OPEN( 99,file ='run.def',status='old',form='formatted')
258c      CALL defrun_new( 99, .TRUE. , clesphy0 )
259      CALL conf_gcm( 99, .TRUE. , clesphy0 )
260c     CLOSE(99)
261
262c  on recalcule eventuellement le pas de temps
263
264      IF(MOD(day_step,iperiod).NE.0)
265     * STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
266
267      IF(MOD(day_step,iphysiq).NE.0)
268     * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
269
270      zdtvr    = daysec/FLOAT(day_step)
271        IF(dtvr.NE.zdtvr) THEN
272         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
273        ENDIF
274
275c  nombre d'etats dans les fichiers demarrage et histoire
276      nbetatdem = nday / iecri
277      nbetatmoy = nday / periodav + 1
278
279      dtvr = zdtvr
280      CALL iniconst
281      CALL inigeom
282
283      CALL inifilr
284c
285c   ......    P.Le Van    ( modif  le 29/04/97 )   ......... 
286c
287      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
288     *                tetagdiv, tetagrot , tetatemp              )
289c+jld
290C     initialisation constantes thermo utilisees dans diagedyn
291      IF (ip_ebil.ge.1 ) THEN
292        CALL suphec
293      ENDIF
294c-jld
295c
296
297      CALL pression ( ip1jmp1, ap, bp, ps, p       )
298      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
299c
300
301c  numero de stockage pour les fichiers de redemarrage:
302
303c-----------------------------------------------------------------------
304c   temps de depart et de fin:
305c   --------------------------
306
307      itau = 0
308      iday = day_ini+itau/day_step
309      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
310         IF(time.GT.1.) THEN
311          time = time-1.
312          iday = iday+1
313         ENDIF
314      itaufin   = nday*day_step
315      itaufinp1 = itaufin +1
316
317      day_end = day_ini + nday
318      PRINT 300, itau,itaufin,day_ini,day_end
319
320      CALL dynredem0("restart.nc", day_end, phis, nqmx)
321
322      ecripar = .TRUE.
323
324      time_step = zdtvr
325      t_ops = iecri * daysec
326      t_wrt = iecri * daysec
327C      CALL inithist(dynhist_file,day_ref,annee_ref,time_step,
328c     .              t_ops, t_wrt, nqmx, histid, histvid)
329
330      t_ops = iperiod * time_step
331      t_wrt = periodav * daysec
332      CALL initdynav(dynhistave_file,day_ref,annee_ref,time_step,
333     .              t_ops, t_wrt, nqmx, histaveid)
334
335      dtav = iperiod*dtvr/daysec
336
337
338c   Quelques initialisations pour les traceurs
339      call initial0(ijp1llm*nqmx,dq)
340      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
341      istphy=istdyn/iphysiq
342
343c-----------------------------------------------------------------------
344c   Debut de l'integration temporelle:
345c   ----------------------------------
346
347   1  CONTINUE
348
349
350
351      if (ok_nudge) then
352        call nudge(itau,ucov,vcov,teta,masse,ps)
353      endif
354c
355c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
356        CALL  test_period ( ucov,vcov,teta,q,p,phis )
357        PRINT *,' ----   Test_period apres continue   OK ! -----', itau
358c     ENDIF
359c
360      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
361      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
362      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
363      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
364      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
365
366      forward = .TRUE.
367      leapf   = .FALSE.
368      dt      =  dtvr
369
370c   ...    P.Le Van .26/04/94  ....
371
372      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
373      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
374
375
376   2  CONTINUE
377
378c-----------------------------------------------------------------------
379
380c   date:
381c   -----
382
383
384c   gestion des appels de la physique et des dissipations:
385c   ------------------------------------------------------
386c
387c   ...    P.Le Van  ( 6/02/95 )  ....
388
389      apphys = .FALSE.
390      statcl = .FALSE.
391      conser = .FALSE.
392      apdiss = .FALSE.
393
394      IF( purmats ) THEN
395         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
396         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
397         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
398     $                              .AND.   physic    ) apphys = .TRUE.
399      ELSE
400         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
401         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
402         IF( MOD(itau+1,iphysiq).EQ.0. AND. physic    ) apphys = .TRUE.
403      END IF
404
405c-----------------------------------------------------------------------
406c   calcul des tendances dynamiques:
407c   --------------------------------
408
409      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
410c
411c
412      CALL caldyn
413     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
414     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
415
416c-----------------------------------------------------------------------
417c   calcul des tendances advection des traceurs (dont l'humidite)
418c   -------------------------------------------------------------
419
420      IF( forward. OR . leapf )  THEN
421c
422        DO 15  iq = 1, nqmx
423c
424         IF( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
425           CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
426
427         ELSE IF( iq.EQ. nqmx )   THEN
428c
429            iapp_tracvl = 5
430c
431cccc     iapp_tracvl est la frequence en pas du groupement des flux
432cccc      de masse pour  Van-Leer dans la routine  tracvl  .
433c
434            CALL vanleer(numvanle,iapp_tracvl,nqmx,q,pbaru,pbarv,
435     *                      p, masse, dq,  iadv(1), teta, pk     )
436c
437c                   ...  Modif  F.Codron  ....
438c
439         ENDIF
440c
441 15     CONTINUE
442C
443         IF (offline) THEN
444Cmaf stokage du flux de masse pour traceurs OFF-LINE
445
446            CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
447     .   time_step, itau)
448
449
450         ENDIF
451c
452      ENDIF
453
454
455c-----------------------------------------------------------------------
456c   integrations dynamique et traceurs:
457c   ----------------------------------
458
459
460       CALL integrd ( iadv, 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
461     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
462     $              finvmaold                                    )
463
464c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
465c
466c-----------------------------------------------------------------------
467c   calcul des tendances physiques:
468c   -------------------------------
469c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
470c
471       IF( purmats )  THEN
472          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
473       ELSE
474          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
475       ENDIF
476c
477c
478       IF( apphys )  THEN
479c
480c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
481c
482
483         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
484         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
485
486           rdaym_ini  = itau * dtvr / daysec
487           rdayvrai   = rdaym_ini  + day_ini
488
489           IF ( ecritphy.LT.1. )  THEN
490             rday_ecri = rdaym_ini
491           ELSE
492             rday_ecri = NINT( rdayvrai )
493           ENDIF
494c
495
496c rajout debug
497c       lafin = .true.
498c+jld
499      IF (ip_ebil.ge.1 ) THEN
500          ztit='bil dyn'
501          CALL diagedyn(ztit,2,1,1,dtphys
502     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
503     s  , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
504      ENDIF
505c-jld
506        CALL calfis( nqmx, lafin ,rdayvrai,rday_ecri,time  ,
507     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
508     $     du,dv,dteta,dq,w,clesphy0, dufi,dvfi,dhfi,dqfi,dpfi  )
509
510c      ajout des tendances physiques:
511c      ------------------------------
512          CALL addfi( nqmx, dtphys, leapf, forward   ,
513     $                  ucov, vcov, teta , q   ,ps ,
514     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
515c
516c+jld
517      IF (ip_ebil.ge.1 ) THEN
518          ztit='bil phys'
519          CALL diagedyn(ztit,2,1,1,dtphys
520     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
521     s  , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
522      ENDIF
523c-jld
524       ENDIF
525
526        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
527        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
528
529c-----------------------------------------------------------------------
530c
531c
532c   dissipation horizontale et verticale  des petites echelles:
533c   ----------------------------------------------------------
534
535      IF(apdiss) THEN
536c+jld
537      IF (ip_ebil.ge.2 ) THEN
538          ztit='avant dissip'
539          CALL diagedyn(ztit,2,2,2,dtvr
540     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
541     s  , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
542      ENDIF
543        CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
544        CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin0  )
545c-jld
546
547         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
548
549         CALL addit( ijp1llm,ucov ,dudis,ucov )
550         CALL addit( ijmllm ,vcov ,dvdis,vcov )
551c+jld
552        CALL covcont  ( llm    , ucov    , vcov , ucont, vcont        )
553        CALL enercin ( vcov   , ucov  , vcont     , ucont  , ecin  )
554C
555C       On rajoute la tendance due a la transform. Ec -> E therm. cree
556C       lors de la dissipation
557        DO l  =  1, llm
558          DO    ij     = 1, ip1jmp1
559            dhecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
560            dhdis(ij,l) = dhdis(ij,l) + dhecdt(ij,l)
561          ENDDO
562        ENDDO
563c-jld
564         CALL addit( ijp1llm,teta ,dhdis,teta )
565c+jld
566c
567         IF (ip_ebil.ge.2 ) THEN
568             ztit='apres dissip'
569             CALL diagedyn(ztit,2,2,2,dtdiss
570     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2),aire
571     s  , d_h_vcol , d_qt, d_qw, d_ql, d_ec)
572         ENDIF
573c-jld
574
575c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
576c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
577c
578
579        DO l  =  1, llm
580          DO ij =  1,iim
581           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
582           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
583          ENDDO
584           tpn  = SSUM(iim,tppn,1)/apoln
585           tps  = SSUM(iim,tpps,1)/apols
586
587          DO ij = 1, iip1
588           teta(  ij    ,l) = tpn
589           teta(ij+ip1jm,l) = tps
590          ENDDO
591        ENDDO
592
593        DO ij =  1,iim
594          tppn(ij)  = aire(  ij    ) * ps (  ij    )
595          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
596        ENDDO
597          tpn  = SSUM(iim,tppn,1)/apoln
598          tps  = SSUM(iim,tpps,1)/apols
599
600        DO ij = 1, iip1
601          ps(  ij    ) = tpn
602          ps(ij+ip1jm) = tps
603        ENDDO
604
605
606      END IF
607
608c ajout debug
609c              IF( lafin ) then 
610c                abort_message = 'Simulation finished'
611c                call abort_gcm(modname,abort_message,0)
612c              ENDIF
613       
614c   ********************************************************************
615c   ********************************************************************
616c   .... fin de l'integration dynamique  et physique pour le pas itau ..
617c   ********************************************************************
618c   ********************************************************************
619
620c   preparation du pas d'integration suivant  ......
621
622      IF ( .NOT.purmats ) THEN
623c       ........................................................
624c       ..............  schema matsuno + leapfrog  ..............
625c       ........................................................
626
627            IF(forward. OR. leapf) THEN
628              itau= itau + 1
629              iday= day_ini+itau/day_step
630              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
631                IF(time.GT.1.) THEN
632                  time = time-1.
633                  iday = iday+1
634                ENDIF
635            ENDIF
636
637
638            IF( itau. EQ. itaufinp1 ) then 
639
640              abort_message = 'Simulation finished'
641              call abort_gcm(modname,abort_message,0)
642            ENDIF
643c-----------------------------------------------------------------------
644c   ecriture du fichier histoire moyenne:
645c   -------------------------------------
646
647            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
648               IF(itau.EQ.itaufin) THEN
649                  iav=1
650               ELSE
651                  iav=0
652               ENDIF
653               CALL writedynav(histaveid, nqmx, itau,vcov ,
654     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
655cccIM cf. FH
656               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
657     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
658
659            ENDIF
660
661c-----------------------------------------------------------------------
662c   ecriture de la bande histoire:
663c   ------------------------------
664
665            IF( MOD(itau,iecri*day_step).EQ.0) THEN
666
667               nbetat = nbetatdem
668       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi        )
669c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
670c     ,                          ucov,teta,phi,q,masse,ps,phis)
671
672
673            ENDIF
674
675            IF(itau.EQ.itaufin) THEN
676
677
678       PRINT *,' Appel test_period avant redem ', itau
679       CALL  test_period ( ucov,vcov,teta,q,p,phis )
680       CALL dynredem1("restart.nc",0.0,
681     .                     vcov,ucov,teta,q,nqmx,masse,ps)
682
683              CLOSE(99)
684            ENDIF
685
686c-----------------------------------------------------------------------
687c   gestion de l'integration temporelle:
688c   ------------------------------------
689
690            IF( MOD(itau,iperiod).EQ.0 )    THEN
691                    GO TO 1
692            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
693
694                   IF( forward )  THEN
695c      fin du pas forward et debut du pas backward
696
697                      forward = .FALSE.
698                        leapf = .FALSE.
699                           GO TO 2
700
701                   ELSE
702c      fin du pas backward et debut du premier pas leapfrog
703
704                        leapf =  .TRUE.
705                        dt  =  2.*dtvr
706                        GO TO 2
707                   END IF
708            ELSE
709
710c      ......   pas leapfrog  .....
711
712                 leapf = .TRUE.
713                 dt  = 2.*dtvr
714                 GO TO 2
715            END IF
716
717      ELSE
718
719c       ........................................................
720c       ..............       schema  matsuno        ...............
721c       ........................................................
722            IF( forward )  THEN
723
724             itau =  itau + 1
725             iday = day_ini+itau/day_step
726             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
727
728                  IF(time.GT.1.) THEN
729                   time = time-1.
730                   iday = iday+1
731                  ENDIF
732
733               forward =  .FALSE.
734               IF( itau. EQ. itaufinp1 ) then 
735                 abort_message = 'Simulation finished'
736                 call abort_gcm(modname,abort_message,0)
737               ENDIF
738               GO TO 2
739
740            ELSE
741
742            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
743               IF(itau.EQ.itaufin) THEN
744                  iav=1
745               ELSE
746                  iav=0
747               ENDIF
748               CALL writedynav(histaveid, nqmx, itau,vcov ,
749     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
750cccIM cf. FH
751               call bilan_dyn (2,dtvr*iperiod,dtvr*day_step*periodav,
752     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
753
754            ENDIF
755
756               IF(MOD(itau,iecri*day_step).EQ.0) THEN
757                  nbetat = nbetatdem
758       CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi       )
759c       CALL writehist( histid, histvid, nqmx, itau,vcov ,
760c     ,                          ucov,teta,phi,q,masse,ps,phis)
761               ENDIF
762
763                 IF(itau.EQ.itaufin)
764     . CALL dynredem1("restart.nc",0.0,
765     .                     vcov,ucov,teta,q,nqmx,masse,ps)
766
767                 forward = .TRUE.
768                 GO TO  1
769
770            ENDIF
771
772      END IF
773
774 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
775     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
776
777      STOP
778      END
Note: See TracBrowser for help on using the repository browser.