source: LMDZ5/trunk/libf/dyn3d/leapfrog.F @ 1987

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

Add updating pressure, mass and Exner function (ie: all variables which depend on surface pressure) after adding physics tendencies (which include a surface pressure tendency).
Note that this change induces slight changes in GCM results with respect to previous svn version of the code, even if surface pressure tendency is zero (because of recomputation of polar values as an average over polar points on the dynamics grid).
EM

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