source: trunk/LMDZ.MARS/libf/dyn3d/gcm.F @ 1242

Last change on this file since 1242 was 1130, checked in by emillour, 11 years ago

Mars GCM:
Series of changes to enable running in parallel (using LMDZ.COMMON dynamics);
Current LMDZ.MARS can still notheless be compiled and run in serial mode
"as previously".
Summary of main changes:

  • Main programs (newstart, start2archive, xvik) that used to be in dyn3d have been moved to phymars.
  • dyn3d/control.h is now module control_mod.F90
  • rearanged input/outputs routines everywhere to handle serial/MPI cases. physdem.F => phyredem.F90 , phyetat0.F => phyetat0.F90 ; all read/write routines for startfi files are gathered in module iostart.F90
  • added parallelism related routines init_phys_lmdz.F90, comgeomphy.F90, dimphy.F90, iniphysiq.F90, mod_grid_phy_lmdz.F90, mod_phys_lmdz_mpi_data.F90, mod_phys_lmdz_mpi_transfert.F90, mod_phys_lmdz_omp_data.F90, mod_phys_lmdz_omp_transfert.F90, mod_phys_lmdz_para.F90, mod_phys_lmdz_transfert_para.F90 in phymars and mod_const_mpi.F90 in dyn3d (for compliance with parallel case)
  • created generic routines 'planetwide_maxval' and 'planetwide_minval', in module "planetwide_mod", that enable obtaining the min and max of a field over the whole planet.

