source: trunk/LMDZ.GENERIC/libf/dyn3d/gcm.F @ 2236

Last change on this file since 2236 was 1594, checked in by emillour, 8 years ago

Mars and Generic "old" dynamics: to further harmonize with LMDZ.COMMON, replace "idissip" with "dissip_period".
EM

File size: 21.6 KB
RevLine 
[135]1      PROGRAM gcm
2
[1576]3      use ioipsl_getincom, only: getin
[1216]4      use infotrac, only: iniadvtrac, nqtot, iadv
[1006]5      use sponge_mod,only: callsponge,mode_sponge,sponge
[1216]6      use control_mod, only: nday, day_step, iperiod, iphysiq,
[1594]7     &                       iconser, dissip_period
[1395]8!      use comgeomphy, only: initcomgeomphy
[1543]9      USE mod_const_mpi, ONLY: COMM_LMDZ
[1403]10      use filtreg_mod, only: inifilr
[1422]11      USE comvert_mod, ONLY: ap,bp
12      USE comconst_mod, ONLY: daysec,dtvr,dtphys,dtdiss,rad,g,r,cpp
13      USE logic_mod, ONLY: ecripar,forward,leapf,apphys,statcl,conser,
[1593]14     .          apdiss,purmats,iflag_phys,apphys
[1422]15      USE temps_mod, ONLY: day_ini,day_end,itaufin,dt
[1523]16      USE iniphysiq_mod, ONLY: iniphysiq
[135]17      IMPLICIT NONE
18
19c      ......   Version  du 10/01/98    ..........
20
21c             avec  coordonnees  verticales hybrides
22c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
23
24c=======================================================================
25c
26c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
27c   -------
28c
29c   Objet:
30c   ------
31c
32c   GCM LMD nouvelle grille
33c
34c=======================================================================
35c
36c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
37c      et possibilite d'appeler une fonction f(y)  a derivee tangente
38c      hyperbolique a la  place de la fonction a derivee sinusoidale.
39
40c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
41c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
42c
43c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
44c
45c-----------------------------------------------------------------------
46c   Declarations:
47c   -------------
48
49#include "dimensions.h"
50#include "paramet.h"
51#include "comdissnew.h"
52#include "comgeom.h"
[1216]53!#include "control.h"
[135]54#include "netcdf.inc"
55#include "tracstoke.h"
[1216]56!#include"advtrac.h"
[135]57
58      INTEGER*4  iday ! jour julien
59      REAL       time ! Heure de la journee en fraction d''1 jour
60      REAL zdtvr
61
62c   variables dynamiques
63      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
64      real, dimension(ip1jmp1,llm) :: teta   ! temperature potentielle
[1216]65      REAL,allocatable :: q(:,:,:)           ! champs advectes
[135]66      REAL ps(ip1jmp1)                       ! pression  au sol
67      REAL pext(ip1jmp1)                     ! pression  extensive
68      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
69      REAL pks(ip1jmp1)                      ! exner au  sol
70      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
71      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
72      REAL masse(ip1jmp1,llm)                ! masse d''air
73      REAL phis(ip1jmp1)                     ! geopotentiel au sol
74      REAL phi(ip1jmp1,llm)                  ! geopotentiel
75      REAL w(ip1jmp1,llm)                    ! vitesse verticale
76
77c variables dynamiques intermediaire pour le transport
78      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
79
80c   variables dynamiques au pas -1
81      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
82
83      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
84      REAL massem1(ip1jmp1,llm)
85
86c   tendances dynamiques
87      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
[1216]88      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
89      REAL,ALLOCATABLE :: dq(:,:,:)
[135]90
91c   tendances de la dissipation
92      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
93      REAL dhdis(ip1jmp1,llm)
94
95c   tendances physiques
96      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
[1216]97      REAL dhfi(ip1jmp1,llm),dpfi(ip1jmp1)
98      REAL,ALLOCATABLE :: dqfi(:,:,:)
[135]99
100c   variables pour le fichier histoire
101      REAL dtav      ! intervalle de temps elementaire
102
103      REAL tppn(iim),tpps(iim),tpn,tps
104c
105!      INTEGER iadv(nqmx) ! indice schema de transport pour le traceur iq
106
107      INTEGER itau,itaufinp1,iav
108
109
110      EXTERNAL caldyn, traceur
[1403]111      EXTERNAL dissip,geopot,iniconst
[135]112      EXTERNAL integrd,SCOPY
113      EXTERNAL inigeom
[993]114      EXTERNAL exner_hyb
[135]115      EXTERNAL defrun_new, test_period
116      REAL  SSUM
117      REAL time_0 , finvmaold(ip1jmp1,llm)
118
119      LOGICAL lafin
120      INTEGER ij,iq,l,ierr,numvanle,iapp_tracvl
121
122      REAL rdayvrai,rdaym_ini,rday_ecri
123!      LOGICAL first
124      REAL beta(ip1jmp1,llm)
125
126      LOGICAL offline  ! Controle du stockage ds "fluxmass"
127      PARAMETER (offline=.false.)
128
129      character*20 modname
130      character*80 abort_message
131
[1576]132! flag to run with or without tracer transport (read from run.def)
133      logical tracer
134     
[135]135C Calendrier
136      LOGICAL true_calendar
137      PARAMETER (true_calendar = .false.)
138
139! flag to set/remove calls to groupeun
140      logical callgroupeun
[253]141      parameter (callgroupeun = .false.)
[135]142
143!     added by RDW for tests without dynamics
144      logical calldyn
145      parameter (calldyn = .true.)
146
[253]147!     added by RW for test
148!      real pmean,airetot
[135]149
[253]150!     added by FF for dissipation / energy conservation tests
151      logical dissip_conservative
152      parameter (dissip_conservative = .false.)
153      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
154!     d (theta)/ d t (pot. temperature) due to the creation
155!     of thermal energy by dissipation
156      REAL dtetaecdt(ip1jmp1,llm)
157      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
158
[1216]159c-----------------------------------------------------------------------
160c    variables pour l'initialisation de la physique :
161c    ------------------------------------------------
[1395]162!      INTEGER ngridmx
163!      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
164!      REAL zcufi(ngridmx),zcvfi(ngridmx)
165!      REAL latfi(ngridmx),lonfi(ngridmx)
166!      REAL airefi(ngridmx)
167!      SAVE latfi, lonfi, airefi
168!      INTEGER i,j
[253]169
[135]170c-----------------------------------------------------------------------
171c   Initialisations:
172c   ----------------
173
174      modname = 'gcm'
175      lafin    = .FALSE.
176
177c-----------------------------------------------------------------------
[1216]178      CALL defrun_new( 99, .TRUE. )
[1576]179      WRITE(*,*) ""
180      WRITE(*,*) "Run with or without tracer transport ?"
181      tracer=.true. ! default value
182      call getin("tracer",tracer)
183      WRITE(*,*)" tracer = ",tracer
[135]184
[1216]185! Initialize tracers
186      CALL iniadvtrac(nqtot,numvanle)
187! Allocation de la tableau q : champs advectes   
188      allocate(q(ip1jmp1,llm,nqtot))
189      allocate(dq(ip1jmp1,llm,nqtot))
190      allocate(dqfi(ip1jmp1,llm,nqtot))
191
[1416]192      CALL dynetat0("start.nc",vcov,ucov,
[135]193     .              teta,q,masse,ps,phis, time_0)
194
195c  on recalcule eventuellement le pas de temps
196
197      IF(MOD(day_step,iperiod).NE.0)
198     * STOP'Il faut choisir un nb de pas par jour multiple de iperiod'
199
200      IF(MOD(day_step,iphysiq).NE.0)
201     * STOP'Il faut choisir un nb de pas par jour multiple de iphysiq'
202
203      zdtvr    = daysec/FLOAT(day_step)
204
205        IF(dtvr.NE.zdtvr) THEN
206         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
207        ENDIF
208
209c  nombre d'etats dans les fichiers demarrage et histoire
210
211      dtvr = zdtvr
212      CALL iniconst
213      CALL inigeom
214      CALL inifilr
215
216c
217c   ......    P.Le Van    ( modif  le 29/04/97 )   ......... 
218c
219
220      CALL inidissip ( lstardis, nitergdiv, nitergrot, niterh   ,
221     *                tetagdiv, tetagrot , tetatemp              )
222c
223
[1216]224c-----------------------------------------------------------------------
225c   Initialisation de la physique :
226c   -------------------------------
227
228!      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
[1395]229!         latfi(1)=rlatu(1)
230!         lonfi(1)=0.
231!         zcufi(1) = cu(1)
232!         zcvfi(1) = cv(1)
233!         DO j=2,jjm
234!            DO i=1,iim
235!               latfi((j-2)*iim+1+i)= rlatu(j)
236!               lonfi((j-2)*iim+1+i)= rlonv(i)
237!               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
238!               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
239!            ENDDO
240!         ENDDO
241!         latfi(ngridmx)= rlatu(jjp1)
242!         lonfi(ngridmx)= 0.
243!         zcufi(ngridmx) = cu(ip1jm+1)
244!         zcvfi(ngridmx) = cv(ip1jm-iim)
245!
246!         ! build airefi(), mesh area on physics grid
247!         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
248!         ! Poles are single points on physics grid
249!         airefi(1)=airefi(1)*iim
250!         airefi(ngridmx)=airefi(ngridmx)*iim
[1216]251
252! Initialisation de la physique: pose probleme quand on tourne
253! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
254! Il faut une cle CPP_PHYS
255!#ifdef CPP_PHYS
[1395]256!         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
[1543]257         CALL iniphysiq(iim,jjm,llm,
258     &                  (jjm-1)*iim+2,comm_lmdz,
259     &                  daysec,day_ini,dtphys,
[1523]260     &                  rlatu,rlatv,rlonu,rlonv,
261     &                  aire,cu,cv,rad,g,r,cpp,
[1593]262     &                iflag_phys)
[1216]263!#endif
264!         call_iniphys=.false.
265!      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
266
[135]267      CALL pression ( ip1jmp1, ap, bp, ps, p       )
268
269      call dump2d(iip1,jjp1,ps,'PRESSION SURFACE')
270      CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
271
[253]272
[135]273c  numero de stockage pour les fichiers de redemarrage:
274
275c-----------------------------------------------------------------------
276c   temps de depart et de fin:
277c   --------------------------
278
279      itau = 0
280      iday = day_ini+itau/day_step
281      time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
282
283
284         IF(time.GT.1.) THEN
285          time = time-1.
286          iday = iday+1
287         ENDIF
288      itaufin   = nday*day_step
289c ********************************
290c      itaufin = 120   ! temporaire !!
291c ********************************
292      itaufinp1 = itaufin +1
293
294      day_end = day_ini + nday
295      PRINT 300, itau,itaufin,day_ini,day_end
296 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
297     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
298
[1416]299      CALL dynredem0("restart.nc",day_end,phis)
[135]300
301      ecripar = .TRUE.
302
303      dtav = iperiod*dtvr/daysec
304
305
306c   Quelques initialisations pour les traceurs
[1216]307      call initial0(ijp1llm*nqtot,dq)
[135]308c     istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
309c     istphy=istdyn/iphysiq
310
311      write(*,*) "gcm: callgroupeun set to:",callgroupeun
312c-----------------------------------------------------------------------
313c   Debut de l'integration temporelle:
314c   ----------------------------------
315
316   1  CONTINUE
317c
318      IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
319        CALL test_period ( ucov,vcov,teta,q,p,phis )
320        PRINT *,' ----   Test_period apres continue   OK ! -----', itau
321      ENDIF
322
323      if (callgroupeun) then
324        call groupeun(jjp1,llm,ucov,.true.)
325        call groupeun(jjm,llm,vcov,.true.)
326        call groupeun(jjp1,llm,teta,.true.)
327        call groupeun(jjp1,llm,masse,.true.)
328        call groupeun(jjp1,1,ps,.false.)
329      endif
330
331      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
332      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
333      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
334      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
335      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
336
337      forward = .TRUE.
338      leapf   = .FALSE.
339      dt      =  dtvr
340
341c   ...    P.Le Van .26/04/94  ....
342
343      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
344      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
345
346
347   2  CONTINUE
348
349c-----------------------------------------------------------------------
350
351c   date:
352c   -----
353
354!      write(*,*) 'GCM: itau=',itau
355
356c   gestion des appels de la physique et des dissipations:
357c   ------------------------------------------------------
358c
359c   ...    P.Le Van  ( 6/02/95 )  ....
360
361      apphys = .FALSE.
362      statcl = .FALSE.
363      conser = .FALSE.
364      apdiss = .FALSE.
365
366      IF( purmats ) THEN
367         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
[1594]368         IF( MOD(itau,dissip_period ).EQ.0
369     &                              .AND..NOT.forward ) apdiss = .TRUE.
[135]370         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
[1593]371     $                        .AND. (iflag_phys.eq.1) ) apphys = .TRUE.
[135]372      ELSE
373         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
[1594]374         IF( MOD(itau+1,dissip_period)  .EQ. 0        ) apdiss = .TRUE.
[1593]375         IF( MOD(itau+1,iphysiq).EQ.0
376     &                     .AND. (iflag_phys.eq.1)    ) apphys = .TRUE.
[135]377      END IF
378
379c-----------------------------------------------------------------------
380c   calcul des tendances dynamiques:
381c   --------------------------------
382
383      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
384
385      if(calldyn)then
386         CALL caldyn
387     $        ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
388     $        phi,conser,du,dv,dteta,dp,w,pbaru,pbarv,time+iday-day_ini)
389      endif
390
391c-----------------------------------------------------------------------
392c   calcul des tendances advection des traceurs (dont l'humidite)
393c   -------------------------------------------------------------
394
[253]395
396
[135]397      if (tracer) then
398       IF( forward. OR . leapf )  THEN
399
[1216]400        DO iq = 1, nqtot
[135]401c
402         IF ( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
403            CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
404
[1216]405         ELSE IF( iq.EQ. nqtot )   THEN
[135]406c
407            iapp_tracvl = 5
408c
409cccc     iapp_tracvl est la frequence en pas du groupement des flux
410cccc      de masse pour  Van-Leer dans la routine  tracvl  .
411c
412
[1216]413            CALL vanleer(numvanle,iapp_tracvl,nqtot,q,pbaru,pbarv,
[135]414     *                      p, masse, dq,  iadv(1), teta, pk     )
415
416c
417c                   ...  Modif  F.Codron  ....
418c
419         ENDIF
420c
421        ENDDO
422C
423c        IF (offline) THEN
424C maf stokage du flux de masse pour traceurs OFF-LINE
425
426c           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
427c    .   time_step, itau)
428
429c        ENDIF
430c
431      ENDIF
432          END IF   ! tracer
433
434
435c-----------------------------------------------------------------------
436c   integrations dynamique et traceurs:
437c   ----------------------------------
438
439          if(calldyn)then
440             CALL integrd(2,vcovm1,ucovm1,tetam1,psm1,massem1,
441     $            dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis,
442     $            finvmaold)
[253]443          else
444             print*,'Currently no dynamics in this GCM...'
[135]445          endif
446
447
448
449
450
451c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
452c
453c-----------------------------------------------------------------------
454c   calcul des tendances physiques:
455c   -------------------------------
456c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
457c
458       IF( purmats )  THEN
459          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
460       ELSE
461          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
462       ENDIF
463c
464c
465       IF( apphys )  THEN
466c
467c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
468c
469
470         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
471         CALL exner_hyb(  ip1jmp1, ps, p,beta,pks, pk, pkf )
472
473           rdaym_ini  = itau * dtvr / daysec
474           rdayvrai   = rdaym_ini  + day_ini
475
[1593]476! Ehouarn: what was this for ??
477!           IF ( ecritphy.LT.1. )  THEN
478!             rday_ecri = rdaym_ini
479!           ELSE
[649]480             rday_ecri = INT(rdaym_ini)+INT(day_ini)
[1593]481!           ENDIF
[135]482c
483
484
485
[1216]486        CALL calfis( nqtot, lafin ,rdayvrai,rday_ecri,time  ,
[135]487     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
[1576]488     $     du,dv,dteta,dq,w, dufi,dvfi,dhfi,dqfi,dpfi)
[135]489
490
491c      ajout des tendances physiques:
492c      ------------------------------
[253]493!        if(1.eq.2)then
[1216]494          CALL addfi( nqtot, dtphys, leapf, forward   ,
[135]495     $                  ucov, vcov, teta , q   ,ps , masse,
496     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
[253]497!       else
498!          print*,'Currently no physics in this GCM...'
499!       endif
[135]500c
501       ENDIF
502
503       CALL pression ( ip1jmp1, ap, bp, ps, p                  )
504
505       CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
506c   ----------------------------------------------------------
507
508c
509c
510c   dissipation horizontale et verticale  des petites echelles:
511c   ----------------------------------------------------------
512
513      IF(apdiss) THEN
514
515c        Sponge layer
516c        ~~~~~~~~~~~~
[1006]517
[135]518         IF (callsponge) THEN
[1006]519           pext(1:ip1jmp1)=ps(1:ip1jmp1)*aire(1:ip1jmp1)
520           CALL sponge(ucov,vcov,teta,pext,dtdiss,mode_sponge)
[135]521         ENDIF
522
523c        Dissipation horizontale
524c        ~~~~~~~~~~~~~~~~~~~~~~~
[253]525
526
527         if(dissip_conservative) then
528!     calculate kinetic energy before dissipation
529            call covcont(llm,ucov,vcov,ucont,vcont)
530            call enercin(vcov,ucov,vcont,ucont,ecin0)
531         end if
532
[135]533         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
534
[993]535         ucov(:,:)=ucov(:,:)+dudis(:,:)
536         vcov(:,:)=vcov(:,:)+dvdis(:,:)
[253]537
538
539         if (dissip_conservative) then
540C           On rajoute la tendance due a la transform. Ec -> E therm. cree
541C           lors de la dissipation
542            call covcont(llm,ucov,vcov,ucont,vcont)
543            call enercin(vcov,ucov,vcont,ucont,ecin)
544            dtetaecdt(:,:) = (ecin0(:,:)-ecin(:,:))/ pk(:,:)
545            dhdis(:,:)     = dhdis(:,:) + dtetaecdt(:,:)
546         endif
547
[993]548         teta(:,:)=teta(:,:)+dhdis(:,:)
[135]549
550
551c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
552c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
553c
554
555        DO l  =  1, llm
556          DO ij =  1,iim
557           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
558           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
559          ENDDO
560
561           tpn  = SSUM(iim,tppn,1)/apoln
562           tps  = SSUM(iim,tpps,1)/apols
563
564          DO ij = 1, iip1
565           teta(  ij    ,l) = tpn
566           teta(ij+ip1jm,l) = tps
567          ENDDO
568        ENDDO
569
570        DO ij =  1,iim
571          tppn(ij)  = aire(  ij    ) * ps (  ij    )
572          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
573        ENDDO
574          tpn  = SSUM(iim,tppn,1)/apoln
575
576          tps  = SSUM(iim,tpps,1)/apols
577
578        DO ij = 1, iip1
579          ps(  ij    ) = tpn
580          ps(ij+ip1jm) = tps
581        ENDDO
582
583
584
585
586
587
588
589      END IF
590       
591c   ********************************************************************
592c   ********************************************************************
593c   .... fin de l'integration dynamique  et physique pour le pas itau ..
594c   ********************************************************************
595c   ********************************************************************
596
597c   preparation du pas d'integration suivant  ......
598
599      IF ( .NOT.purmats ) THEN
600c       ........................................................
601c       ..............  schema matsuno + leapfrog  ..............
602c       ........................................................
603
604            IF(forward. OR. leapf) THEN
605              itau= itau + 1
606              iday= day_ini+itau/day_step
607              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
608                IF(time.GT.1.) THEN
609                  time = time-1.
610                  iday = iday+1
611                ENDIF
612            ENDIF
613
614
615            IF( itau. EQ. itaufinp1 ) then 
616              abort_message = 'Simulation finished'
617              call abort_gcm(modname,abort_message,0)
618            ENDIF
619c-----------------------------------------------------------------------
620c   ecriture du fichier histoire moyenne:
621c   -------------------------------------
622
623c           IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
624c              IF(itau.EQ.itaufin) THEN
625c                 iav=1
626c              ELSE
627c                 iav=0
628c              ENDIF
[1216]629c              CALL writedynav(histaveid, nqtot, itau,vcov ,
[135]630c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
631c           ENDIF
632
633c-----------------------------------------------------------------------
634
635
636            IF(itau.EQ.itaufin) THEN
637
638
639       PRINT *,' Appel test_period avant redem ', itau
640       CALL test_period ( ucov,vcov,teta,q,p,phis )
641       CALL dynredem1("restart.nc",0.0,
[1416]642     .                     vcov,ucov,teta,q,masse,ps)
[135]643
644              CLOSE(99)
645            ENDIF
646
647c-----------------------------------------------------------------------
648c   gestion de l'integration temporelle:
649c   ------------------------------------
650
651            IF( MOD(itau,iperiod).EQ.0 )    THEN
652                    GO TO 1
653            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
654
655                   IF( forward )  THEN
656c      fin du pas forward et debut du pas backward
657
658                      forward = .FALSE.
659                        leapf = .FALSE.
660                           GO TO 2
661
662                   ELSE
663c      fin du pas backward et debut du premier pas leapfrog
664
665                        leapf =  .TRUE.
666                        dt  =  2.*dtvr
667                        GO TO 2
668                   END IF
669            ELSE
670
671c      ......   pas leapfrog  .....
672
673                 leapf = .TRUE.
674                 dt  = 2.*dtvr
675                 GO TO 2
676            END IF
677
678      ELSE
679
680c       ........................................................
681c       ..............       schema  matsuno        ...............
682c       ........................................................
683            IF( forward )  THEN
684
685             itau =  itau + 1
686             iday = day_ini+itau/day_step
687             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
688
689                  IF(time.GT.1.) THEN
690                   time = time-1.
691                   iday = iday+1
692                  ENDIF
693
694               forward =  .FALSE.
695               IF( itau. EQ. itaufinp1 ) then 
696                 abort_message = 'Simulation finished'
697                 call abort_gcm(modname,abort_message,0)
698               ENDIF
699               GO TO 2
700
701            ELSE
702
703            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
704               IF(itau.EQ.itaufin) THEN
705                  iav=1
706               ELSE
707                  iav=0
708               ENDIF
[1216]709c              CALL writedynav(histaveid, nqtot, itau,vcov ,
[135]710c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
711
712            ENDIF
713
714
715                 IF(itau.EQ.itaufin)
716     . CALL dynredem1("restart.nc",0.0,
[1216]717     .                     vcov,ucov,teta,q,nqtot,masse,ps)
[135]718
719                 forward = .TRUE.
720                 GO TO  1
721
722
723            ENDIF
724
725      END IF
726
727      STOP
728      END
Note: See TracBrowser for help on using the repository browser.