source: LMDZ5/branches/LF-private/libf/dyn3d/leapfrog.F @ 4400

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

Improved sponge layer:

  • Sponge tendencies are now computed analytically, instead than using a Forward Euler approximation.
  • Sponge tendencies are added within top_bound, and the sponge is applied after physics tendencies have been taken into account.

These changes imply that GCM results (when using sponge layer) will be differentwrt bench test cases using previous revisions.
EM

  • 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 1793 2013-07-18 07:13:18Z dcugnet $
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.