source: LMDZ5/branches/testing/libf/dyn3d/leapfrog.F @ 1795

Last change on this file since 1795 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794


Testing release based on r1794

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 25.4 KB
Line 
1!
2! $Id: leapfrog.F 1795 2013-07-18 08:20:28Z emillour $
3!
4c
5c
6      SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
7     &                    time_0)
8
9
10cIM : pour sortir les param. du modele dans un fis. netcdf 110106
11#ifdef CPP_IOIPSL
12      use IOIPSL
13#endif
14      USE infotrac
15      USE guide_mod, ONLY : guide_main
16      USE write_field
17      USE control_mod
18      IMPLICIT NONE
19
20c      ......   Version  du 10/01/98    ..........
21
22c             avec  coordonnees  verticales hybrides
23c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
24
25c=======================================================================
26c
27c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
28c   -------
29c
30c   Objet:
31c   ------
32c
33c   GCM LMD nouvelle grille
34c
35c=======================================================================
36c
37c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
38c      et possibilite d'appeler une fonction f(y)  a derivee tangente
39c      hyperbolique a la  place de la fonction a derivee sinusoidale.
40
41c  ... Possibilite de choisir le shema pour l'advection de
42c        q  , en modifiant iadv dans traceur.def  (10/02) .
43c
44c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
45c      Pour Van-Leer iadv=10
46c
47c-----------------------------------------------------------------------
48c   Declarations:
49c   -------------
50
51#include "dimensions.h"
52#include "paramet.h"
53#include "comconst.h"
54#include "comdissnew.h"
55#include "comvert.h"
56#include "comgeom.h"
57#include "logic.h"
58#include "temps.h"
59#include "ener.h"
60#include "description.h"
61#include "serre.h"
62!#include "com_io_dyn.h"
63#include "iniprint.h"
64#include "academic.h"
65
66! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
67! #include "clesphys.h"
68
69      INTEGER         longcles
70      PARAMETER     ( longcles = 20 )
71      REAL  clesphy0( longcles )
72
73      real zqmin,zqmax
74
75c   variables dynamiques
76      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
77      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
78      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
79      REAL ps(ip1jmp1)                       ! pression  au sol
80      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
81      REAL pks(ip1jmp1)                      ! exner au  sol
82      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
83      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
84      REAL masse(ip1jmp1,llm)                ! masse d'air
85      REAL phis(ip1jmp1)                     ! geopotentiel au sol
86      REAL phi(ip1jmp1,llm)                  ! geopotentiel
87      REAL w(ip1jmp1,llm)                    ! vitesse verticale
88
89c variables dynamiques intermediaire pour le transport
90      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
91
92c   variables dynamiques au pas -1
93      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
94      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
95      REAL massem1(ip1jmp1,llm)
96
97c   tendances dynamiques
98      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
99      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
100
101c   tendances de la dissipation
102      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
103      REAL dtetadis(ip1jmp1,llm)
104
105c   tendances physiques
106      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
107      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
108
109c   variables pour le fichier histoire
110      REAL dtav      ! intervalle de temps elementaire
111
112      REAL tppn(iim),tpps(iim),tpn,tps
113c
114      INTEGER itau,itaufinp1,iav
115!      INTEGER  iday ! jour julien
116      REAL       time
117
118      REAL  SSUM
119      REAL time_0
120!     REAL finvmaold(ip1jmp1,llm)
121
122cym      LOGICAL  lafin
123      LOGICAL :: lafin=.false.
124      INTEGER ij,iq,l
125      INTEGER ik
126
127      real time_step, t_wrt, t_ops
128
129!      REAL rdayvrai,rdaym_ini
130! jD_cur: jour julien courant
131! jH_cur: heure julienne courante
132      REAL :: jD_cur, jH_cur
133      INTEGER :: an, mois, jour
134      REAL :: secondes
135
136      LOGICAL first,callinigrads
137cIM : pour sortir les param. du modele dans un fis. netcdf 110106
138      save first
139      data first/.true./
140      real dt_cum
141      character*10 infile
142      integer zan, tau0, thoriid
143      integer nid_ctesGCM
144      save nid_ctesGCM
145      real degres
146      real rlong(iip1), rlatg(jjp1)
147      real zx_tmp_2d(iip1,jjp1)
148      integer ndex2d(iip1*jjp1)
149      logical ok_sync
150      parameter (ok_sync = .true.)
151      logical physic
152
153      data callinigrads/.true./
154      character*10 string10
155
156      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
157      REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
158
159c+jld variables test conservation energie
160      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
161C     Tendance de la temp. potentiel d (theta)/ d t due a la
162C     tansformation d'energie cinetique en energie thermique
163C     cree par la dissipation
164      REAL dtetaecdt(ip1jmp1,llm)
165      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
166      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
167      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
168      CHARACTER*15 ztit
169!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
170!IM   SAVE      ip_ebil_dyn
171!IM   DATA      ip_ebil_dyn/0/
172c-jld
173
174      character*80 dynhist_file, dynhistave_file
175      character(len=*),parameter :: modname="leapfrog"
176      character*80 abort_message
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      logical , parameter :: flag_verif = .false.
189     
190
191      integer itau_w   ! pas de temps ecriture = itap + itau_phy
192
193
194      itaufin   = nday*day_step
195      itaufinp1 = itaufin +1
196      itau = 0
197      physic=.true.
198      if (iflag_phys==0.or.iflag_phys==2) physic=.false.
199
200c      iday = day_ini+itau/day_step
201c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
202c         IF(time.GT.1.) THEN
203c          time = time-1.
204c          iday = iday+1
205c         ENDIF
206
207
208c-----------------------------------------------------------------------
209c   On initialise la pression et la fonction d'Exner :
210c   --------------------------------------------------
211
212      dq(:,:,:)=0.
213      CALL pression ( ip1jmp1, ap, bp, ps, p       )
214      if (pressure_exner) then
215        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
216      else
217        CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
218      endif
219
220c-----------------------------------------------------------------------
221c   Debut de l'integration temporelle:
222c   ----------------------------------
223
224   1  CONTINUE ! Matsuno Forward step begins here
225
226      jD_cur = jD_ref + day_ini - day_ref +                             &
227     &          itau/day_step
228      jH_cur = jH_ref + start_time +                                    &
229     &          mod(itau,day_step)/float(day_step)
230      jD_cur = jD_cur + int(jH_cur)
231      jH_cur = jH_cur - int(jH_cur)
232
233
234#ifdef CPP_IOIPSL
235      if (ok_guide) then
236        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
237      endif
238#endif
239
240
241c
242c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
243c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
244c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
245c     ENDIF
246c
247
248! Save fields obtained at previous time step as '...m1'
249      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
250      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
251      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
252      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
253      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
254
255      forward = .TRUE.
256      leapf   = .FALSE.
257      dt      =  dtvr
258
259c   ...    P.Le Van .26/04/94  ....
260! Ehouarn: finvmaold is actually not used
261!      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
262!      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
263
264   2  CONTINUE ! Matsuno backward or leapfrog step begins here
265
266c-----------------------------------------------------------------------
267
268c   date:
269c   -----
270
271
272c   gestion des appels de la physique et des dissipations:
273c   ------------------------------------------------------
274c
275c   ...    P.Le Van  ( 6/02/95 )  ....
276
277      apphys = .FALSE.
278      statcl = .FALSE.
279      conser = .FALSE.
280      apdiss = .FALSE.
281
282      IF( purmats ) THEN
283      ! Purely Matsuno time stepping
284         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
285         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
286     s        apdiss = .TRUE.
287         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
288     s          .and. physic                        ) apphys = .TRUE.
289      ELSE
290      ! Leapfrog/Matsuno time stepping
291         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
292         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
293     s        apdiss = .TRUE.
294         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
295      END IF
296
297! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
298!          supress dissipation step
299      if (llm.eq.1) then
300        apdiss=.false.
301      endif
302
303c-----------------------------------------------------------------------
304c   calcul des tendances dynamiques:
305c   --------------------------------
306
307      ! compute geopotential phi()
308      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
309
310      time = jD_cur + jH_cur
311      CALL caldyn
312     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
313     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
314
315
316c-----------------------------------------------------------------------
317c   calcul des tendances advection des traceurs (dont l'humidite)
318c   -------------------------------------------------------------
319
320      IF( forward. OR . leapf )  THEN
321! Ehouarn: NB: at this point p with ps are not synchronized
322!              (whereas mass and ps are...)
323         CALL caladvtrac(q,pbaru,pbarv,
324     *        p, masse, dq,  teta,
325     .        flxw, pk)
326         
327         IF (offline) THEN
328Cmaf stokage du flux de masse pour traceurs OFF-LINE
329
330#ifdef CPP_IOIPSL
331           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
332     .   dtvr, itau)
333#endif
334
335
336         ENDIF ! of IF (offline)
337c
338      ENDIF ! of IF( forward. OR . leapf )
339
340
341c-----------------------------------------------------------------------
342c   integrations dynamique et traceurs:
343c   ----------------------------------
344
345
346       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
347     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
348!     $              finvmaold                                    )
349
350
351c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
352c
353c-----------------------------------------------------------------------
354c   calcul des tendances physiques:
355c   -------------------------------
356c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
357c
358       IF( purmats )  THEN
359          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
360       ELSE
361          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
362       ENDIF
363c
364c
365       IF( apphys )  THEN
366c
367c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
368c
369
370         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
371         if (pressure_exner) then
372           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
373         else
374           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
375         endif
376
377!           rdaym_ini  = itau * dtvr / daysec
378!           rdayvrai   = rdaym_ini  + day_ini
379!           jD_cur = jD_ref + day_ini - day_ref
380!     $        + int (itau * dtvr / daysec)
381!           jH_cur = jH_ref +                                            &
382!     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
383           jD_cur = jD_ref + day_ini - day_ref +                        &
384     &          itau/day_step
385
386           IF (planet_type .eq."generic") THEN
387              ! AS: we make jD_cur to be pday
388              jD_cur = int(day_ini + itau/day_step)
389           ENDIF
390
391           jH_cur = jH_ref + start_time +                               &
392     &              mod(itau,day_step)/float(day_step)
393           jD_cur = jD_cur + int(jH_cur)
394           jH_cur = jH_cur - int(jH_cur)
395!         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
396!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
397!         write(lunout,*)'current date = ',an, mois, jour, secondes
398
399c rajout debug
400c       lafin = .true.
401
402
403c   Inbterface avec les routines de phylmd (phymars ... )
404c   -----------------------------------------------------
405
406c+jld
407
408c  Diagnostique de conservation de l'énergie : initialisation
409         IF (ip_ebil_dyn.ge.1 ) THEN
410          ztit='bil dyn'
411! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
412           IF (planet_type.eq."earth") THEN
413#ifdef CPP_EARTH
414            CALL diagedyn(ztit,2,1,1,dtphys
415     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
416#endif
417           ENDIF
418         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
419c-jld
420#ifdef CPP_IOIPSL
421cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
422cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
423c        IF (first) THEN
424c         first=.false.
425c#include "ini_paramLMDZ_dyn.h"
426c        ENDIF
427c
428c#include "write_paramLMDZ_dyn.h"
429c
430#endif
431! #endif of #ifdef CPP_IOIPSL
432         CALL calfis( lafin , jD_cur, jH_cur,
433     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
434     $               du,dv,dteta,dq,
435     $               flxw,
436     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
437
438c      ajout des tendances physiques:
439c      ------------------------------
440          CALL addfi( dtphys, leapf, forward   ,
441     $                  ucov, vcov, teta , q   ,ps ,
442     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
443
444         IF (ok_strato) THEN
445           CALL top_bound( vcov,ucov,teta,masse,dtphys)
446         ENDIF
447       
448c
449c  Diagnostique de conservation de l'énergie : difference
450         IF (ip_ebil_dyn.ge.1 ) THEN
451          ztit='bil phys'
452          IF (planet_type.eq."earth") THEN
453           CALL diagedyn(ztit,2,1,1,dtphys
454     &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
455          ENDIF
456         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
457
458       ENDIF ! of IF( apphys )
459
460      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
461!   Academic case : Simple friction and Newtonan relaxation
462!   -------------------------------------------------------
463        DO l=1,llm   
464          DO ij=1,ip1jmp1
465           teta(ij,l)=teta(ij,l)-dtvr*
466     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
467          ENDDO
468        ENDDO ! of DO l=1,llm
469       
470        if (planet_type.eq."giant") then
471          ! add an intrinsic heat flux at the base of the atmosphere
472          teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
473        endif
474
475        call friction(ucov,vcov,dtvr)
476       
477        ! Sponge layer (if any)
478        IF (ok_strato) THEN
479!          dufi(:,:)=0.
480!          dvfi(:,:)=0.
481!          dtetafi(:,:)=0.
482!          dqfi(:,:,:)=0.
483!          dpfi(:)=0.
484!          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
485           CALL top_bound( vcov,ucov,teta,masse,dtvr)
486!          CALL addfi( dtvr, leapf, forward   ,
487!     $                  ucov, vcov, teta , q   ,ps ,
488!     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
489        ENDIF ! of IF (ok_strato)
490      ENDIF ! of IF (iflag_phys.EQ.2)
491
492
493c-jld
494
495        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
496        if (pressure_exner) then
497          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
498        else
499          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
500        endif
501
502
503c-----------------------------------------------------------------------
504c   dissipation horizontale et verticale  des petites echelles:
505c   ----------------------------------------------------------
506
507      IF(apdiss) THEN
508
509
510c   calcul de l'energie cinetique avant dissipation
511        call covcont(llm,ucov,vcov,ucont,vcont)
512        call enercin(vcov,ucov,vcont,ucont,ecin0)
513
514c   dissipation
515        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
516        ucov=ucov+dudis
517        vcov=vcov+dvdis
518c       teta=teta+dtetadis
519
520
521c------------------------------------------------------------------------
522        if (dissip_conservative) then
523C       On rajoute la tendance due a la transform. Ec -> E therm. cree
524C       lors de la dissipation
525            call covcont(llm,ucov,vcov,ucont,vcont)
526            call enercin(vcov,ucov,vcont,ucont,ecin)
527            dtetaecdt= (ecin0-ecin)/ pk
528c           teta=teta+dtetaecdt
529            dtetadis=dtetadis+dtetaecdt
530        endif
531        teta=teta+dtetadis
532c------------------------------------------------------------------------
533
534
535c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
536c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
537c
538
539        DO l  =  1, llm
540          DO ij =  1,iim
541           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
542           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
543          ENDDO
544           tpn  = SSUM(iim,tppn,1)/apoln
545           tps  = SSUM(iim,tpps,1)/apols
546
547          DO ij = 1, iip1
548           teta(  ij    ,l) = tpn
549           teta(ij+ip1jm,l) = tps
550          ENDDO
551        ENDDO
552
553        if (1 == 0) then
554!!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
555!!!                     2) should probably not be here anyway
556!!! but are kept for those who would want to revert to previous behaviour
557           DO ij =  1,iim
558             tppn(ij)  = aire(  ij    ) * ps (  ij    )
559             tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
560           ENDDO
561             tpn  = SSUM(iim,tppn,1)/apoln
562             tps  = SSUM(iim,tpps,1)/apols
563
564           DO ij = 1, iip1
565             ps(  ij    ) = tpn
566             ps(ij+ip1jm) = tps
567           ENDDO
568        endif ! of if (1 == 0)
569
570      END IF ! of IF(apdiss)
571
572c ajout debug
573c              IF( lafin ) then 
574c                abort_message = 'Simulation finished'
575c                call abort_gcm(modname,abort_message,0)
576c              ENDIF
577       
578c   ********************************************************************
579c   ********************************************************************
580c   .... fin de l'integration dynamique  et physique pour le pas itau ..
581c   ********************************************************************
582c   ********************************************************************
583
584c   preparation du pas d'integration suivant  ......
585
586      IF ( .NOT.purmats ) THEN
587c       ........................................................
588c       ..............  schema matsuno + leapfrog  ..............
589c       ........................................................
590
591            IF(forward. OR. leapf) THEN
592              itau= itau + 1
593c              iday= day_ini+itau/day_step
594c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
595c                IF(time.GT.1.) THEN
596c                  time = time-1.
597c                  iday = iday+1
598c                ENDIF
599            ENDIF
600
601
602            IF( itau. EQ. itaufinp1 ) then 
603              if (flag_verif) then
604                write(79,*) 'ucov',ucov
605                write(80,*) 'vcov',vcov
606                write(81,*) 'teta',teta
607                write(82,*) 'ps',ps
608                write(83,*) 'q',q
609                WRITE(85,*) 'q1 = ',q(:,:,1)
610                WRITE(86,*) 'q3 = ',q(:,:,3)
611              endif
612
613              abort_message = 'Simulation finished'
614
615              call abort_gcm(modname,abort_message,0)
616            ENDIF
617c-----------------------------------------------------------------------
618c   ecriture du fichier histoire moyenne:
619c   -------------------------------------
620
621            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
622               IF(itau.EQ.itaufin) THEN
623                  iav=1
624               ELSE
625                  iav=0
626               ENDIF
627               
628               IF (ok_dynzon) THEN
629#ifdef CPP_IOIPSL
630                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
631     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
632#endif
633               END IF
634               IF (ok_dyn_ave) THEN
635#ifdef CPP_IOIPSL
636                 CALL writedynav(itau,vcov,
637     &                 ucov,teta,pk,phi,q,masse,ps,phis)
638#endif
639               ENDIF
640
641            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
642
643c-----------------------------------------------------------------------
644c   ecriture de la bande histoire:
645c   ------------------------------
646
647            IF( MOD(itau,iecri).EQ.0) THEN
648             ! Ehouarn: output only during LF or Backward Matsuno
649             if (leapf.or.(.not.leapf.and.(.not.forward))) then
650              CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
651              unat=0.
652              do l=1,llm
653                unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
654                vnat(:,l)=vcov(:,l)/cv(:)
655              enddo
656#ifdef CPP_IOIPSL
657              if (ok_dyn_ins) then
658!               write(lunout,*) "leapfrog: call writehist, itau=",itau
659               CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
660!               call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
661!               call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
662!              call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
663!               call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
664!               call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
665              endif ! of if (ok_dyn_ins)
666#endif
667! For some Grads outputs of fields
668              if (output_grads_dyn) then
669#include "write_grads_dyn.h"
670              endif
671             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
672            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
673
674            IF(itau.EQ.itaufin) THEN
675
676
677!              if (planet_type.eq."earth") then
678! Write an Earth-format restart file
679                CALL dynredem1("restart.nc",start_time,
680     &                         vcov,ucov,teta,q,masse,ps)
681!              endif ! of if (planet_type.eq."earth")
682
683              CLOSE(99)
684              !!! Ehouarn: Why not stop here and now?
685            ENDIF ! of IF (itau.EQ.itaufin)
686
687c-----------------------------------------------------------------------
688c   gestion de l'integration temporelle:
689c   ------------------------------------
690
691            IF( MOD(itau,iperiod).EQ.0 )    THEN
692                    GO TO 1
693            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
694
695                   IF( forward )  THEN
696c      fin du pas forward et debut du pas backward
697
698                      forward = .FALSE.
699                        leapf = .FALSE.
700                           GO TO 2
701
702                   ELSE
703c      fin du pas backward et debut du premier pas leapfrog
704
705                        leapf =  .TRUE.
706                        dt  =  2.*dtvr
707                        GO TO 2
708                   END IF ! of IF (forward)
709            ELSE
710
711c      ......   pas leapfrog  .....
712
713                 leapf = .TRUE.
714                 dt  = 2.*dtvr
715                 GO TO 2
716            END IF ! of IF (MOD(itau,iperiod).EQ.0)
717                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
718
719      ELSE ! of IF (.not.purmats)
720
721c       ........................................................
722c       ..............       schema  matsuno        ...............
723c       ........................................................
724            IF( forward )  THEN
725
726             itau =  itau + 1
727c             iday = day_ini+itau/day_step
728c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
729c
730c                  IF(time.GT.1.) THEN
731c                   time = time-1.
732c                   iday = iday+1
733c                  ENDIF
734
735               forward =  .FALSE.
736               IF( itau. EQ. itaufinp1 ) then 
737                 abort_message = 'Simulation finished'
738                 call abort_gcm(modname,abort_message,0)
739               ENDIF
740               GO TO 2
741
742            ELSE ! of IF(forward) i.e. backward step
743
744              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
745               IF(itau.EQ.itaufin) THEN
746                  iav=1
747               ELSE
748                  iav=0
749               ENDIF
750
751               IF (ok_dynzon) THEN
752#ifdef CPP_IOIPSL
753                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
754     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
755#endif
756               ENDIF
757               IF (ok_dyn_ave) THEN
758#ifdef CPP_IOIPSL
759                 CALL writedynav(itau,vcov,
760     &                 ucov,teta,pk,phi,q,masse,ps,phis)
761#endif
762               ENDIF
763
764              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
765
766              IF(MOD(itau,iecri         ).EQ.0) THEN
767c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
768                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
769                unat=0.
770                do l=1,llm
771                  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
772                  vnat(:,l)=vcov(:,l)/cv(:)
773                enddo
774#ifdef CPP_IOIPSL
775              if (ok_dyn_ins) then
776!                write(lunout,*) "leapfrog: call writehist (b)",
777!     &                        itau,iecri
778                CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
779              endif ! of if (ok_dyn_ins)
780#endif
781! For some Grads outputs
782                if (output_grads_dyn) then
783#include "write_grads_dyn.h"
784                endif
785
786              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
787
788              IF(itau.EQ.itaufin) THEN
789!                if (planet_type.eq."earth") then
790                  CALL dynredem1("restart.nc",start_time,
791     &                           vcov,ucov,teta,q,masse,ps)
792!                endif ! of if (planet_type.eq."earth")
793              ENDIF ! of IF(itau.EQ.itaufin)
794
795              forward = .TRUE.
796              GO TO  1
797
798            ENDIF ! of IF (forward)
799
800      END IF ! of IF(.not.purmats)
801
802      STOP
803      END
Note: See TracBrowser for help on using the repository browser.