LMDZ
leapfrog.F
Go to the documentation of this file.
1 !
2 ! $Id: leapfrog.F 2375 2015-10-18 06:38:58Z emillour $
3 !
4 c
5 c
6  SUBROUTINE leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
7 
8 
9 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
10 #ifdef CPP_IOIPSL
11  use ioipsl
12 #endif
13  USE infotrac, ONLY: nqtot,ok_iso_verif
14  USE guide_mod, ONLY : guide_main
15  USE write_field, ONLY: writefield
16  USE control_mod, ONLY: nday, day_step, planet_type, offline,
17  & iconser, iphysiq, iperiod, dissip_period,
18  & iecri, ip_ebil_dyn, ok_dynzon, ok_dyn_ins,
19  & periodav, ok_dyn_ave, output_grads_dyn
20  use exner_hyb_m, only: exner_hyb
21  use exner_milieu_m, only: exner_milieu
22 
23  IMPLICIT NONE
24 
25 c ...... Version du 10/01/98 ..........
26 
27 c avec coordonnees verticales hybrides
28 c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
29 
30 c=======================================================================
31 c
32 c Auteur: P. Le Van /L. Fairhead/F.Hourdin
33 c -------
34 c
35 c Objet:
36 c ------
37 c
38 c GCM LMD nouvelle grille
39 c
40 c=======================================================================
41 c
42 c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
43 c et possibilite d'appeler une fonction f(y) a derivee tangente
44 c hyperbolique a la place de la fonction a derivee sinusoidale.
45 
46 c ... Possibilite de choisir le shema pour l'advection de
47 c q , en modifiant iadv dans traceur.def (10/02) .
48 c
49 c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
50 c Pour Van-Leer iadv=10
51 c
52 c-----------------------------------------------------------------------
53 c Declarations:
54 c -------------
55 
56 #include "dimensions.h"
57 #include "paramet.h"
58 #include "comconst.h"
59 #include "comdissnew.h"
60 #include "comvert.h"
61 #include "comgeom.h"
62 #include "logic.h"
63 #include "temps.h"
64 #include "ener.h"
65 #include "description.h"
66 #include "serre.h"
67 !#include "com_io_dyn.h"
68 #include "iniprint.h"
69 #include "academic.h"
70 
71  REAL,INTENT(IN) :: time_0 ! not used
72 
73 c dynamical variables:
74  REAL,INTENT(INOUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
75  REAL,INTENT(INOUT) :: vcov(ip1jm,llm) ! meridional covariant wind
76  REAL,INTENT(INOUT) :: teta(ip1jmp1,llm) ! potential temperature
77  REAL,INTENT(INOUT) :: ps(ip1jmp1) ! surface pressure (Pa)
78  REAL,INTENT(INOUT) :: masse(ip1jmp1,llm) ! air mass
79  REAL,INTENT(INOUT) :: phis(ip1jmp1) ! geopotentiat at the surface
80  REAL,INTENT(INOUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers
81 
82  REAL p (ip1jmp1,llmp1 ) ! interlayer pressure
83  REAL pks(ip1jmp1) ! exner at the surface
84  REAL pk(ip1jmp1,llm) ! exner at mid-layer
85  REAL pkf(ip1jmp1,llm) ! filtered exner at mid-layer
86  REAL phi(ip1jmp1,llm) ! geopotential
87  REAL w(ip1jmp1,llm) ! vertical velocity
88 
89  real zqmin,zqmax
90 
91 c variables dynamiques intermediaire pour le transport
92  REAL pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
93 
94 c variables dynamiques au pas -1
95  REAL vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
96  REAL tetam1(ip1jmp1,llm),psm1(ip1jmp1)
97  REAL massem1(ip1jmp1,llm)
98 
99 c tendances dynamiques
100  REAL dv(ip1jm,llm),du(ip1jmp1,llm)
101  REAL dteta(ip1jmp1,llm),dq(ip1jmp1,llm,nqtot),dp(ip1jmp1)
102 
103 c tendances de la dissipation
104  REAL dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
105  REAL dtetadis(ip1jmp1,llm)
106 
107 c tendances physiques
108  REAL dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
109  REAL dtetafi(ip1jmp1,llm),dqfi(ip1jmp1,llm,nqtot),dpfi(ip1jmp1)
110 
111 c variables pour le fichier histoire
112  REAL dtav ! intervalle de temps elementaire
113 
114  REAL tppn(iim),tpps(iim),tpn,tps
115 c
116  INTEGER itau,itaufinp1,iav
117 ! INTEGER iday ! jour julien
118  REAL time
119 
120  REAL SSUM
121 ! REAL finvmaold(ip1jmp1,llm)
122 
123 cym LOGICAL lafin
124  LOGICAL :: lafin=.false.
125  INTEGER ij,iq,l
126  INTEGER ik
127 
128  real time_step, t_wrt, t_ops
129 
130 ! REAL rdayvrai,rdaym_ini
131 ! jD_cur: jour julien courant
132 ! jH_cur: heure julienne courante
133  REAL :: jD_cur, jH_cur
134  INTEGER :: an, mois, jour
135  REAL :: secondes
136 
137  LOGICAL first,callinigrads
138 cIM : pour sortir les param. du modele dans un fis. netcdf 110106
139  save first
140  data first/.true./
141  real dt_cum
142  character*10 infile
143  integer zan, tau0, thoriid
144  integer nid_ctesGCM
145  save nid_ctesgcm
146  real degres
147  real rlong(iip1), rlatg(jjp1)
148  real zx_tmp_2d(iip1,jjp1)
149  integer ndex2d(iip1*jjp1)
150  logical ok_sync
151  parameter(ok_sync = .true.)
152  logical physic
153 
154  data callinigrads/.true./
155  character*10 string10
156 
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  if (nday>=0) then
195  itaufin = nday*day_step
196  else
197  itaufin = -nday
198  endif
199  itaufinp1 = itaufin +1
200  itau = 0
201  physic=.true.
202  if (iflag_phys==0.or.iflag_phys==2) physic=.false.
203 
204 c iday = day_ini+itau/day_step
205 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
206 c IF(time.GT.1.) THEN
207 c time = time-1.
208 c iday = iday+1
209 c ENDIF
210 
211 
212 c-----------------------------------------------------------------------
213 c On initialise la pression et la fonction d'Exner :
214 c --------------------------------------------------
215 
216  dq(:,:,:)=0.
217  CALL pression ( ip1jmp1, ap, bp, ps, p )
218  if (pressure_exner) then
219  CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
220  else
221  CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
222  endif
223 
224 c-----------------------------------------------------------------------
225 c Debut de l'integration temporelle:
226 c ----------------------------------
227 
228  1 CONTINUE ! Matsuno Forward step begins here
229 
230 c date: (NB: date remains unchanged for Backward step)
231 c -----
232 
233  jd_cur = jd_ref + day_ini - day_ref + &
234  & (itau+1)/day_step
235  jh_cur = jh_ref + start_time + &
236  & mod(itau+1,day_step)/float(day_step)
237  jd_cur = jd_cur + int(jh_cur)
238  jh_cur = jh_cur - int(jh_cur)
239 
240  if (ok_iso_verif) then
241  call check_isotopes_seq(q,ip1jmp1,'leapfrog 321')
242  endif !if (ok_iso_verif) then
243 
244 #ifdef CPP_IOIPSL
245  if (ok_guide) then
246  call guide_main(itau,ucov,vcov,teta,q,masse,ps)
247  endif
248 #endif
249 
250 
251 c
252 c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
253 c CALL test_period ( ucov,vcov,teta,q,p,phis )
254 c PRINT *,' ---- Test_period apres continue OK ! -----', itau
255 c ENDIF
256 c
257 
258 ! Save fields obtained at previous time step as '...m1'
259  CALL scopy( ijmllm ,vcov , 1, vcovm1 , 1 )
260  CALL scopy( ijp1llm,ucov , 1, ucovm1 , 1 )
261  CALL scopy( ijp1llm,teta , 1, tetam1 , 1 )
262  CALL scopy( ijp1llm,masse, 1, massem1, 1 )
263  CALL scopy( ip1jmp1, ps , 1, psm1 , 1 )
264 
265  forward = .true.
266  leapf = .false.
267  dt = dtvr
268 
269 c ... P.Le Van .26/04/94 ....
270 ! Ehouarn: finvmaold is actually not used
271 ! CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
272 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
273 
274  if (ok_iso_verif) then
275  call check_isotopes_seq(q,ip1jmp1,'leapfrog 400')
276  endif !if (ok_iso_verif) then
277 
278  2 CONTINUE ! Matsuno backward or leapfrog step begins here
279 
280 c-----------------------------------------------------------------------
281 
282 c date: (NB: only leapfrog step requires recomputing date)
283 c -----
284 
285  IF (leapf) THEN
286  jd_cur = jd_ref + day_ini - day_ref +
287  & (itau+1)/day_step
288  jh_cur = jh_ref + start_time +
289  & mod(itau+1,day_step)/float(day_step)
290  jd_cur = jd_cur + int(jh_cur)
291  jh_cur = jh_cur - int(jh_cur)
292  ENDIF
293 
294 
295 c gestion des appels de la physique et des dissipations:
296 c ------------------------------------------------------
297 c
298 c ... P.Le Van ( 6/02/95 ) ....
299 
300  apphys = .false.
301  statcl = .false.
302  conser = .false.
303  apdiss = .false.
304 
305  IF( purmats ) THEN
306  ! Purely Matsuno time stepping
307  IF( mod(itau,iconser) .EQ.0.AND. forward ) conser = .true.
308  IF( mod(itau,dissip_period ).EQ.0.AND..NOT.forward )
309  s apdiss = .true.
310  IF( mod(itau,iphysiq ).EQ.0.AND..NOT.forward
311  s .and. physic ) apphys = .true.
312  ELSE
313  ! Leapfrog/Matsuno time stepping
314  IF( mod(itau ,iconser) .EQ. 0 ) conser = .true.
315  IF( mod(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
316  s apdiss = .true.
317  IF( mod(itau+1,iphysiq).EQ.0.AND.physic ) apphys=.true.
318  END IF
319 
320 ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
321 ! supress dissipation step
322  if (llm.eq.1) then
323  apdiss=.false.
324  endif
325 
326 
327  if (ok_iso_verif) then
328  call check_isotopes_seq(q,ip1jmp1,'leapfrog 589')
329  endif !if (ok_iso_verif) then
330 
331 c-----------------------------------------------------------------------
332 c calcul des tendances dynamiques:
333 c --------------------------------
334 
335  ! compute geopotential phi()
336  CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
337 
338  time = jd_cur + jh_cur
339  CALL caldyn
340  $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
341  $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
342 
343 
344 c-----------------------------------------------------------------------
345 c calcul des tendances advection des traceurs (dont l'humidite)
346 c -------------------------------------------------------------
347 
348  if (ok_iso_verif) then
349  call check_isotopes_seq(q,ip1jmp1,
350  & 'leapfrog 686: avant caladvtrac')
351  endif !if (ok_iso_verif) then
352 
353  IF( forward. or . leapf ) THEN
354 ! Ehouarn: NB: fields sent to advtrac are those at the beginning of the time step
355  CALL caladvtrac(q,pbaru,pbarv,
356  * p, masse, dq, teta,
357  . flxw, pk)
358  !write(*,*) 'caladvtrac 346'
359 
360 
361  IF (offline) THEN
362 Cmaf stokage du flux de masse pour traceurs OFF-LINE
363 
364 #ifdef CPP_IOIPSL
365  CALL fluxstokenc(pbaru,pbarv,masse,teta,phi,phis,
366  . dtvr, itau)
367 #endif
368 
369 
370  ENDIF ! of IF (offline)
371 c
372  ENDIF ! of IF( forward. OR . leapf )
373 
374 
375 c-----------------------------------------------------------------------
376 c integrations dynamique et traceurs:
377 c ----------------------------------
378 
379  if (ok_iso_verif) then
380  write(*,*) 'leapfrog 720'
381  call check_isotopes_seq(q,ip1jmp1,'leapfrog 756')
382  endif !if (ok_iso_verif) then
383 
384  CALL integrd ( nqtot,vcovm1,ucovm1,tetam1,psm1,massem1 ,
385  $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
386 ! $ finvmaold )
387 
388  if (ok_iso_verif) then
389  write(*,*) 'leapfrog 724'
390  call check_isotopes_seq(q,ip1jmp1,'leapfrog 762')
391  endif !if (ok_iso_verif) then
392 
393 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
394 c
395 c-----------------------------------------------------------------------
396 c calcul des tendances physiques:
397 c -------------------------------
398 c ######## P.Le Van ( Modif le 6/02/95 ) ###########
399 c
400  IF( purmats ) THEN
401  IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .true.
402  ELSE
403  IF( itau+1. eq. itaufin ) lafin = .true.
404  ENDIF
405 c
406 c
407  IF( apphys ) THEN
408 c
409 c ....... Ajout P.Le Van ( 17/04/96 ) ...........
410 c
411 
412  CALL pression ( ip1jmp1, ap, bp, ps, p )
413  if (pressure_exner) then
414  CALL exner_hyb( ip1jmp1, ps, p,pks, pk, pkf )
415  else
416  CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
417  endif
418 
419 ! Appel a geopot ajoute le 2014/05/08 pour garantir la convergence numerique
420 ! avec dyn3dmem
421  CALL geopot ( ip1jmp1, teta , pk , pks, phis , phi )
422 
423 ! rdaym_ini = itau * dtvr / daysec
424 ! rdayvrai = rdaym_ini + day_ini
425 ! jD_cur = jD_ref + day_ini - day_ref
426 ! $ + int (itau * dtvr / daysec)
427 ! jH_cur = jH_ref + &
428 ! & (itau * dtvr / daysec - int(itau * dtvr / daysec))
429  jd_cur = jd_ref + day_ini - day_ref + &
430  & (itau+1)/day_step
431 
432  IF (planet_type .eq."generic") THEN
433  ! AS: we make jD_cur to be pday
434  jd_cur = int(day_ini + itau/day_step)
435  ENDIF
436 
437  jh_cur = jh_ref + start_time + &
438  & mod(itau+1,day_step)/float(day_step)
439  jd_cur = jd_cur + int(jh_cur)
440  jh_cur = jh_cur - int(jh_cur)
441 ! write(lunout,*)'itau, jD_cur = ', itau, jD_cur, jH_cur
442 ! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
443 ! write(lunout,*)'current date = ',an, mois, jour, secondes
444 
445 c rajout debug
446 c lafin = .true.
447 
448 
449 c Inbterface avec les routines de phylmd (phymars ... )
450 c -----------------------------------------------------
451 
452 c+jld
453 
454 c Diagnostique de conservation de l'énergie : initialisation
455  IF (ip_ebil_dyn.ge.1 ) THEN
456  ztit='bil dyn'
457 ! Ehouarn: be careful, diagedyn is Earth-specific!
458  IF (planet_type.eq."earth") THEN
459  CALL diagedyn(ztit,2,1,1,dtphys
460  & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
461  ENDIF
462  ENDIF ! of IF (ip_ebil_dyn.ge.1 )
463 c-jld
464 #ifdef CPP_IOIPSL
465 cIM decommenter les 6 lignes suivantes pour sortir quelques parametres dynamiques de LMDZ
466 cIM uncomment next 6 lines to get some parameters for LMDZ dynamics
467 c IF (first) THEN
468 c first=.false.
469 c#include "ini_paramLMDZ_dyn.h"
470 c ENDIF
471 c
472 c#include "write_paramLMDZ_dyn.h"
473 c
474 #endif
475 ! #endif of #ifdef CPP_IOIPSL
476 #ifdef CPP_PHYS
477  CALL calfis( lafin , jd_cur, jh_cur,
478  $ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
479  $ du,dv,dteta,dq,
480  $ flxw,dufi,dvfi,dtetafi,dqfi,dpfi )
481 #endif
482 c ajout des tendances physiques:
483 c ------------------------------
484  CALL addfi( dtphys, leapf, forward ,
485  $ ucov, vcov, teta , q ,ps ,
486  $ dufi, dvfi, dtetafi , dqfi ,dpfi )
487  ! since addfi updates ps(), also update p(), masse() and pk()
488  CALL pression (ip1jmp1,ap,bp,ps,p)
489  CALL massdair(p,masse)
490  if (pressure_exner) then
491  CALL exner_hyb(ip1jmp1,ps,p,pks,pk,pkf)
492  else
493  CALL exner_milieu(ip1jmp1,ps,p,pks,pk,pkf)
494  endif
495 
496  IF (ok_strato) THEN
497  CALL top_bound( vcov,ucov,teta,masse,dtphys)
498  ENDIF
499 
500 c
501 c Diagnostique de conservation de l'énergie : difference
502  IF (ip_ebil_dyn.ge.1 ) THEN
503  ztit='bil phys'
504  IF (planet_type.eq."earth") THEN
505  CALL diagedyn(ztit,2,1,1,dtphys
506  & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
507  ENDIF
508  ENDIF ! of IF (ip_ebil_dyn.ge.1 )
509 
510  ENDIF ! of IF( apphys )
511 
512  IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
513 ! Academic case : Simple friction and Newtonan relaxation
514 ! -------------------------------------------------------
515  DO l=1,llm
516  DO ij=1,ip1jmp1
517  teta(ij,l)=teta(ij,l)-dtvr*
518  & (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
519  ENDDO
520  ENDDO ! of DO l=1,llm
521 
522  if (planet_type.eq."giant") then
523  ! add an intrinsic heat flux at the base of the atmosphere
524  teta(:,1)=teta(:,1)+dtvr*aire(:)*ihf/cpp/masse(:,1)
525  endif
526 
527  call friction(ucov,vcov,dtvr)
528 
529  ! Sponge layer (if any)
530  IF (ok_strato) THEN
531 ! dufi(:,:)=0.
532 ! dvfi(:,:)=0.
533 ! dtetafi(:,:)=0.
534 ! dqfi(:,:,:)=0.
535 ! dpfi(:)=0.
536 ! CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
537  CALL top_bound( vcov,ucov,teta,masse,dtvr)
538 ! CALL addfi( dtvr, leapf, forward ,
539 ! $ ucov, vcov, teta , q ,ps ,
540 ! $ dufi, dvfi, dtetafi , dqfi ,dpfi )
541  ENDIF ! of IF (ok_strato)
542  ENDIF ! of IF (iflag_phys.EQ.2)
543 
544 
545 c-jld
546 
547  CALL pression ( ip1jmp1, ap, bp, ps, p )
548  if (pressure_exner) then
549  CALL exner_hyb( ip1jmp1, ps, p, pks, pk, pkf )
550  else
551  CALL exner_milieu( ip1jmp1, ps, p, pks, pk, pkf )
552  endif
553  CALL massdair(p,masse)
554 
555  if (ok_iso_verif) then
556  call check_isotopes_seq(q,ip1jmp1,'leapfrog 1196')
557  endif !if (ok_iso_verif) then
558 
559 c-----------------------------------------------------------------------
560 c dissipation horizontale et verticale des petites echelles:
561 c ----------------------------------------------------------
562 
563  IF(apdiss) THEN
564 
565 
566 c calcul de l'energie cinetique avant dissipation
567  call covcont(llm,ucov,vcov,ucont,vcont)
568  call enercin(vcov,ucov,vcont,ucont,ecin0)
569 
570 c dissipation
571  CALL dissip(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
572  ucov=ucov+dudis
573  vcov=vcov+dvdis
574 c teta=teta+dtetadis
575 
576 
577 c------------------------------------------------------------------------
578  if (dissip_conservative) then
579 C On rajoute la tendance due a la transform. Ec -> E therm. cree
580 C lors de la dissipation
581  call covcont(llm,ucov,vcov,ucont,vcont)
582  call enercin(vcov,ucov,vcont,ucont,ecin)
583  dtetaecdt= (ecin0-ecin)/ pk
584 c teta=teta+dtetaecdt
585  dtetadis=dtetadis+dtetaecdt
586  endif
587  teta=teta+dtetadis
588 c------------------------------------------------------------------------
589 
590 
591 c ....... P. Le Van ( ajout le 17/04/96 ) ...........
592 c ... Calcul de la valeur moyenne, unique de h aux poles .....
593 c
594 
595  DO l = 1, llm
596  DO ij = 1,iim
597  tppn(ij) = aire( ij ) * teta( ij ,l)
598  tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
599  ENDDO
600  tpn = ssum(iim,tppn,1)/apoln
601  tps = ssum(iim,tpps,1)/apols
602 
603  DO ij = 1, iip1
604  teta( ij ,l) = tpn
605  teta(ij+ip1jm,l) = tps
606  ENDDO
607  ENDDO
608 
609  if (1 == 0) then
610 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
611 !!! 2) should probably not be here anyway
612 !!! but are kept for those who would want to revert to previous behaviour
613  DO ij = 1,iim
614  tppn(ij) = aire( ij ) * ps( ij )
615  tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
616  ENDDO
617  tpn = ssum(iim,tppn,1)/apoln
618  tps = ssum(iim,tpps,1)/apols
619 
620  DO ij = 1, iip1
621  ps( ij ) = tpn
622  ps(ij+ip1jm) = tps
623  ENDDO
624  endif ! of if (1 == 0)
625 
626  END IF ! of IF(apdiss)
627 
628 c ajout debug
629 c IF( lafin ) then
630 c abort_message = 'Simulation finished'
631 c call abort_gcm(modname,abort_message,0)
632 c ENDIF
633 
634 c ********************************************************************
635 c ********************************************************************
636 c .... fin de l'integration dynamique et physique pour le pas itau ..
637 c ********************************************************************
638 c ********************************************************************
639 
640 c preparation du pas d'integration suivant ......
641 
642  if (ok_iso_verif) then
643  call check_isotopes_seq(q,ip1jmp1,'leapfrog 1509')
644  endif !if (ok_iso_verif) then
645 
646  IF ( .NOT.purmats ) THEN
647 c ........................................................
648 c .............. schema matsuno + leapfrog ..............
649 c ........................................................
650 
651  IF(forward. or. leapf) THEN
652  itau= itau + 1
653 c iday= day_ini+itau/day_step
654 c time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
655 c IF(time.GT.1.) THEN
656 c time = time-1.
657 c iday = iday+1
658 c ENDIF
659  ENDIF
660 
661 
662  IF( itau. eq. itaufinp1 ) then
663  if (flag_verif) then
664  write(79,*) 'ucov',ucov
665  write(80,*) 'vcov',vcov
666  write(81,*) 'teta',teta
667  write(82,*) 'ps',ps
668  write(83,*) 'q',q
669  WRITE(85,*) 'q1 = ',q(:,:,1)
670  WRITE(86,*) 'q3 = ',q(:,:,3)
671  endif
672 
673  abort_message = 'Simulation finished'
674 
675  call abort_gcm(modname,abort_message,0)
676  ENDIF
677 c-----------------------------------------------------------------------
678 c ecriture du fichier histoire moyenne:
679 c -------------------------------------
680 
681  IF(mod(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
682  IF(itau.EQ.itaufin) THEN
683  iav=1
684  ELSE
685  iav=0
686  ENDIF
687 
688  IF (ok_dynzon) THEN
689 #ifdef CPP_IOIPSL
690  CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
691  & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
692 #endif
693  END IF
694  IF (ok_dyn_ave) THEN
695 #ifdef CPP_IOIPSL
696  CALL writedynav(itau,vcov,
697  & ucov,teta,pk,phi,q,masse,ps,phis)
698 #endif
699  ENDIF
700 
701  ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
702 
703  if (ok_iso_verif) then
704  call check_isotopes_seq(q,ip1jmp1,'leapfrog 1584')
705  endif !if (ok_iso_verif) then
706 
707 c-----------------------------------------------------------------------
708 c ecriture de la bande histoire:
709 c ------------------------------
710 
711  IF( mod(itau,iecri).EQ.0) THEN
712  ! Ehouarn: output only during LF or Backward Matsuno
713  if (leapf.or.(.not.leapf.and.(.not.forward))) then
714  CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
715  unat=0.
716  do l=1,llm
717  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
718  vnat(:,l)=vcov(:,l)/cv(:)
719  enddo
720 #ifdef CPP_IOIPSL
721  if (ok_dyn_ins) then
722 ! write(lunout,*) "leapfrog: call writehist, itau=",itau
723  CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
724 ! call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
725 ! call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
726 ! call WriteField('teta',reshape(teta,(/iip1,jmp1,llm/)))
727 ! call WriteField('ps',reshape(ps,(/iip1,jmp1/)))
728 ! call WriteField('masse',reshape(masse,(/iip1,jmp1,llm/)))
729  endif ! of if (ok_dyn_ins)
730 #endif
731 ! For some Grads outputs of fields
732  if (output_grads_dyn) then
733 #include "write_grads_dyn.h"
734  endif
735  endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
736  ENDIF ! of IF(MOD(itau,iecri).EQ.0)
737 
738  IF(itau.EQ.itaufin) THEN
739 
740 
741 ! if (planet_type.eq."earth") then
742 ! Write an Earth-format restart file
743  CALL dynredem1("restart.nc",start_time,
744  & vcov,ucov,teta,q,masse,ps)
745 ! endif ! of if (planet_type.eq."earth")
746 
747  CLOSE(99)
748  !!! Ehouarn: Why not stop here and now?
749  ENDIF ! of IF (itau.EQ.itaufin)
750 
751 c-----------------------------------------------------------------------
752 c gestion de l'integration temporelle:
753 c ------------------------------------
754 
755  IF( mod(itau,iperiod).EQ.0 ) THEN
756  GO TO 1
757  ELSE IF ( mod(itau-1,iperiod). eq. 0 ) THEN
758 
759  IF( forward ) THEN
760 c fin du pas forward et debut du pas backward
761 
762  forward = .false.
763  leapf = .false.
764  GO TO 2
765 
766  ELSE
767 c fin du pas backward et debut du premier pas leapfrog
768 
769  leapf = .true.
770  dt = 2.*dtvr
771  GO TO 2
772  END IF ! of IF (forward)
773  ELSE
774 
775 c ...... pas leapfrog .....
776 
777  leapf = .true.
778  dt = 2.*dtvr
779  GO TO 2
780  END IF ! of IF (MOD(itau,iperiod).EQ.0)
781  ! ELSEIF (MOD(itau-1,iperiod).EQ.0)
782 
783  ELSE ! of IF (.not.purmats)
784 
785  if (ok_iso_verif) then
786  call check_isotopes_seq(q,ip1jmp1,'leapfrog 1664')
787  endif !if (ok_iso_verif) then
788 
789 c ........................................................
790 c .............. schema matsuno ...............
791 c ........................................................
792  IF( forward ) THEN
793 
794  itau = itau + 1
795 c iday = day_ini+itau/day_step
796 c time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
797 c
798 c IF(time.GT.1.) THEN
799 c time = time-1.
800 c iday = iday+1
801 c ENDIF
802 
803  forward = .false.
804  IF( itau. eq. itaufinp1 ) then
805  abort_message = 'Simulation finished'
806  call abort_gcm(modname,abort_message,0)
807  ENDIF
808  GO TO 2
809 
810  ELSE ! of IF(forward) i.e. backward step
811 
812  if (ok_iso_verif) then
813  call check_isotopes_seq(q,ip1jmp1,'leapfrog 1698')
814  endif !if (ok_iso_verif) then
815 
816  IF(mod(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
817  IF(itau.EQ.itaufin) THEN
818  iav=1
819  ELSE
820  iav=0
821  ENDIF
822 
823  IF (ok_dynzon) THEN
824 #ifdef CPP_IOIPSL
825  CALL bilan_dyn(2,dtvr*iperiod,dtvr*day_step*periodav,
826  & ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
827 #endif
828  ENDIF
829  IF (ok_dyn_ave) THEN
830 #ifdef CPP_IOIPSL
831  CALL writedynav(itau,vcov,
832  & ucov,teta,pk,phi,q,masse,ps,phis)
833 #endif
834  ENDIF
835 
836  ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
837 
838  IF(mod(itau,iecri ).EQ.0) THEN
839 c IF(MOD(itau,iecri*day_step).EQ.0) THEN
840  CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
841  unat=0.
842  do l=1,llm
843  unat(iip2:ip1jm,l)=ucov(iip2:ip1jm,l)/cu(iip2:ip1jm)
844  vnat(:,l)=vcov(:,l)/cv(:)
845  enddo
846 #ifdef CPP_IOIPSL
847  if (ok_dyn_ins) then
848 ! write(lunout,*) "leapfrog: call writehist (b)",
849 ! & itau,iecri
850  CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
851  endif ! of if (ok_dyn_ins)
852 #endif
853 ! For some Grads outputs
854  if (output_grads_dyn) then
855 #include "write_grads_dyn.h"
856  endif
857 
858  ENDIF ! of IF(MOD(itau,iecri ).EQ.0)
859 
860  IF(itau.EQ.itaufin) THEN
861 ! if (planet_type.eq."earth") then
862  CALL dynredem1("restart.nc",start_time,
863  & vcov,ucov,teta,q,masse,ps)
864 ! endif ! of if (planet_type.eq."earth")
865  ENDIF ! of IF(itau.EQ.itaufin)
866 
867  forward = .true.
868  GO TO 1
869 
870  ENDIF ! of IF (forward)
871 
872  END IF ! of IF(.not.purmats)
873 
874  stop
875  END
fluxstokenc
subroutine fluxstokenc(pbaru, pbarv, masse, teta, phi, phis, time_step, itau)
Definition: fluxstokenc.F:6
diagedyn
subroutine diagedyn(tit, iprt, idiag, idiag2, dtime, ucov, vcov, ps, p, pk, teta, q, ql)
Definition: diagedyn.F:8
iip2
!$Header iip2
Definition: paramet.h:14
ip1jmp1
!$Header llmm1 INTEGER ip1jmp1
Definition: paramet.h:14
calfis
subroutine calfis(lafin, jD_cur, jH_cur, pucov, pvcov, pteta, pq, pmasse, pps, pp, ppk, pphis, pphi, pducov, pdvcov, pdteta, pdq, flxw, pdufi, pdvfi, pdhfi, pdqfi, pdpsfi)
Definition: calfis.F:28
apols
!$Header!CDK comgeom COMMON comgeom apols
Definition: comgeom.h:8
exner_milieu_m
Definition: exner_milieu_m.F90:1
write_field::writefield
Definition: write_field.F90:15
covcont
subroutine covcont(klevel, ucov, vcov, ucont, vcont)
Definition: covcont.F90:2
integrd
subroutine integrd(nq, vcovm1, ucovm1, tetam1, psm1, massem1, dv, du, dteta, dq, dp, vcov, ucov, teta, q, ps, masse, phis)
Definition: integrd.F:8
guide_mod::guide_main
subroutine guide_main(itau, ucov, vcov, teta, q, masse, ps)
Definition: guide_mod.F90:315
llmp1
!$Header llmp1
Definition: paramet.h:14
bp
!$Id bp(llm+1)
exner_hyb_m::exner_hyb
subroutine exner_hyb(ngrid, ps, p, pks, pk, pkf)
Definition: exner_hyb_m.F90:8
infotrac::ok_iso_verif
logical, save ok_iso_verif
Definition: infotrac.F90:44
caldyn
subroutine caldyn(itau, ucov, vcov, teta, ps, masse, pk, pkf, phis, phi, conser, du, dv, dteta, dp, w, pbaru, pbarv, time)
Definition: caldyn.F:7
addfi
subroutine addfi(pdt, leapf, forward, pucov, pvcov, pteta, pq, pps, pdufi, pdvfi, pdhfi, pdqfi, pdpfi)
Definition: addfi.F:7
knewt_t
!$Id knewt_t
Definition: academic.h:4
dt
!$Id calend INTEGER itaufin INTEGER itau_phy INTEGER day_ref REAL dt
Definition: temps.h:15
exner_milieu_m::exner_milieu
subroutine exner_milieu(ngrid, ps, p, pks, pk, pkf)
Definition: exner_milieu_m.F90:8
control_mod::planet_type
character(len=10), save planet_type
Definition: control_mod.F90:32
abort_gcm
subroutine abort_gcm(modname, message, ierr)
Definition: abort_gcm.F:7
llm
!$Id Turb_fcg_gcssold get_uvd hqturb_gcssold endif!large scale llm day day1 day day1 *dt_toga endif!time annee_ref dt_toga u_toga vq_toga w_prof vq_prof llm day day1 day day1 *dt_dice endif!time annee_ref dt_dice swup_dice vg_dice omega_dice tg_prof vg_profd w_profd omega_profd!do llm!print llm l llm
Definition: 1D_interp_cases.h:103
aire
!$Header!CDK comgeom COMMON comgeom aire
Definition: comgeom.h:25
iflag_phys
!$Id ysinus ok_gradsfile hybrid COMMON logici iflag_phys
Definition: logic.h:10
control_mod::day_step
integer, save day_step
Definition: control_mod.F90:15
geopot
subroutine geopot(ngrid, teta, pk, pks, phis, phi)
Definition: geopot.F:5
knewt_g
!$Id knewt_g
Definition: academic.h:4
infotrac::nqtot
integer, save nqtot
Definition: infotrac.F90:6
scopy
subroutine scopy(n, sx, incx, sy, incy)
Definition: cray.F:9
ok_strato
!$Id ysinus ok_strato
Definition: logic.h:10
writehist
subroutine writehist(time, vcov, ucov, teta, phi, q, masse, ps, phis)
Definition: writehist.F:5
statcl
!$Id && statcl
Definition: logic.h:10
ip1jm
!$Header llmm1 INTEGER ip1jm
Definition: paramet.h:14
pression
subroutine pression(ngrid, ap, bp, ps, p)
Definition: pression.F90:2
conser
!$Id conser
Definition: logic.h:10
day_ini
!$Id && day_ini
Definition: temps.h:15
dynredem1
subroutine dynredem1(fichnom, time, vcov, ucov, teta, q, masse, ps)
Definition: dynredem.F90:160
dtphys
!$Id mode_top_bound COMMON comconstr dtphys
Definition: comconst.h:7
apoln
!$Header!CDK comgeom COMMON comgeom apoln
Definition: comgeom.h:8
false
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL false
Definition: calcul_STDlev.h:26
control_mod::offline
logical, save offline
Definition: control_mod.F90:30
ijmllm
!$Header llmm1 INTEGER ijp1llm INTEGER ijmllm
Definition: paramet.h:14
dissip
subroutine dissip(vcov, ucov, teta, p, dv, du, dh)
Definition: dissip.F:5
ok_guide
!$Id ysinus ok_guide
Definition: logic.h:10
day_ref
!$Id day_ref
Definition: temps.h:15
jjp1
!$Header jjp1
Definition: paramet.h:14
caladvtrac
subroutine caladvtrac(q, pbaru, pbarv, p, masse, dq, teta, flxw, pk)
Definition: caladvtrac.F:9
apdiss
!$Id apdiss
Definition: logic.h:10
cpp
!$Id mode_top_bound COMMON comconstr cpp
Definition: comconst.h:7
writedynav
subroutine writedynav(time, vcov, ucov, teta, ppk, phi, q, masse, ps, phis)
Definition: writedynav.F90:4
parameter
!$Header!integer nvarmx parameter(nfmx=10, imx=200, jmx=150, lmx=200, nvarmx=1000) real xd(imx
apphys
!$Id apphys
Definition: logic.h:10
exner_hyb_m
Definition: exner_hyb_m.F90:1
write_field
Definition: write_field.F90:4
true
!$Id itapm1 ENDIF!IM on interpole les champs sur les niveaux STD de pression!IM a chaque pas de temps de la physique c!positionnement de l argument logique a false c!pour ne pas recalculer deux fois la meme chose!c!a cet effet un appel a plevel_new a ete deplace c!a la fin de la serie d appels c!la boucle DO nlevSTD a ete internalisee c!dans d ou la creation de cette routine c c!CALL ulevSTD CALL &zphi philevSTD CALL &zx_rh rhlevSTD!DO klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon klev DO klon du jour ou toutes les read_climoz CALL true
Definition: calcul_STDlev.h:139
guide_mod
Definition: guide_mod.F90:4
cu
!$Id!Parameters for parameters that control the rate of approach!to quasi equilibrium noff nlm real tlcrit real entp real sigd real coeffs real dtmax real cu real betad real damp real delta COMMON cvparam nlm tlcrit sigd coeffs cu
Definition: cvparam.h:12
purmats
!$Id ysinus ok_gradsfile hybrid COMMON logici iflag_trac LOGICAL purmats
Definition: logic.h:10
dtvr
!$Id mode_top_bound COMMON comconstr dtvr
Definition: comconst.h:7
check_isotopes_seq
subroutine check_isotopes_seq(q, ip1jmp1, err_msg)
Definition: check_isotopes.F:2
forward
!$Id forward
Definition: logic.h:10
leapf
!$Id leapf
Definition: logic.h:10
leapfrog
subroutine leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
Definition: leapfrog.F:7
iim
c c zjulian c cym CALL iim cym klev iim
Definition: ini_bilKP_ave.h:24
friction
subroutine friction(ucov, vcov, pdt)
Definition: friction.F:6
itaufin
!$Id itaufin
Definition: temps.h:15
massdair
subroutine massdair(p, masse)
Definition: massdair.F:5
control_mod::nday
integer, save nday
Definition: control_mod.F90:14
cv
!$Header!CDK comgeom COMMON comgeom cv
Definition: comgeom.h:25
enercin
subroutine enercin(vcov, ucov, vcont, ucont, ecin)
Definition: enercin.F90:2
top_bound
subroutine top_bound(vcov, ucov, teta, masse, dt)
Definition: top_bound.F:5
bilan_dyn
subroutine bilan_dyn(ntrac, dt_app, dt_cum, ps, masse, pk, flux_u, flux_v, teta, phi, ucov, vcov, trac)
Definition: bilan_dyn.F:6
start_time
!$Id start_time
Definition: temps.h:15
control_mod
Definition: control_mod.F90:5