EM
File size: 22.0 KB
Line 
1      PROGRAM gcm
2
3      use infotrac, only: iniadvtrac, nqtot, iadv
4      use control_mod, only: day_step, iperiod, iphysiq, ndynstep,
5     &                       nday_r, idissip, iconser, ecritstart,
6     &                       ecritphy
7      use comgeomphy, only: initcomgeomphy
8      IMPLICIT NONE
9
10c      ......   Version  du 10/01/98    ..........
11
12c             avec  coordonnees  verticales hybrides
13c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
14
15c=======================================================================
16c
17c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
18c   -------
19c
20c   Objet:
21c   ------
22c
23c   GCM LMD nouvelle grille
24c
25c=======================================================================
26c
27c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
28c      et possibilite d'appeler une fonction f(y)  a derivee tangente
29c      hyperbolique a la  place de la fonction a derivee sinusoidale.
30
31c  ... Possibilite de choisir le shema de Van-leer pour l'advection de
32c        q  , en faisant iadv = 3  dans   traceur  (29/04/97) .
33c
34c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
35c
36c-----------------------------------------------------------------------
37c   Declarations:
38c   -------------
39
40#include "dimensions.h"
41#include "paramet.h"
42#include "comconst.h"
43#include "comdissnew.h"
44#include "comvert.h"
45#include "comgeom.h"
46#include "logic.h"
47#include "temps.h"
48!#include "control.h"
49#include "ener.h"
50#include "netcdf.inc"
51#include "description.h"
52#include "serre.h"
53#include "tracstoke.h"
54#include "sponge.h"
55!#include"advtrac.h"
56
57      INTEGER*4  iday ! jour julien
58      REAL       time ! Heure de la journee en fraction d''1 jour
59      REAL zdtvr
60
61c   variables dynamiques
62      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
63      real, dimension(ip1jmp1,llm) :: teta   ! temperature potentielle
64      REAL,allocatable :: q(:,:,:)           ! champs advectes
65      REAL ps(ip1jmp1)                       ! pression  au sol
66      REAL pext(ip1jmp1)                     ! pression  extensive
67      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
68      REAL pks(ip1jmp1)                      ! exner au  sol
69      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
70      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
71      REAL masse(ip1jmp1,llm)                ! masse d''air
72      REAL phis(ip1jmp1)                     ! geopotentiel au sol
73      REAL phi(ip1jmp1,llm)                  ! geopotentiel
74      REAL w(ip1jmp1,llm)                    ! vitesse verticale
75
76c variables dynamiques intermediaire pour le transport
77      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
78
79c   variables dynamiques au pas -1
80      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
81
82      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
83      REAL massem1(ip1jmp1,llm)
84
85c   tendances dynamiques
86      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
87      REAL dteta(ip1jmp1,llm),dp(ip1jmp1)
88      REAL,ALLOCATABLE :: dq(:,:,:)
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),dpfi(ip1jmp1)
97      REAL,ALLOCATABLE :: dqfi(:,:,:)
98
99c   variables pour le fichier histoire
100      REAL dtav      ! intervalle de temps elementaire
101
102      REAL tppn(iim),tpps(iim),tpn,tps
103c
104
105      INTEGER itau,itaufinp1,iav
106
107
108      EXTERNAL caldyn, traceur
109      EXTERNAL dissip,geopot,iniconst,inifilr
110      EXTERNAL integrd,SCOPY
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,ierr,numvanle,iapp_tracvl
119
120      REAL rdayvrai,rdaym_ini,rday_ecri
121!      LOGICAL first
122      REAL beta(ip1jmp1,llm)
123
124      LOGICAL offline  ! Controle du stockage ds "fluxmass"
125      PARAMETER (offline=.false.)
126
127      character*20 modname
128      character*80 abort_message
129
130      LOGICAL tracer
131          data tracer/.true./
132!      INTEGER nq
133
134C Calendrier
135      LOGICAL true_calendar
136      PARAMETER (true_calendar = .false.)
137
138! flag to set/remove calls to groupeun
139      logical callgroupeun
140      parameter (callgroupeun = .false.)
141
142c-----------------------------------------------------------------------
143c    variables pour l'initialisation de la physique :
144c    ------------------------------------------------
145      INTEGER ngridmx
146      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
147      REAL zcufi(ngridmx),zcvfi(ngridmx)
148      REAL latfi(ngridmx),lonfi(ngridmx)
149      REAL airefi(ngridmx)
150      SAVE latfi, lonfi, airefi
151      INTEGER i,j
152
153c-----------------------------------------------------------------------
154c   Initialisations:
155c   ----------------
156
157      modname = 'gcm'
158      descript = 'Run GCM LMDZ'
159      lafin    = .FALSE.
160
161c-----------------------------------------------------------------------
162      CALL defrun_new( 99, .TRUE. )
163
164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165! FH 2008/05/02
166! A nettoyer. On ne veut qu'une ou deux routines d'interface
167! dynamique -> physique pour l'initialisation
168!#ifdef CPP_PHYS
169      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
170      call initcomgeomphy
171!#endif
172!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173
174! Initialize tracers
175      call iniadvtrac(nqtot,numvanle)
176! Allocation de la tableau q : champs advectes   
177      allocate(q(ip1jmp1,llm,nqtot))
178      allocate(dq(ip1jmp1,llm,nqtot))
179      allocate(dqfi(ip1jmp1,llm,nqtot))
180
181      CALL dynetat0("start.nc",nqtot,vcov,ucov,
182     .              teta,q,masse,ps,phis,time_0)
183
184      ! in case time_0 (because of roundoffs) is close to zero,
185      ! set it to zero to avoid roundoff propagation issues
186      if ((time_0.gt.0.).and.(time_0.lt.(1./day_step))) then
187        write(*,*)"GCM: In start.nc, time=",time_0
188        write(*,*)"     but day_step=",day_step
189        write(*,*)"     and 1./day_step=",1./day_step
190        write(*,*)"     fix this drift by setting time=0"
191        time_0=0.
192      endif
193
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/REAL(day_step)
204        IF(dtvr.NE.zdtvr) THEN
205         PRINT*,'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
206        ENDIF
207
208c  nombre d'etats dans les fichiers demarrage et histoire
209
210      dtvr = zdtvr
211      CALL iniconst
212      CALL inigeom
213
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
224c-----------------------------------------------------------------------
225c   Initialisation de la physique :
226c   -------------------------------
227
228!      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
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
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
256!         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
257         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys,
258     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
259     &                1)
260!     &                iflag_phys)
261!#endif
262!         call_iniphys=.false.
263!      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
264!        call inifis(ngridmx,llm,nqtot,day_ini,daysec,dtphys,
265!     .                latfi,lonfi,airefi,rad,g,r,cpp)
266
267      call dump2d(iip1,jjp1,ps,'PRESSION SURFACE')
268
269c
270c  numero de stockage pour les fichiers de redemarrage:
271
272c-----------------------------------------------------------------------
273c   temps de depart et de fin:
274c   --------------------------
275
276      itau = 0
277      iday = day_ini+itau/day_step
278      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
279         IF(time.GT.1.) THEN
280          time = time-1.
281          iday = iday+1
282         ENDIF
283      if (ndynstep .gt. 0) then
284         itaufin = ndynstep
285      else
286         itaufin=nint(nday_r*day_step) ! nint() to avoid problematic roundoffs
287      endif
288      ! check that this is compatible with call sequence dyn/phys/dissip
289      ! i.e. that itaufin is a multiple of iphysiq and idissip
290      if ((modulo(itaufin,iphysiq).ne.0).or.
291     &    (modulo(itaufin,idissip).ne.0)) then
292        if (ndynstep .gt. 0) then
293       write(*,'(A,I5)')
294     &  "gcm: Problem: incompatibility between ndynstep=",ndynstep
295        else
296       write(*,'((A,F9.2),2(A,I5))')
297     &  "gcm: Problem: incompatibility between nday=",nday_r,
298     &  " day_step=",day_step," which imply itaufin=",itaufin
299        endif
300        write(*,'(2(A,I5))')
301     &   "  whereas iphysiq=",iphysiq," and idissip=",
302     &  idissip
303        stop
304      endif
305!      write(*,*)"gcm: itaufin=",itaufin
306c ********************************
307c      itaufin = 120   ! temporaire !!
308c ********************************
309      itaufinp1 = itaufin +1
310
311      if (ndynstep .gt. 0) then
312        day_end = day_ini
313     &          + floor(float(ndynstep)/float(day_step)+time_0)
314      else
315        day_end = day_ini + floor(nday_r+time_0)
316      endif
317      PRINT 300, itau,itaufin,day_ini,day_end
318 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x, 
319     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
320
321      CALL dynredem0("restart.nc",day_ini,anne_ini,phis,nqtot)
322
323      ecripar = .TRUE.
324
325      dtav = iperiod*dtvr/daysec
326
327
328c   Quelques initialisations pour les traceurs
329      dq(:,:,:)=0
330c     istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
331c     istphy=istdyn/iphysiq
332
333      write(*,*) "gcm: callgroupeun set to:",callgroupeun
334c-----------------------------------------------------------------------
335c   Debut de l'integration temporelle:
336c   ----------------------------------
337
338   1  CONTINUE
339c
340c TN 09/2012. To ensure "1+1=2" in dynamical core :
341c update atmospheric pressure IN the main loop
342      CALL pression ( ip1jmp1, ap, bp, ps, p       )
343      CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
344
345      IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
346        CALL test_period ( ucov,vcov,teta,q,p,phis )
347        write(*,*)' GCM ---- Test_period apres continue   OK ! -----',
348     &            ' itau: ',itau
349      ENDIF
350
351      if (callgroupeun) then
352        call groupeun(jjp1,llm,ucov,.true.)
353        call groupeun(jjm,llm,vcov,.true.)
354        call groupeun(jjp1,llm,teta,.true.)
355        call groupeun(jjp1,llm,masse,.true.)
356        call groupeun(jjp1,1,ps,.false.)
357      endif
358
359      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
360      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
361      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
362      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
363      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
364
365      forward = .TRUE.
366      leapf   = .FALSE.
367      dt      =  dtvr
368
369c   ...    P.Le Van .26/04/94  ....
370
371      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
372      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
373
374
375   2  CONTINUE
376
377c-----------------------------------------------------------------------
378
379c   date:
380c   -----
381
382!      write(*,*) 'GCM: itau=',itau
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
416
417c-----------------------------------------------------------------------
418c   calcul des tendances advection des traceurs (dont l'humidite)
419c   -------------------------------------------------------------
420
421      if (tracer) then
422       IF( forward. OR . leapf )  THEN
423
424        DO iq = 1, nqtot
425c
426         IF ( iadv(iq).EQ.1.OR.iadv(iq).EQ.2 )  THEN
427            CALL traceur( iq,iadv,q,teta,pk,w, pbaru, pbarv, dq )
428
429         ELSE IF( iq.EQ. nqtot )   THEN
430c
431            iapp_tracvl = 5
432c
433cccc     iapp_tracvl est la frequence en pas du groupement des flux
434cccc      de masse pour  Van-Leer dans la routine  tracvl  .
435c
436
437            CALL vanleer(numvanle,iapp_tracvl,nqtot,q,pbaru,pbarv,
438     *                      p, masse, dq,  iadv(1), teta, pk     )
439
440c
441c                   ...  Modif  F.Codron  ....
442c
443         ENDIF
444c
445        ENDDO
446C
447c        IF (offline) THEN
448C maf stokage du flux de masse pour traceurs OFF-LINE
449
450c           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
451c    .   time_step, itau)
452
453c        ENDIF
454c
455      ENDIF
456          END IF   ! tracer
457
458
459c-----------------------------------------------------------------------
460c   integrations dynamique et traceurs:
461c   ----------------------------------
462
463       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
464     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
465     $              finvmaold )
466
467c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
468c
469c-----------------------------------------------------------------------
470c   calcul des tendances physiques:
471c   -------------------------------
472c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
473c
474       IF( purmats )  THEN
475          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
476       ELSE
477          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
478       ENDIF
479c
480c
481       IF( apphys )  THEN
482c
483c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
484c
485
486         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
487         CALL exner_hyb(  ip1jmp1, ps, p,beta,pks, pk, pkf )
488
489           rdaym_ini  = itau * dtvr / daysec + time_0
490           rdayvrai   = rdaym_ini  + day_ini
491
492           IF ( ecritphy.LT.1. )  THEN
493             rday_ecri = rdaym_ini
494           ELSE
495             rday_ecri = INT( rdayvrai )
496           ENDIF
497c
498        CALL calfis( nqtot, lafin ,rdayvrai,rday_ecri,time  ,
499     $                 ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
500     $     du,dv,dteta,dq,w, dufi,dvfi,dhfi,dqfi,dpfi,tracer)
501
502
503c      ajout des tendances physiques:
504c      ------------------------------
505          CALL addfi( nqtot, dtphys, leapf, forward   ,
506     $                  ucov, vcov, teta , q   ,ps , masse,
507     $                 dufi, dvfi, dhfi , dqfi ,dpfi  )
508c
509       ENDIF
510
511       CALL pression ( ip1jmp1, ap, bp, ps, p                  )
512
513       CALL exner_hyb( ip1jmp1, ps, p,beta, pks, pk, pkf )
514c   ----------------------------------------------------------
515
516c
517c
518c   dissipation horizontale et verticale  des petites echelles:
519c   ----------------------------------------------------------
520
521
522      IF(apdiss) THEN
523
524c        Sponge layer
525c        ~~~~~~~~~~~~
526         DO ij=1, ip1jmp1
527            pext(ij)=ps(ij)*aire(ij)
528         ENDDO
529         IF (callsponge) THEN
530            CALL sponge(ucov,vcov,teta,pext,dtdiss,mode_sponge)
531         ENDIF
532
533c        Dissipation horizontale
534c        ~~~~~~~~~~~~~~~~~~~~~~~
535         CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dhdis)
536
537         CALL addit( ijp1llm,ucov ,dudis,ucov )
538         CALL addit( ijmllm ,vcov ,dvdis,vcov )
539         CALL addit( ijp1llm,teta ,dhdis,teta )
540
541
542c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
543c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
544c
545
546        DO l  =  1, llm
547          DO ij =  1,iim
548           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
549           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
550          ENDDO
551
552           tpn  = SSUM(iim,tppn,1)/apoln
553           tps  = SSUM(iim,tpps,1)/apols
554
555          DO ij = 1, iip1
556           teta(  ij    ,l) = tpn
557           teta(ij+ip1jm,l) = tps
558          ENDDO
559        ENDDO
560
561        DO ij =  1,iim
562          tppn(ij)  = aire(  ij    ) * ps (  ij    )
563          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
564        ENDDO
565          tpn  = SSUM(iim,tppn,1)/apoln
566
567          tps  = SSUM(iim,tpps,1)/apols
568
569        DO ij = 1, iip1
570          ps(  ij    ) = tpn
571          ps(ij+ip1jm) = tps
572        ENDDO
573
574
575      END IF
576       
577c   ********************************************************************
578c   ********************************************************************
579c   .... fin de l'integration dynamique  et physique pour le pas itau ..
580c   ********************************************************************
581c   ********************************************************************
582
583c   preparation du pas d'integration suivant  ......
584
585      IF ( .NOT.purmats ) THEN
586c       ........................................................
587c       ..............  schema matsuno + leapfrog  ..............
588c       ........................................................
589
590            IF(forward. OR. leapf) THEN
591              itau= itau + 1
592              iday= day_ini+itau/day_step
593              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
594                !IF(time.GT.1.) THEN
595                IF(time.GE.1.) THEN
596                  time = time-1.
597                  iday = iday+1
598                ENDIF
599            ENDIF
600
601
602            IF( itau. EQ. itaufinp1 ) then 
603              abort_message = 'Simulation finished'
604              call abort_gcm(modname,abort_message,0)
605            ENDIF
606c-----------------------------------------------------------------------
607c   ecriture du fichier histoire moyenne:
608c   -------------------------------------
609
610c           IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
611c              IF(itau.EQ.itaufin) THEN
612c                 iav=1
613c              ELSE
614c                 iav=0
615c              ENDIF
616c              CALL writedynav(histaveid, nqtot, itau,vcov ,
617c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
618c           ENDIF
619
620c-----------------------------------------------------------------------
621
622            IF( (ecritstart.GT.0) .and. (MOD(itau,ecritstart).EQ.0)
623     .        .or. (itau.EQ.itaufin) ) THEN
624           
625       !write(*,*)' GCM: Appel test_period avant redem ; itau=',itau
626       CALL test_period ( ucov,vcov,teta,q,p,phis )
627       
628       write(*,'(A,I7,A,F12.5)')
629     .  'GCM:    Ecriture du fichier restart   ; itau=  ',itau,
630     .  ' date=',REAL(itau)/REAL(day_step)
631       CALL dynredem1("restart.nc",REAL(itau)/REAL(day_step),
632     .                vcov,ucov,teta,q,nqtot,masse,ps)
633     
634      CLOSE(99)
635     
636            ENDIF
637
638
639c-----------------------------------------------------------------------
640c   gestion de l'integration temporelle:
641c   ------------------------------------
642
643            IF( MOD(itau,iperiod).EQ.0 )    THEN
644                    GO TO 1
645            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
646
647                   IF( forward )  THEN
648c      fin du pas forward et debut du pas backward
649
650                      forward = .FALSE.
651                        leapf = .FALSE.
652                           GO TO 2
653
654                   ELSE
655c      fin du pas backward et debut du premier pas leapfrog
656
657                        leapf =  .TRUE.
658                        dt  =  2.*dtvr
659                        GO TO 2
660                   END IF
661            ELSE
662
663c      ......   pas leapfrog  .....
664
665                 leapf = .TRUE.
666                 dt  = 2.*dtvr
667                 GO TO 2
668            END IF
669
670      ELSE
671
672c       ........................................................
673c       ..............       schema  matsuno        ...............
674c       ........................................................
675            IF( forward )  THEN
676
677             itau =  itau + 1
678             iday = day_ini+itau/day_step
679             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
680
681                  IF(time.GE.1.) THEN
682                   time = time-1.
683                   iday = iday+1
684                  ENDIF
685
686               forward =  .FALSE.
687               IF( itau. EQ. itaufinp1 ) then 
688                 abort_message = 'Simulation finished'
689                 call abort_gcm(modname,abort_message,0)
690               ENDIF
691               GO TO 2
692
693            ELSE
694
695            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
696               IF(itau.EQ.itaufin) THEN
697                  iav=1
698               ELSE
699                  iav=0
700               ENDIF
701c              CALL writedynav(histaveid, nqtot, itau,vcov ,
702c    ,                          ucov,teta,pk,phi,q,masse,ps,phis)
703
704            ENDIF
705
706
707            IF( (ecritstart.GT.0) .and. (MOD(itau,ecritstart).EQ.0)
708     .      .or. (itau.EQ.itaufin) ) THEN
709           
710              CALL dynredem1("restart.nc",
711     .                REAL(itau)/REAL(day_step),
712     .                vcov,ucov,teta,q,nqtot,masse,ps)
713     
714              forward = .TRUE.
715              GO TO  1
716             
717            ENDIF
718
719
720            ENDIF
721
722      END IF
723
724      STOP
725      END
Note: See TracBrowser for help on using the repository browser.