My Project
 All Classes Files Functions Variables Macros
leapfrog.F
Go to the documentation of this file.
1 !
2 ! $Id: leapfrog.F 1676 2012-11-08 20:03:30Z aslmd $
3 !
4 c
5 c
6  SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
7  & time_0)
8 
9 
10 cIM : 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 
20 c ...... Version du 10/01/98 ..........
21 
22 c avec coordonnees verticales hybrides
23 c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
24 
25 c=======================================================================
26 c
27 c Auteur: P. Le Van /L. Fairhead/F.Hourdin
28 c -------
29 c
30 c Objet:
31 c ------
32 c
33 c GCM LMD nouvelle grille
34 c
35 c=======================================================================
36 c
37 c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
38 c et possibilite d'appeler une fonction f(y) a derivee tangente
39 c hyperbolique a la place de la fonction a derivee sinusoidale.
40 
41 c ... Possibilite de choisir le shema pour l'advection de
42 c q , en modifiant iadv dans traceur.def (10/02) .
43 c
44 c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
45 c Pour Van-Leer iadv=10
46 c
47 c-----------------------------------------------------------------------
48 c Declarations:
49 c -------------
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 
75 c 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 
89 c variables dynamiques intermediaire pour le transport
90  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
91 
92 c 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 
97 c tendances dynamiques
98  REAL dv(ip1jm,llm),du(ip1jmp1,llm)
99  REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
100 
101 c tendances de la dissipation
102  REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
103  REAL dtetadis(ip1jmp1,llm)
104 
105 c tendances physiques
106  REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
107  REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
108 
109 c variables pour le fichier histoire
110  REAL dtav ! intervalle de temps elementaire
111 
112  REAL tppn(iim),tpps(iim),tpn,tps
113 c
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 
122 cym 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
137 cIM : 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 
159 c+jld variables test conservation energie
160  REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
161 C Tendance de la temp. potentiel d (theta)/ d t due a la
162 C tansformation d'energie cinetique en energie thermique
163 C 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/
172 c-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 
200 c iday = day_ini+itau/day_step
201 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
202 c IF(time.GT.1.) THEN
203 c time = time-1.
204 c iday = iday+1
205 c ENDIF
206 
207 
208 c-----------------------------------------------------------------------
209 c On initialise la pression et la fonction d'Exner :
210 c --------------------------------------------------
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 
220 c-----------------------------------------------------------------------
221 c Debut de l'integration temporelle:
222 c ----------------------------------
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 
241 c
242 c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
243 c CALL test_period ( ucov,vcov,teta,q,p,phis )
244 c PRINT *,' ---- Test_period apres continue OK ! -----', itau
245 c ENDIF
246 c
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 
259 c ... 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 
266 c-----------------------------------------------------------------------
267 
268 c date:
269 c -----
270 
271 
272 c gestion des appels de la physique et des dissipations:
273 c ------------------------------------------------------
274 c
275 c ... 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 
303 c-----------------------------------------------------------------------
304 c calcul des tendances dynamiques:
305 c --------------------------------
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 
316 c-----------------------------------------------------------------------
317 c calcul des tendances advection des traceurs (dont l'humidite)
318 c -------------------------------------------------------------
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
328 Cmaf 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)
337 c
338  ENDIF ! of IF( forward. OR . leapf )
339 
340 
341 c-----------------------------------------------------------------------
342 c integrations dynamique et traceurs:
343 c ----------------------------------
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 
351 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
352 c
353 c-----------------------------------------------------------------------
354 c calcul des tendances physiques:
355 c -------------------------------
356 c ######## P.Le Van ( Modif le 6/02/95 ) ###########
357 c
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
363 c
364 c
365  IF( apphys ) THEN
366 c
367 c ....... Ajout P.Le Van ( 17/04/96 ) ...........
368 c
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 
399 c rajout debug
400 c lafin = .true.
401 
402 
403 c Inbterface avec les routines de phylmd (phymars ... )
404 c -----------------------------------------------------
405 
406 c+jld
407 
408 c 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 )
419 c-jld
420 #ifdef CPP_IOIPSL
421 cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
422 cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
423 c IF (first) THEN
424 c first=.false.
425 c#include "ini_paramLMDZ_dyn.h"
426 c ENDIF
427 c
428 c#include "write_paramLMDZ_dyn.h"
429 c
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 
438  IF (ok_strato) THEN
439  CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
440  ENDIF
441 
442 c ajout des tendances physiques:
443 c ------------------------------
444  CALL addfi( dtphys, leapf, forward ,
445  $ ucov, vcov, teta , q ,ps ,
446  $ dufi, dvfi, dtetafi , dqfi ,dpfi )
447 c
448 c Diagnostique de conservation de l'énergie : difference
449  IF (ip_ebil_dyn.ge.1 ) THEN
450  ztit='bil phys'
451  IF (planet_type.eq."earth") THEN
452  CALL diagedyn(ztit,2,1,1,dtphys
453  & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
454  ENDIF
455  ENDIF ! of IF (ip_ebil_dyn.ge.1 )
456 
457  ENDIF ! of IF( apphys )
458 
459  IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
460 ! Academic case : Simple friction and Newtonan relaxation
461 ! -------------------------------------------------------
462  DO l=1,llm
463  DO ij=1,ip1jmp1
464  teta(ij,l)=teta(ij,l)-dtvr*
465  & (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
466  ENDDO
467  ENDDO ! of DO l=1,llm
468 
469  if (planet_type.eq."giant") then
470  ! add an intrinsic heat flux at the base of the atmosphere
471  teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
472  endif
473 
474  call friction(ucov,vcov,dtvr)
475 
476  ! Sponge layer (if any)
477  IF (ok_strato) THEN
478  dufi(:,:)=0.
479  dvfi(:,:)=0.
480  dtetafi(:,:)=0.
481  dqfi(:,:,:)=0.
482  dpfi(:)=0.
483  CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
484  CALL addfi( dtvr, leapf, forward ,
485  $ ucov, vcov, teta , q ,ps ,
486  $ dufi, dvfi, dtetafi , dqfi ,dpfi )
487  ENDIF ! of IF (ok_strato)
488  ENDIF ! of IF (iflag_phys.EQ.2)
489 
490 
491 c-jld
492 
493  CALL pression( ip1jmp1, ap, bp, ps, p )
494  if (pressure_exner) then
495  CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
496  else
497  CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
498  endif
499 
500 
501 c-----------------------------------------------------------------------
502 c dissipation horizontale et verticale des petites echelles:
503 c ----------------------------------------------------------
504 
505  IF(apdiss) THEN
506 
507 
508 c calcul de l'energie cinetique avant dissipation
509  call covcont(llm,ucov,vcov,ucont,vcont)
510  call enercin(vcov,ucov,vcont,ucont,ecin0)
511 
512 c dissipation
513  CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
514  ucov=ucov+dudis
515  vcov=vcov+dvdis
516 c teta=teta+dtetadis
517 
518 
519 c------------------------------------------------------------------------
520  if (dissip_conservative) then
521 C On rajoute la tendance due a la transform. Ec -> E therm. cree
522 C lors de la dissipation
523  call covcont(llm,ucov,vcov,ucont,vcont)
524  call enercin(vcov,ucov,vcont,ucont,ecin)
525  dtetaecdt= (ecin0-ecin)/ pk
526 c teta=teta+dtetaecdt
527  dtetadis=dtetadis+dtetaecdt
528  endif
529  teta=teta+dtetadis
530 c------------------------------------------------------------------------
531 
532 
533 c ....... P. Le Van ( ajout le 17/04/96 ) ...........
534 c ... Calcul de la valeur moyenne, unique de h aux poles .....
535 c
536 
537  DO l = 1, llm
538  DO ij = 1,iim
539  tppn(ij) = aire( ij ) * teta( ij ,l)
540  tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
541  ENDDO
542  tpn = ssum(iim,tppn,1)/apoln
543  tps = ssum(iim,tpps,1)/apols
544 
545  DO ij = 1, iip1
546  teta( ij ,l) = tpn
547  teta(ij+ip1jm,l) = tps
548  ENDDO
549  ENDDO
550 
551  if (1 == 0) then
552 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
553 !!! 2) should probably not be here anyway
554 !!! but are kept for those who would want to revert to previous behaviour
555  DO ij = 1,iim
556  tppn(ij) = aire( ij ) * ps( ij )
557  tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
558  ENDDO
559  tpn = ssum(iim,tppn,1)/apoln
560  tps = ssum(iim,tpps,1)/apols
561 
562  DO ij = 1, iip1
563  ps( ij ) = tpn
564  ps(ij+ip1jm) = tps
565  ENDDO
566  endif ! of if (1 == 0)
567 
568  END IF ! of IF(apdiss)
569 
570 c ajout debug
571 c IF( lafin ) then
572 c abort_message = 'Simulation finished'
573 c call abort_gcm(modname,abort_message,0)
574 c ENDIF
575 
576 c ********************************************************************
577 c ********************************************************************
578 c .... fin de l'integration dynamique et physique pour le pas itau ..
579 c ********************************************************************
580 c ********************************************************************
581 
582 c preparation du pas d'integration suivant ......
583 
584  IF ( .NOT.purmats ) THEN
585 c ........................................................
586 c .............. schema matsuno + leapfrog ..............
587 c ........................................................
588 
589  IF(forward. or. leapf) THEN
590  itau= itau + 1
591 c iday= day_ini+itau/day_step
592 c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
593 c IF(time.GT.1.) THEN
594 c time = time-1.
595 c iday = iday+1
596 c ENDIF
597  ENDIF
598 
599 
600  IF( itau. eq. itaufinp1 ) then
601  if (flag_verif) then
602  write(79,*) 'ucov',ucov
603  write(80,*) 'vcov',vcov
604  write(81,*) 'teta',teta
605  write(82,*) 'ps',ps
606  write(83,*) 'q',q
607  WRITE(85,*) 'q1 = ',q(:,:,1)
608  WRITE(86,*) 'q3 = ',q(:,:,3)
609  endif
610 
611  abort_message = 'Simulation finished'
612 
613  call abort_gcm(modname,abort_message,0)
614  ENDIF
615 c-----------------------------------------------------------------------
616 c ecriture du fichier histoire moyenne:
617 c -------------------------------------
618 
619  IF(mod(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
620  IF(itau.EQ.itaufin) THEN
621  iav=1
622  ELSE
623  iav=0
624  ENDIF
625 
626  IF (ok_dynzon) THEN
627 #ifdef CPP_IOIPSL
628  CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
629  & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
630 #endif
631  END IF
632  IF (ok_dyn_ave) THEN
633 #ifdef CPP_IOIPSL
634  CALL writedynav(itau,vcov,
635  & ucov,teta,pk,phi,q,masse,ps,phis)
636 #endif
637  ENDIF
638 
639  ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
640 
641 c-----------------------------------------------------------------------
642 c ecriture de la bande histoire:
643 c ------------------------------
644 
645  IF( mod(itau,iecri).EQ.0) THEN
646  ! Ehouarn: output only during LF or Backward Matsuno
647  if (leapf.or.(.not.leapf.and.(.not.forward))) then
648  CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
649  unat=0.
650  do l=1,llm
651  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
652  vnat(:,l)=vcov(:,l)/cv(:)
653  enddo
654 #ifdef CPP_IOIPSL
655  if (ok_dyn_ins) then
656 ! write(lunout,*) "leapfrog: call writehist, itau=",itau
657  CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
658 ! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
659 ! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
660 ! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
661 ! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
662 ! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
663  endif ! of if (ok_dyn_ins)
664 #endif
665 ! For some Grads outputs of fields
666  if (output_grads_dyn) then
667 #include "write_grads_dyn.h"
668  endif
669  endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
670  ENDIF ! of IF(MOD(itau,iecri).EQ.0)
671 
672  IF(itau.EQ.itaufin) THEN
673 
674 
675 ! if (planet_type.eq."earth") then
676 ! Write an Earth-format restart file
677  CALL dynredem1("restart.nc",start_time,
678  & vcov,ucov,teta,q,masse,ps)
679 ! endif ! of if (planet_type.eq."earth")
680 
681  CLOSE(99)
682  !!! Ehouarn: Why not stop here and now?
683  ENDIF ! of IF (itau.EQ.itaufin)
684 
685 c-----------------------------------------------------------------------
686 c gestion de l'integration temporelle:
687 c ------------------------------------
688 
689  IF( mod(itau,iperiod).EQ.0 ) THEN
690  go to 1
691  ELSE IF ( mod(itau-1,iperiod). eq. 0 ) THEN
692 
693  IF( forward ) THEN
694 c fin du pas forward et debut du pas backward
695 
696  forward = .false.
697  leapf = .false.
698  go to 2
699 
700  ELSE
701 c fin du pas backward et debut du premier pas leapfrog
702 
703  leapf = .true.
704  dt = 2.*dtvr
705  go to 2
706  END IF ! of IF (forward)
707  ELSE
708 
709 c ...... pas leapfrog .....
710 
711  leapf = .true.
712  dt = 2.*dtvr
713  go to 2
714  END IF ! of IF (MOD(itau,iperiod).EQ.0)
715  ! ELSEIF (MOD(itau-1,iperiod).EQ.0)
716 
717  ELSE ! of IF (.not.purmats)
718 
719 c ........................................................
720 c .............. schema matsuno ...............
721 c ........................................................
722  IF( forward ) THEN
723 
724  itau = itau + 1
725 c iday = day_ini+itau/day_step
726 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
727 c
728 c IF(time.GT.1.) THEN
729 c time = time-1.
730 c iday = iday+1
731 c ENDIF
732 
733  forward = .false.
734  IF( itau. eq. itaufinp1 ) then
735  abort_message = 'Simulation finished'
736  call abort_gcm(modname,abort_message,0)
737  ENDIF
738  go to 2
739 
740  ELSE ! of IF(forward) i.e. backward step
741 
742  IF(mod(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
743  IF(itau.EQ.itaufin) THEN
744  iav=1
745  ELSE
746  iav=0
747  ENDIF
748 
749  IF (ok_dynzon) THEN
750 #ifdef CPP_IOIPSL
751  CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
752  & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
753 #endif
754  ENDIF
755  IF (ok_dyn_ave) THEN
756 #ifdef CPP_IOIPSL
757  CALL writedynav(itau,vcov,
758  & ucov,teta,pk,phi,q,masse,ps,phis)
759 #endif
760  ENDIF
761 
762  ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
763 
764  IF(mod(itau,iecri ).EQ.0) THEN
765 c IF(MOD(itau,iecri*day_step).EQ.0) THEN
766  CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
767  unat=0.
768  do l=1,llm
769  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
770  vnat(:,l)=vcov(:,l)/cv(:)
771  enddo
772 #ifdef CPP_IOIPSL
773  if (ok_dyn_ins) then
774 ! write(lunout,*) "leapfrog: call writehist (b)",
775 ! & itau,iecri
776  CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
777  endif ! of if (ok_dyn_ins)
778 #endif
779 ! For some Grads outputs
780  if (output_grads_dyn) then
781 #include "write_grads_dyn.h"
782  endif
783 
784  ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
785 
786  IF(itau.EQ.itaufin) THEN
787 ! if (planet_type.eq."earth") then
788  CALL dynredem1("restart.nc",start_time,
789  & vcov,ucov,teta,q,masse,ps)
790 ! endif ! of if (planet_type.eq."earth")
791  ENDIF ! of IF(itau.EQ.itaufin)
792 
793  forward = .true.
794  go to 1
795 
796  ENDIF ! of IF (forward)
797 
798  END IF ! of IF(.not.purmats)
799 
800  stop
801  END