source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/leapfrog.F @ 3990

Last change on this file since 3990 was 1529, checked in by Laurent Fairhead, 13 years ago

Necessary modifications to use the model in an aquaplanet or terraplanet configuration


Modifications nécessaires pour l'utilisation du moèle en configuration aquaplanète
ou terraplanète

Dans la dynamique:
Changement pour le paramètre iflag_phys:

auparavant : 0 pas de physique, 1 phyique, 2 rappel
ici on ajoute iflag_phys>=100 pour les aqua et terra

Du coup, dans leapfrog.F et gcm.F on appelle la physique si
iflag_phys=1.or.iflag_phys>=100
Dans iniacademic, on initialise si iflag_phys>=2 au lieu de =2
Dans gcm.F, on appelle en plus iniaqua (sous une clef CPP_EARTH)
Dans iniacademic, on met ps=108080 pour les aqua et terra pour répondre
à une specification CFMIP.

Dans la physique:
On ajoute phyaqua.F qui contient :
iniaqua : initialise les startphy.nc et limit.nc pour la phyique
zenang_an : calcule un ensolleillement moyen sur l'année en fonction de

la latiude (A. Campoy).

profil_sst : qui calcule différents profils latitudinaux de SSTs.
writelim : écrit le fichier limit.nc

Dans physiq.F
On ajoute la possibilité d'appeller un calcul d'ensoleillement moyen
sur l'année quand solarlong0=1000.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 24.4 KB
Line 
1!
2! $Id: leapfrog.F 1529 2011-05-26 15:17:33Z 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      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
89
90c variables dynamiques intermediaire pour le transport
91      REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
92
93c   variables dynamiques au pas -1
94      REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
95      REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
96      REAL massem1(ip1jmp1,llm)
97
98c   tendances dynamiques
99      REAL dv(ip1jm,llm),du(ip1jmp1,llm)
100      REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
101
102c   tendances de la dissipation
103      REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
104      REAL dtetadis(ip1jmp1,llm)
105
106c   tendances physiques
107      REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
108      REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
109
110c   variables pour le fichier histoire
111      REAL dtav      ! intervalle de temps elementaire
112
113      REAL tppn(iim),tpps(iim),tpn,tps
114c
115      INTEGER itau,itaufinp1,iav
116!      INTEGER  iday ! jour julien
117      REAL       time
118
119      REAL  SSUM
120      REAL time_0 , 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 (disvert_type==1) then
215        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
216      else ! we assume that we are in the disvert_type==2 case
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
225
226      jD_cur = jD_ref + day_ini - day_ref + int (itau * dtvr / daysec)
227      jH_cur = jH_ref +                                                 &
228     &          (itau * dtvr / daysec - int(itau * dtvr / daysec))
229
230
231#ifdef CPP_IOIPSL
232      if (ok_guide) then
233        call guide_main(itau,ucov,vcov,teta,q,masse,ps)
234      endif
235#endif
236
237
238c
239c     IF( MOD( itau, 10* day_step ).EQ.0 )  THEN
240c       CALL  test_period ( ucov,vcov,teta,q,p,phis )
241c       PRINT *,' ----   Test_period apres continue   OK ! -----', itau
242c     ENDIF
243c
244
245! Save fields obtained at previous time step as '...m1'
246      CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
247      CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
248      CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
249      CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
250      CALL SCOPY( ip1jmp1, ps  , 1,   psm1 , 1 )
251
252      forward = .TRUE.
253      leapf   = .FALSE.
254      dt      =  dtvr
255
256c   ...    P.Le Van .26/04/94  ....
257
258      CALL SCOPY   ( ijp1llm,   masse, 1, finvmaold,     1 )
259      CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
260
261   2  CONTINUE
262
263c-----------------------------------------------------------------------
264
265c   date:
266c   -----
267
268
269c   gestion des appels de la physique et des dissipations:
270c   ------------------------------------------------------
271c
272c   ...    P.Le Van  ( 6/02/95 )  ....
273
274      apphys = .FALSE.
275      statcl = .FALSE.
276      conser = .FALSE.
277      apdiss = .FALSE.
278
279      IF( purmats ) THEN
280      ! Purely Matsuno time stepping
281         IF( MOD(itau,iconser) .EQ.0.AND.  forward    ) conser = .TRUE.
282         IF( MOD(itau,dissip_period ).EQ.0.AND..NOT.forward )
283     s        apdiss = .TRUE.
284         IF( MOD(itau,iphysiq ).EQ.0.AND..NOT.forward
285     s          .and. physic                        ) apphys = .TRUE.
286      ELSE
287      ! Leapfrog/Matsuno time stepping
288         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
289         IF( MOD(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
290     s        apdiss = .TRUE.
291         IF( MOD(itau+1,iphysiq).EQ.0.AND.physic       ) apphys=.TRUE.
292      END IF
293
294! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
295!          supress dissipation step
296      if (llm.eq.1) then
297        apdiss=.false.
298      endif
299
300c-----------------------------------------------------------------------
301c   calcul des tendances dynamiques:
302c   --------------------------------
303
304      CALL geopot  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
305
306      time = jD_cur + jH_cur
307      CALL caldyn
308     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
309     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
310
311
312c-----------------------------------------------------------------------
313c   calcul des tendances advection des traceurs (dont l'humidite)
314c   -------------------------------------------------------------
315
316      IF( forward. OR . leapf )  THEN
317
318         CALL caladvtrac(q,pbaru,pbarv,
319     *        p, masse, dq,  teta,
320     .        flxw, pk)
321         
322         IF (offline) THEN
323Cmaf stokage du flux de masse pour traceurs OFF-LINE
324
325#ifdef CPP_IOIPSL
326           CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
327     .   dtvr, itau)
328#endif
329
330
331         ENDIF ! of IF (offline)
332c
333      ENDIF ! of IF( forward. OR . leapf )
334
335
336c-----------------------------------------------------------------------
337c   integrations dynamique et traceurs:
338c   ----------------------------------
339
340
341       CALL integrd ( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
342     $         dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis ,
343     $              finvmaold                                    )
344
345
346c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
347c
348c-----------------------------------------------------------------------
349c   calcul des tendances physiques:
350c   -------------------------------
351c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
352c
353       IF( purmats )  THEN
354          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
355       ELSE
356          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
357       ENDIF
358c
359c
360       IF( apphys )  THEN
361c
362c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
363c
364
365         CALL pression (  ip1jmp1, ap, bp, ps,  p      )
366         if (disvert_type==1) then
367           CALL exner_hyb(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
368         else ! we assume that we are in the disvert_type==2 case
369           CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
370         endif
371
372!           rdaym_ini  = itau * dtvr / daysec
373!           rdayvrai   = rdaym_ini  + day_ini
374           jD_cur = jD_ref + day_ini - day_ref
375     $        + int (itau * dtvr / daysec)
376           jH_cur = jH_ref +                                            &
377     &              (itau * dtvr / daysec - int(itau * dtvr / daysec))
378!         write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
379!         call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
380!         write(lunout,*)'current date = ',an, mois, jour, secondes
381
382c rajout debug
383c       lafin = .true.
384
385
386c   Inbterface avec les routines de phylmd (phymars ... )
387c   -----------------------------------------------------
388
389c+jld
390
391c  Diagnostique de conservation de l'énergie : initialisation
392         IF (ip_ebil_dyn.ge.1 ) THEN
393          ztit='bil dyn'
394! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
395           IF (planet_type.eq."earth") THEN
396            CALL diagedyn(ztit,2,1,1,dtphys
397     &    , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
398           ENDIF
399         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
400c-jld
401#ifdef CPP_IOIPSL
402cIM : pour sortir les param. du modele dans un fis. netcdf 110106
403         IF (first) THEN
404          first=.false.
405#include "ini_paramLMDZ_dyn.h"
406         ENDIF
407c
408#include "write_paramLMDZ_dyn.h"
409c
410#endif
411! #endif of #ifdef CPP_IOIPSL
412         CALL calfis( lafin , jD_cur, jH_cur,
413     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
414     $               du,dv,dteta,dq,
415     $               flxw,
416     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
417
418         IF (ok_strato) THEN
419           CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
420         ENDIF
421       
422c      ajout des tendances physiques:
423c      ------------------------------
424          CALL addfi( dtphys, leapf, forward   ,
425     $                  ucov, vcov, teta , q   ,ps ,
426     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
427c
428c  Diagnostique de conservation de l'énergie : difference
429         IF (ip_ebil_dyn.ge.1 ) THEN
430          ztit='bil phys'
431          IF (planet_type.eq."earth") THEN
432           CALL diagedyn(ztit,2,1,1,dtphys
433     &     , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
434          ENDIF
435         ENDIF ! of IF (ip_ebil_dyn.ge.1 )
436
437       ENDIF ! of IF( apphys )
438
439      IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
440!   Academic case : Simple friction and Newtonan relaxation
441!   -------------------------------------------------------
442        DO l=1,llm   
443          DO ij=1,ip1jmp1
444           teta(ij,l)=teta(ij,l)-dtvr*
445     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
446          ENDDO
447        ENDDO ! of DO l=1,llm
448       
449        if (planet_type.eq."giant") then
450          ! add an intrinsic heat flux at the base of the atmosphere
451          teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
452        endif
453
454        call friction(ucov,vcov,dtvr)
455       
456        ! Sponge layer (if any)
457        IF (ok_strato) THEN
458          dufi(:,:)=0.
459          dvfi(:,:)=0.
460          dtetafi(:,:)=0.
461          dqfi(:,:,:)=0.
462          dpfi(:)=0.
463          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
464          CALL addfi( dtvr, leapf, forward   ,
465     $                  ucov, vcov, teta , q   ,ps ,
466     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
467        ENDIF ! of IF (ok_strato)
468      ENDIF ! of IF (iflag_phys.EQ.2)
469
470
471c-jld
472
473        CALL pression ( ip1jmp1, ap, bp, ps, p                  )
474        if (disvert_type==1) then
475          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
476        else ! we assume that we are in the disvert_type==2 case
477          CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
478        endif
479
480
481c-----------------------------------------------------------------------
482c   dissipation horizontale et verticale  des petites echelles:
483c   ----------------------------------------------------------
484
485      IF(apdiss) THEN
486
487
488c   calcul de l'energie cinetique avant dissipation
489        call covcont(llm,ucov,vcov,ucont,vcont)
490        call enercin(vcov,ucov,vcont,ucont,ecin0)
491
492c   dissipation
493        CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
494        ucov=ucov+dudis
495        vcov=vcov+dvdis
496c       teta=teta+dtetadis
497
498
499c------------------------------------------------------------------------
500        if (dissip_conservative) then
501C       On rajoute la tendance due a la transform. Ec -> E therm. cree
502C       lors de la dissipation
503            call covcont(llm,ucov,vcov,ucont,vcont)
504            call enercin(vcov,ucov,vcont,ucont,ecin)
505            dtetaecdt= (ecin0-ecin)/ pk
506c           teta=teta+dtetaecdt
507            dtetadis=dtetadis+dtetaecdt
508        endif
509        teta=teta+dtetadis
510c------------------------------------------------------------------------
511
512
513c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
514c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
515c
516
517        DO l  =  1, llm
518          DO ij =  1,iim
519           tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
520           tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
521          ENDDO
522           tpn  = SSUM(iim,tppn,1)/apoln
523           tps  = SSUM(iim,tpps,1)/apols
524
525          DO ij = 1, iip1
526           teta(  ij    ,l) = tpn
527           teta(ij+ip1jm,l) = tps
528          ENDDO
529        ENDDO
530
531        DO ij =  1,iim
532          tppn(ij)  = aire(  ij    ) * ps (  ij    )
533          tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
534        ENDDO
535          tpn  = SSUM(iim,tppn,1)/apoln
536          tps  = SSUM(iim,tpps,1)/apols
537
538        DO ij = 1, iip1
539          ps(  ij    ) = tpn
540          ps(ij+ip1jm) = tps
541        ENDDO
542
543
544      END IF ! of IF(apdiss)
545
546c ajout debug
547c              IF( lafin ) then 
548c                abort_message = 'Simulation finished'
549c                call abort_gcm(modname,abort_message,0)
550c              ENDIF
551       
552c   ********************************************************************
553c   ********************************************************************
554c   .... fin de l'integration dynamique  et physique pour le pas itau ..
555c   ********************************************************************
556c   ********************************************************************
557
558c   preparation du pas d'integration suivant  ......
559
560      IF ( .NOT.purmats ) THEN
561c       ........................................................
562c       ..............  schema matsuno + leapfrog  ..............
563c       ........................................................
564
565            IF(forward. OR. leapf) THEN
566              itau= itau + 1
567c              iday= day_ini+itau/day_step
568c              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
569c                IF(time.GT.1.) THEN
570c                  time = time-1.
571c                  iday = iday+1
572c                ENDIF
573            ENDIF
574
575
576            IF( itau. EQ. itaufinp1 ) then 
577              if (flag_verif) then
578                write(79,*) 'ucov',ucov
579                write(80,*) 'vcov',vcov
580                write(81,*) 'teta',teta
581                write(82,*) 'ps',ps
582                write(83,*) 'q',q
583                WRITE(85,*) 'q1 = ',q(:,:,1)
584                WRITE(86,*) 'q3 = ',q(:,:,3)
585              endif
586
587              abort_message = 'Simulation finished'
588
589              call abort_gcm(modname,abort_message,0)
590            ENDIF
591c-----------------------------------------------------------------------
592c   ecriture du fichier histoire moyenne:
593c   -------------------------------------
594
595            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
596               IF(itau.EQ.itaufin) THEN
597                  iav=1
598               ELSE
599                  iav=0
600               ENDIF
601               
602               IF (ok_dynzon) THEN
603#ifdef CPP_IOIPSL
604                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
605     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
606#endif
607               END IF
608               IF (ok_dyn_ave) THEN
609#ifdef CPP_IOIPSL
610                 CALL writedynav(itau,vcov,
611     &                 ucov,teta,pk,phi,q,masse,ps,phis)
612#endif
613               ENDIF
614
615            ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
616
617c-----------------------------------------------------------------------
618c   ecriture de la bande histoire:
619c   ------------------------------
620
621            IF( MOD(itau,iecri).EQ.0) THEN
622             ! Ehouarn: output only during LF or Backward Matsuno
623             if (leapf.or.(.not.leapf.and.(.not.forward))) then
624              nbetat = nbetatdem
625              CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
626              unat=0.
627              do l=1,llm
628                unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
629                vnat(:,l)=vcov(:,l)/cv(:)
630              enddo
631#ifdef CPP_IOIPSL
632              if (ok_dyn_ins) then
633!               write(lunout,*) "leapfrog: call writehist, itau=",itau
634               CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
635!               call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
636!               call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
637!              call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
638!               call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
639!               call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
640              endif ! of if (ok_dyn_ins)
641#endif
642! For some Grads outputs of fields
643              if (output_grads_dyn) then
644#include "write_grads_dyn.h"
645              endif
646             endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
647            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
648
649            IF(itau.EQ.itaufin) THEN
650
651
652!              if (planet_type.eq."earth") then
653! Write an Earth-format restart file
654                CALL dynredem1("restart.nc",0.0,
655     &                         vcov,ucov,teta,q,masse,ps)
656!              endif ! of if (planet_type.eq."earth")
657
658              CLOSE(99)
659            ENDIF ! of IF (itau.EQ.itaufin)
660
661c-----------------------------------------------------------------------
662c   gestion de l'integration temporelle:
663c   ------------------------------------
664
665            IF( MOD(itau,iperiod).EQ.0 )    THEN
666                    GO TO 1
667            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
668
669                   IF( forward )  THEN
670c      fin du pas forward et debut du pas backward
671
672                      forward = .FALSE.
673                        leapf = .FALSE.
674                           GO TO 2
675
676                   ELSE
677c      fin du pas backward et debut du premier pas leapfrog
678
679                        leapf =  .TRUE.
680                        dt  =  2.*dtvr
681                        GO TO 2
682                   END IF ! of IF (forward)
683            ELSE
684
685c      ......   pas leapfrog  .....
686
687                 leapf = .TRUE.
688                 dt  = 2.*dtvr
689                 GO TO 2
690            END IF ! of IF (MOD(itau,iperiod).EQ.0)
691                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
692
693      ELSE ! of IF (.not.purmats)
694
695c       ........................................................
696c       ..............       schema  matsuno        ...............
697c       ........................................................
698            IF( forward )  THEN
699
700             itau =  itau + 1
701c             iday = day_ini+itau/day_step
702c             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
703c
704c                  IF(time.GT.1.) THEN
705c                   time = time-1.
706c                   iday = iday+1
707c                  ENDIF
708
709               forward =  .FALSE.
710               IF( itau. EQ. itaufinp1 ) then 
711                 abort_message = 'Simulation finished'
712                 call abort_gcm(modname,abort_message,0)
713               ENDIF
714               GO TO 2
715
716            ELSE ! of IF(forward) i.e. backward step
717
718              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
719               IF(itau.EQ.itaufin) THEN
720                  iav=1
721               ELSE
722                  iav=0
723               ENDIF
724
725               IF (ok_dynzon) THEN
726#ifdef CPP_IOIPSL
727                 CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
728     &                 ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
729#endif
730               ENDIF
731               IF (ok_dyn_ave) THEN
732#ifdef CPP_IOIPSL
733                 CALL writedynav(itau,vcov,
734     &                 ucov,teta,pk,phi,q,masse,ps,phis)
735#endif
736               ENDIF
737
738              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
739
740              IF(MOD(itau,iecri         ).EQ.0) THEN
741c              IF(MOD(itau,iecri*day_step).EQ.0) THEN
742                nbetat = nbetatdem
743                CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
744                unat=0.
745                do l=1,llm
746                  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
747                  vnat(:,l)=vcov(:,l)/cv(:)
748                enddo
749#ifdef CPP_IOIPSL
750              if (ok_dyn_ins) then
751!                write(lunout,*) "leapfrog: call writehist (b)",
752!     &                        itau,iecri
753                CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
754              endif ! of if (ok_dyn_ins)
755#endif
756! For some Grads outputs
757                if (output_grads_dyn) then
758#include "write_grads_dyn.h"
759                endif
760
761              ENDIF ! of IF(MOD(itau,iecri         ).EQ.0)
762
763              IF(itau.EQ.itaufin) THEN
764!                if (planet_type.eq."earth") then
765                  CALL dynredem1("restart.nc",0.0,
766     &                           vcov,ucov,teta,q,masse,ps)
767!                endif ! of if (planet_type.eq."earth")
768              ENDIF ! of IF(itau.EQ.itaufin)
769
770              forward = .TRUE.
771              GO TO  1
772
773            ENDIF ! of IF (forward)
774
775      END IF ! of IF(.not.purmats)
776
777      STOP
778      END
Note: See TracBrowser for help on using the repository browser.