source: trunk/libf/dyn3d/leapfrog.F @ 91

Last change on this file since 91 was 53, checked in by aslmd, 14 years ago

modele principal: commit mineur:

M libf/dyn3d/leapfrog.F
M libf/dyn3d/comconst.h
M libf/dyn3d/iniacademic.F
M libf/dyn3d/conf_planete.F90
propagation des changements du commit 52 a la dynamique sequentielle

File size: 26.4 KB
RevLine 
[1]1!
[7]2! $Id: leapfrog.F 1446 2010-10-22 09:27:25Z emillour $
[1]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      INTEGER nbetatmoy, nbetatdem,nbetat
75
76c   variables dynamiques
77      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
78      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
79      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
80      REAL ps(ip1jmp1)                       ! pression  au sol
81      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
82      REAL pks(ip1jmp1)                      ! exner au  sol
83      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
84      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
85      REAL masse(ip1jmp1,llm)                ! masse d'air
86      REAL phis(ip1jmp1)                     ! geopotentiel au sol
87      REAL phi(ip1jmp1,llm)                  ! geopotentiel
88      REAL w(ip1jmp1,llm)                    ! vitesse verticale
[5]89! ADAPTATION GCM POUR CP(T)
90      REAL temp(ip1jmp1,llm)                 ! temperature 
91      REAL tsurpk(ip1jmp1,llm)               ! cpp*T/pk 
[1]92
93c variables dynamiques intermediaire pour le transport
94      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
95
96c   variables dynamiques au pas -1
97      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
98      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
99      REAL massem1(ip1jmp1,llm)
100
[6]101c   tendances dynamiques en */s
[1]102      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
103      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
104
[6]105c   tendances de la dissipation en */s
[1]106      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
107      REAL dtetadis(ip1jmp1,llm)
108
[6]109c   tendances de la couche superieure */s
110      REAL dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
111      REAL dtetatop(ip1jmp1,llm)
112      REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
113
114c   tendances physiques */s
[1]115      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
116      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
117
[37]118      real :: duspg(ip1jmp1,llm) ! for bilan_dyn
119
[1]120c   variables pour le fichier histoire
121      REAL dtav      ! intervalle de temps elementaire
122
123      REAL tppn(iim),tpps(iim),tpn,tps
124c
125      INTEGER itau,itaufinp1,iav
126!      INTEGER  iday ! jour julien
127      REAL       time
128
129      REAL  SSUM
130      REAL time_0 , finvmaold(ip1jmp1,llm)
131
132cym      LOGICAL  lafin
133      LOGICAL :: lafin=.false.
134      INTEGER ij,iq,l
135      INTEGER ik
136
137      real time_step, t_wrt, t_ops
138
139!      REAL rdayvrai,rdaym_ini
140! jD_cur: jour julien courant
141! jH_cur: heure julienne courante
142      REAL :: jD_cur, jH_cur
143      INTEGER :: an, mois, jour
144      REAL :: secondes
145
146      LOGICAL first,callinigrads
147cIM : pour sortir les param. du modele dans un fis. netcdf 110106
148      save first
149      data first/.true./
150      real dt_cum
151      character*10 infile
152      integer zan, tau0, thoriid
153      integer nid_ctesGCM
154      save nid_ctesGCM
155      real degres
156      real rlong(iip1), rlatg(jjp1)
157      real zx_tmp_2d(iip1,jjp1)
158      integer ndex2d(iip1*jjp1)
159      logical ok_sync
160      parameter (ok_sync = .true.)
161
162      data callinigrads/.true./
163      character*10 string10
164
165      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
166      REAL :: flxw(ip1jmp1,llm)  ! flux de masse verticale
167
168c+jld variables test conservation energie
169      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
170C     Tendance de la temp. potentiel d (theta)/ d t due a la
171C     tansformation d'energie cinetique en energie thermique
172C     cree par la dissipation
173      REAL dtetaecdt(ip1jmp1,llm)
174      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
175      REAL vnat(ip1jm,llm),unat(ip1jmp1,llm)
176      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
177      CHARACTER*15 ztit
178!IM   INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
179!IM   SAVE      ip_ebil_dyn
180!IM   DATA      ip_ebil_dyn/0/
181c-jld
182
[37]183      integer :: itau_w ! for write_paramLMDZ_dyn.h
184
[1]185      character*80 dynhist_file, dynhistave_file
186      character(len=20) :: modname
187      character*80 abort_message
188
189      logical dissip_conservative
190      save dissip_conservative
191      data dissip_conservative/.true./
192
193      INTEGER testita
194      PARAMETER (testita = 9)
195
196      logical , parameter :: flag_verif = .false.
197     
[37]198      ! for CP(T)
199      real :: dtec
200      real,external :: cpdet
201      real :: ztetaec(ip1jmp1,llm)
[1]202
203      itaufin   = nday*day_step
204      itaufinp1 = itaufin +1
205      modname="leapfrog"
206     
[6]207c INITIALISATIONS
208        dudis(:,:)   =0.
209        dvdis(:,:)   =0.
210        dtetadis(:,:)=0.
211        dutop(:,:)   =0.
212        dvtop(:,:)   =0.
213        dtetatop(:,:)=0.
214        dqtop(:,:,:) =0.
215        dptop(:)     =0.
216        dufi(:,:)   =0.
217        dvfi(:,:)   =0.
218        dtetafi(:,:)=0.
219        dqfi(:,:,:) =0.
220        dpfi(:)     =0.
[1]221
222      itau = 0
223c      iday = day_ini+itau/day_step
224c      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
225c         IF(time.GT.1.) THEN
226c          time = time-1.
227c          iday = iday+1
228c         ENDIF
229
230
231c-----------------------------------------------------------------------
232c   On initialise la pression et la fonction d'Exner :
233c   --------------------------------------------------
234
[7]235      dq(:,:,:)=0.
[1]236      CALL pression ( ip1jmp1, ap, bp, ps, p       )
237      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
238
239c-----------------------------------------------------------------------
240c   Debut de l'integration temporelle:
241c   ----------------------------------
242
243   1  CONTINUE
244
245      jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
246      jH_cur = jH_ref +                                                 &
247     &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
248
249
250#ifdef CPP_IOIPSL
251      if (ok_guide) then
252        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
253      endif
254#endif
255
256
257c
258c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
259c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
260c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
261c     ENDIF
262c
263
264! Save fields obtained at previous time step as '...m1'
265      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
266      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
267      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
268      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
269      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
270
271      forward = .TRUE.
272      leapf   = .FALSE.
273      dt      =  dtvr
274
275c   ...    P.Le Van .26/04/94  ....
276
277      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
278      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
279
280   2  CONTINUE
281
282c-----------------------------------------------------------------------
283
284c   date:
285c   -----
286
287
288c   gestion des appels de la physique et des dissipations:
289c   ------------------------------------------------------
290c
291c   ...    P.Le Van  ( 6/02/95 )  ....
292
293      apphys = .FALSE.
294      statcl = .FALSE.
295      conser = .FALSE.
296      apdiss = .FALSE.
297
298      IF( purmats ) THEN
299      ! Purely Matsuno time stepping
300         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
301         IF( MOD(itau,idissip ).EQ.0.AND..NOT.forward ) apdiss = .TRUE.
302         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
303     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
304      ELSE
305      ! Leapfrog/Matsuno time stepping
306         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
307         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
308         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
309      END IF
310
311! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
312!          supress dissipation step
313      if (llm.eq.1) then
314        apdiss=.false.
315      endif
316
317c-----------------------------------------------------------------------
318c   calcul des tendances dynamiques:
319c   --------------------------------
320
[5]321! ADAPTATION GCM POUR CP(T)
322      call tpot2t(ijp1llm,teta,temp,pk)
323      tsurpk = cpp*temp/pk
324      CALL geopot  ( ip1jmp1, tsurpk  , pk , pks,  phis  , phi   )
[1]325
326      time = jD_cur + jH_cur
327      CALL caldyn
[5]328     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,tsurpk,phis ,
[1]329     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
330
331
332c-----------------------------------------------------------------------
333c   calcul des tendances advection des traceurs (dont l'humidite)
334c   -------------------------------------------------------------
335
336      IF( forward. OR . leapf )  THEN
337
338         CALL caladvtrac(q,pbaru,pbarv,
339     *        p, masse, dq,  teta,
340     .        flxw, pk)
341         
342         IF (offline) THEN
343Cmaf stokage du flux de masse pour traceurs OFF-LINE
344
345#ifdef CPP_IOIPSL
346           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
347     .   dtvr, itau)
348#endif
349
350
351         ENDIF ! of IF (offline)
352c
353      ENDIF ! of IF( forward. OR . leapf )
354
355
356c-----------------------------------------------------------------------
357c   integrations dynamique et traceurs:
358c   ----------------------------------
359
360
361       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
362     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
363     $              finvmaold                                    )
364
365
366c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
367c
368c-----------------------------------------------------------------------
369c   calcul des tendances physiques:
370c   -------------------------------
371c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
372c
373       IF( purmats )  THEN
374          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
375       ELSE
376          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
377       ENDIF
378c
379c
380       IF( apphys )  THEN
381c
382c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
383c
384
385         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
386         CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
387
388!           rdaym_ini  = itau * dtvr / daysec
389!           rdayvrai   = rdaym_ini  + day_ini
390           jD_cur = jD_ref + day_ini - day_ref
391     $        + int (itau * dtvr / daysec)
392           jH_cur = jH_ref +                                            &
393     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
394!         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
395!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
396!         write(lunout,*)'current date = ',an, mois, jour, secondes
397
398c rajout debug
399c       lafin = .true.
400
401
[6]402c   Interface avec les routines de phylmd (phymars ... )
[1]403c   -----------------------------------------------------
404
405c+jld
406
407c  Diagnostique de conservation de l'énergie : initialisation
408         IF (ip_ebil_dyn.ge.1 ) THEN
409          ztit='bil dyn'
410! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
411           IF (planet_type.eq."earth") THEN
412            CALL diagedyn(ztit,2,1,1,dtphys
413     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
414           ENDIF
415         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
416c-jld
417#ifdef CPP_IOIPSL
418cIM : pour sortir les param. du modele dans un fis. netcdf 110106
419         IF (first) THEN
420#include "ini_paramLMDZ_dyn.h"
421         ENDIF
422c
423#include "write_paramLMDZ_dyn.h"
424c
425#endif
426! #endif of #ifdef CPP_IOIPSL
[6]427
[1]428         CALL calfis( lafin , jD_cur, jH_cur,
429     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
430     $               du,dv,dteta,dq,
431     $               flxw,
432     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
433
[6]434c      ajout des tendances physiques:
435c      ------------------------------
436          CALL addfi( dtphys, leapf, forward   ,
437     $                  ucov, vcov, teta , q   ,ps ,
438     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
439
440c      Couche superieure :
441c      -------------------
[1]442         IF (ok_strato) THEN
[8]443           CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
[1]444         ENDIF
[6]445c dqtop=0, dptop=0
[1]446       
447c      ajout des tendances physiques:
448c      ------------------------------
449          CALL addfi( dtphys, leapf, forward   ,
450     $                  ucov, vcov, teta , q   ,ps ,
[8]451     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
[1]452c
453c  Diagnostique de conservation de l'énergie : difference
454         IF (ip_ebil_dyn.ge.1 ) THEN
455          ztit='bil phys'
456          IF (planet_type.eq."earth") THEN
457           CALL diagedyn(ztit,2,1,1,dtphys
458     &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
459          ENDIF
460         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
461
462       ENDIF ! of IF( apphys )
463
464      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
465!   Academic case : Simple friction and Newtonan relaxation
466!   -------------------------------------------------------
467        DO l=1,llm   
468          DO ij=1,ip1jmp1
469           teta(ij,l)=teta(ij,l)-dtvr*
470     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
471          ENDDO
472        ENDDO ! of DO l=1,llm
473          call friction(ucov,vcov,dtvr)
[53]474   
475       if (planet_type.eq."giant") then
476          ! Intrinsic heat flux
477          ! Aymeric -- for giant planets
478          if (ihf .gt. 1.e-6) then
479          !print *, '**** INTRINSIC HEAT FLUX ****', ihf
480            DO ij=1,ip1jmp1
481              teta(ij,1) = teta(ij,1)
482     &        + dtvr * aire(ij) * ihf / cpp / masse(ij,1)
483            ENDDO
484          !print *, '**** d teta '
485          !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
486          endif
487       endif
488
489   
[1]490        ! Sponge layer (if any)
491        IF (ok_strato) THEN
[7]492          dutop(:,:)=0.
493          dvtop(:,:)=0.
494          dtetatop(:,:)=0.
495          dqtop(:,:,:)=0.
496          dptop(:)=0.
[6]497          CALL top_bound(vcov,ucov,teta,masse,dutop,dvtop,dtetatop)
[1]498          CALL addfi( dtvr, leapf, forward   ,
499     $                  ucov, vcov, teta , q   ,ps ,
[6]500     $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
[1]501        ENDIF ! of IF (ok_strato)
502      ENDIF ! of IF (iflag_phys.EQ.2)
503
504
505c-jld
506
507        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
508        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
509
510
511c-----------------------------------------------------------------------
512c   dissipation horizontale et verticale  des petites echelles:
513c   ----------------------------------------------------------
514
515      IF(apdiss) THEN
516
517
518c   calcul de l'energie cinetique avant dissipation
519        call covcont(llm,ucov,vcov,ucont,vcont)
520        call enercin(vcov,ucov,vcont,ucont,ecin0)
521
522c   dissipation
[5]523! ADAPTATION GCM POUR CP(T)
524        call tpot2t(ijp1llm,teta,temp,pk)
525
[1]526        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
527        ucov=ucov+dudis
528        vcov=vcov+dvdis
[6]529        dudis=dudis/dtdiss   ! passage en (m/s)/s
530        dvdis=dvdis/dtdiss   ! passage en (m/s)/s
[1]531
532c------------------------------------------------------------------------
533        if (dissip_conservative) then
534C       On rajoute la tendance due a la transform. Ec -> E therm. cree
535C       lors de la dissipation
536            call covcont(llm,ucov,vcov,ucont,vcont)
537            call enercin(vcov,ucov,vcont,ucont,ecin)
[5]538! ADAPTATION GCM POUR CP(T)
539            do ij=1,ip1jmp1
540              do l=1,llm
541                dtec = (ecin0(ij,l)-ecin(ij,l))/cpdet(temp(ij,l))
542                temp(ij,l) = temp(ij,l) + dtec
543              enddo
544            enddo
545            call t2tpot(ijp1llm,temp,ztetaec,pk)
546            dtetaecdt=ztetaec-teta
[1]547            dtetadis=dtetadis+dtetaecdt
548        endif
549        teta=teta+dtetadis
[6]550        dtetadis=dtetadis/dtdiss   ! passage en K/s
[1]551c------------------------------------------------------------------------
552
553
554c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
555c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
556c
557
558        DO l  =  1, llm
559          DO ij =  1,iim
560           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
561           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
562          ENDDO
563           tpn  = SSUM(iim,tppn,1)/apoln
564           tps  = SSUM(iim,tpps,1)/apols
565
566          DO ij = 1, iip1
567           teta(  ij    ,l) = tpn
568           teta(ij+ip1jm,l) = tps
569          ENDDO
570        ENDDO
571
572        DO ij =  1,iim
573          tppn(ij)  = aire(  ij    ) * ps (  ij    )
574          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
575        ENDDO
576          tpn  = SSUM(iim,tppn,1)/apoln
577          tps  = SSUM(iim,tpps,1)/apols
578
579        DO ij = 1, iip1
580          ps(  ij    ) = tpn
581          ps(ij+ip1jm) = tps
582        ENDDO
583
584
585      END IF ! of IF(apdiss)
586
587c ajout debug
588c              IF( lafin ) then 
589c                abort_message = 'Simulation finished'
590c                call abort_gcm(modname,abort_message,0)
591c              ENDIF
592       
593c   ********************************************************************
594c   ********************************************************************
595c   .... fin de l'integration dynamique  et physique pour le pas itau ..
596c   ********************************************************************
597c   ********************************************************************
598
599c   preparation du pas d'integration suivant  ......
600
601      IF ( .NOT.purmats ) THEN
602c       ........................................................
603c       ..............  schema matsuno + leapfrog  ..............
604c       ........................................................
605
606            IF(forward. OR. leapf) THEN
607              itau= itau + 1
608c              iday= day_ini+itau/day_step
609c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
610c                IF(time.GT.1.) THEN
611c                  time = time-1.
612c                  iday = iday+1
613c                ENDIF
614            ENDIF
615
616
617            IF( itau. EQ. itaufinp1 ) then 
618              if (flag_verif) then
619                write(79,*) 'ucov',ucov
620                write(80,*) 'vcov',vcov
621                write(81,*) 'teta',teta
622                write(82,*) 'ps',ps
623                write(83,*) 'q',q
624                WRITE(85,*) 'q1 = ',q(:,:,1)
625                WRITE(86,*) 'q3 = ',q(:,:,3)
626              endif
627
628              abort_message = 'Simulation finished'
629
630              call abort_gcm(modname,abort_message,0)
631            ENDIF
632c-----------------------------------------------------------------------
633c   ecriture du fichier histoire moyenne:
634c   -------------------------------------
635
636            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
637               IF(itau.EQ.itaufin) THEN
638                  iav=1
639               ELSE
640                  iav=0
641               ENDIF
642               
643               IF (ok_dynzon) THEN
644#ifdef CPP_IOIPSL
[6]645c les traceurs ne sont pas sortis, trop lourd.
646c Peut changer eventuellement si besoin.
647                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
648     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
649     &                 du,dudis,duspg,dufi)
[1]650#endif
651               END IF
652               IF (ok_dyn_ave) THEN
653#ifdef CPP_IOIPSL
654                 CALL writedynav(itau,vcov,
655     &                 ucov,teta,pk,phi,q,masse,ps,phis)
656#endif
657               ENDIF
658
659            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
660
661c-----------------------------------------------------------------------
662c   ecriture de la bande histoire:
663c   ------------------------------
664
665            IF( MOD(itau,iecri).EQ.0) THEN
666             ! Ehouarn: output only during LF or Backward Matsuno
667             if (leapf.or.(.not.leapf.and.(.not.forward))) then
668              nbetat = nbetatdem
[5]669! ADAPTATION GCM POUR CP(T)
670              call tpot2t(ijp1llm,teta,temp,pk)
671              tsurpk = cpp*temp/pk
672              CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1]673              unat=0.
674              do l=1,llm
675                unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
676                vnat(:,l)=vcov(:,l)/cv(:)
677              enddo
678#ifdef CPP_IOIPSL
679              if (ok_dyn_ins) then
680!               write(lunout,*) "leapfrog: call writehist, itau=",itau
681               CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
682!               call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
683!               call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
684!              call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
685!               call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
686!               call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
687              endif ! of if (ok_dyn_ins)
688#endif
689! For some Grads outputs of fields
690              if (output_grads_dyn) then
691#include "write_grads_dyn.h"
692              endif
693             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
694            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
695
696            IF(itau.EQ.itaufin) THEN
697
698
[6]699              if (planet_type.eq."mars") then
700! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
701                abort_message = 'dynredem1_mars A FAIRE'
702                call abort_gcm(modname,abort_message,0)
703              else
[1]704                CALL dynredem1("restart.nc",0.0,
705     &                         vcov,ucov,teta,q,masse,ps)
[6]706              endif ! of if (planet_type.eq."mars")
[1]707
708              CLOSE(99)
709            ENDIF ! of IF (itau.EQ.itaufin)
710
711c-----------------------------------------------------------------------
712c   gestion de l'integration temporelle:
713c   ------------------------------------
714
715            IF( MOD(itau,iperiod).EQ.0 )    THEN
716                    GO TO 1
717            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
718
719                   IF( forward )  THEN
720c      fin du pas forward et debut du pas backward
721
722                      forward = .FALSE.
723                        leapf = .FALSE.
724                           GO TO 2
725
726                   ELSE
727c      fin du pas backward et debut du premier pas leapfrog
728
729                        leapf =  .TRUE.
730                        dt  =  2.*dtvr
731                        GO TO 2
732                   END IF ! of IF (forward)
733            ELSE
734
735c      ......   pas leapfrog  .....
736
737                 leapf = .TRUE.
738                 dt  = 2.*dtvr
739                 GO TO 2
740            END IF ! of IF (MOD(itau,iperiod).EQ.0)
741                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
742
743      ELSE ! of IF (.not.purmats)
744
745c       ........................................................
746c       ..............       schema  matsuno        ...............
747c       ........................................................
748            IF( forward )  THEN
749
750             itau =  itau + 1
751c             iday = day_ini+itau/day_step
752c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
753c
754c                  IF(time.GT.1.) THEN
755c                   time = time-1.
756c                   iday = iday+1
757c                  ENDIF
758
759               forward =  .FALSE.
760               IF( itau. EQ. itaufinp1 ) then 
761                 abort_message = 'Simulation finished'
762                 call abort_gcm(modname,abort_message,0)
763               ENDIF
764               GO TO 2
765
766            ELSE ! of IF(forward) i.e. backward step
767
768              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
769               IF(itau.EQ.itaufin) THEN
770                  iav=1
771               ELSE
772                  iav=0
773               ENDIF
774
775               IF (ok_dynzon) THEN
776#ifdef CPP_IOIPSL
[6]777c les traceurs ne sont pas sortis, trop lourd.
778c Peut changer eventuellement si besoin.
779                 CALL bilan_dyn(dtvr*iperiod,dtvr*day_step*periodav,
780     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,
781     &                 du,dudis,duspg,dufi)
[1]782#endif
783               ENDIF
784               IF (ok_dyn_ave) THEN
785#ifdef CPP_IOIPSL
786                 CALL writedynav(itau,vcov,
787     &                 ucov,teta,pk,phi,q,masse,ps,phis)
788#endif
789               ENDIF
790
791              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
792
793              IF(MOD(itau,iecri         ).EQ.0) THEN
794c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
795                nbetat = nbetatdem
[5]796! ADAPTATION GCM POUR CP(T)
797                call tpot2t(ijp1llm,teta,temp,pk)
798                tsurpk = cpp*temp/pk
799                CALL geopot(ip1jmp1,tsurpk,pk,pks,phis,phi)
[1]800                unat=0.
801                do l=1,llm
802                  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
803                  vnat(:,l)=vcov(:,l)/cv(:)
804                enddo
805#ifdef CPP_IOIPSL
806              if (ok_dyn_ins) then
807!                write(lunout,*) "leapfrog: call writehist (b)",
808!     &                        itau,iecri
809                CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
810              endif ! of if (ok_dyn_ins)
811#endif
812! For some Grads outputs
813                if (output_grads_dyn) then
814#include "write_grads_dyn.h"
815                endif
816
817              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
818
819              IF(itau.EQ.itaufin) THEN
[6]820                if (planet_type.eq."mars") then
821! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem1_mars
822                  abort_message = 'dynredem1_mars A FAIRE'
823                  call abort_gcm(modname,abort_message,0)
824                else
[1]825                  CALL dynredem1("restart.nc",0.0,
[6]826     &                         vcov,ucov,teta,q,masse,ps)
827                endif ! of if (planet_type.eq."mars")
[1]828              ENDIF ! of IF(itau.EQ.itaufin)
829
830              forward = .TRUE.
831              GO TO  1
832
833            ENDIF ! of IF (forward)
834
835      END IF ! of IF(.not.purmats)
836
[6]837      first=.false.
838
[1]839      STOP
840      END
Note: See TracBrowser for help on using the repository browser.