My Project
 All Classes Files Functions Variables Macros
leapfrog_p.F
Go to the documentation of this file.
1 !
2 ! $Id: leapfrog_p.F 1676 2012-11-08 20:03:30Z aslmd $
3 !
4 c
5 c
6 
7  SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
8  & time_0)
9 
10  USE misc_mod
11  USE parallel
12  USE times
13  USE mod_hallo
14  USE bands
15  USE write_field
16  USE write_field_p
17  USE vampir
18  USE timer_filtre, ONLY : print_filtre_timer
19  USE infotrac
20  USE guide_p_mod, ONLY : guide_main
21  USE getparam
22  USE control_mod
23 
24  IMPLICIT NONE
25 
26 c ...... Version du 10/01/98 ..........
27 
28 c avec coordonnees verticales hybrides
29 c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
30 
31 c=======================================================================
32 c
33 c Auteur: P. Le Van /L. Fairhead/F.Hourdin
34 c -------
35 c
36 c Objet:
37 c ------
38 c
39 c GCM LMD nouvelle grille
40 c
41 c=======================================================================
42 c
43 c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
44 c et possibilite d'appeler une fonction f(y) a derivee tangente
45 c hyperbolique a la place de la fonction a derivee sinusoidale.
46 
47 c ... Possibilite de choisir le shema pour l'advection de
48 c q , en modifiant iadv dans traceur.def (10/02) .
49 c
50 c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
51 c Pour Van-Leer iadv=10
52 c
53 c-----------------------------------------------------------------------
54 c Declarations:
55 c -------------
56 
57 #include "dimensions.h"
58 #include "paramet.h"
59 #include "comconst.h"
60 #include "comdissnew.h"
61 #include "comvert.h"
62 #include "comgeom.h"
63 #include "logic.h"
64 #include "temps.h"
65 #include "ener.h"
66 #include "description.h"
67 #include "serre.h"
68 !#include "com_io_dyn.h"
69 #include "iniprint.h"
70 #include "academic.h"
71 
72  INTEGER longcles
73  parameter( longcles = 20 )
74  REAL clesphy0( longcles )
75 
76  real zqmin,zqmax
77 
78 c variables dynamiques
79  REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
80  REAL :: teta(ip1jmp1,llm) ! temperature potentielle
81  REAL :: q(ip1jmp1,llm,nqtot) ! champs advectes
82  REAL :: ps(ip1jmp1) ! pression au sol
83  REAL,SAVE :: p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
84  REAL,SAVE :: pks(ip1jmp1) ! exner au sol
85  REAL,SAVE :: pk(ip1jmp1,llm) ! exner au milieu des couches
86  REAL,SAVE :: pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
87  REAL :: masse(ip1jmp1,llm) ! masse d'air
88  REAL :: phis(ip1jmp1) ! geopotentiel au sol
89  REAL,SAVE :: phi(ip1jmp1,llm) ! geopotentiel
90  REAL,SAVE :: w(ip1jmp1,llm) ! vitesse verticale
91 
92 c variables dynamiques intermediaire pour le transport
93  REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
94 
95 c variables dynamiques au pas -1
96  REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
97  REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
98  REAL,SAVE :: massem1(ip1jmp1,llm)
99 
100 c tendances dynamiques
101  REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
102  REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
103  REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
104 
105 c tendances de la dissipation
106  REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
107  REAL,SAVE :: dtetadis(ip1jmp1,llm)
108 
109 c tendances physiques
110  REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
111  REAL,SAVE :: dtetafi(ip1jmp1,llm)
112  REAL,SAVE :: dpfi(ip1jmp1)
113  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
114 
115 c variables pour le fichier histoire
116  REAL dtav ! intervalle de temps elementaire
117 
118  REAL tppn(iim),tpps(iim),tpn,tps
119 c
120  INTEGER itau,itaufinp1,iav
121 ! INTEGER iday ! jour julien
122  REAL time
123 
124  REAL ssum
125  REAL time_0
126 ! REAL,SAVE :: finvmaold(ip1jmp1,llm)
127 
128 cym LOGICAL lafin
129  LOGICAL :: lafin
130  INTEGER ij,iq,l
131  INTEGER ik
132 
133  real time_step, t_wrt, t_ops
134 
135 ! jD_cur: jour julien courant
136 ! jH_cur: heure julienne courante
137  REAL :: jd_cur, jh_cur
138  INTEGER :: an, mois, jour
139  REAL :: secondes
140 
141  logical :: physic
142  LOGICAL first,callinigrads
143 
144  data callinigrads/.true./
145  character*10 string10
146 
147  REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
148  REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
149 
150 c+jld variables test conservation energie
151  REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
152 C Tendance de la temp. potentiel d (theta)/ d t due a la
153 C tansformation d'energie cinetique en energie thermique
154 C cree par la dissipation
155  REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
156  REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
157  REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
158  REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
159  CHARACTER*15 ztit
160 ! INTEGER ip_ebil_dyn ! PRINT level for energy conserv. diag.
161 ! SAVE ip_ebil_dyn
162 ! DATA ip_ebil_dyn/0/
163 c-jld
164 
165  character*80 dynhist_file, dynhistave_file
166  character(len=*),parameter :: modname="leapfrog"
167  character*80 abort_message
168 
169 
170  logical,PARAMETER :: dissip_conservative=.true.
171 
172  INTEGER testita
173  parameter(testita = 9)
174 
175  logical , parameter :: flag_verif = .false.
176 
177 c declaration liees au parallelisme
178  INTEGER :: ierr
179  LOGICAL :: firstcaldyn
180  LOGICAL :: firstphysic
181  INTEGER :: ijb,ije,j,i
182  type(request) :: testrequest
183  type(request) :: request_dissip
184  type(request) :: request_physic
185  REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
186  REAL,SAVE :: dtetafi_tmp(iip1,llm)
187  REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
188  REAL,SAVE :: dpfi_tmp(iip1)
189 
190  INTEGER :: true_itau
191  INTEGER :: iapptrac
192  INTEGER :: adjustcount
193 ! INTEGER :: var_time
194  LOGICAL :: ok_start_timer=.false.
195  LOGICAL, SAVE :: firstcall=.true.
196 
197 c$OMP MASTER
198  itcount=0
199 c$OMP END MASTER
200  true_itau=0
201  firstcaldyn=.true.
202  firstphysic=.true.
203  iapptrac=0
204  adjustcount = 0
205  lafin=.false.
206 
207  itaufin = nday*day_step
208  itaufinp1 = itaufin +1
209 
210  itau = 0
211  physic=.true.
212  if (iflag_phys==0.or.iflag_phys==2) physic=.false.
213 ! iday = day_ini+itau/day_step
214 ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
215 ! IF(time.GT.1.) THEN
216 ! time = time-1.
217 ! iday = iday+1
218 ! ENDIF
219 
220 c Allocate variables depending on dynamic variable nqtot
221 c$OMP MASTER
222  IF (firstcall) THEN
223  firstcall=.false.
224  ALLOCATE(dq(ip1jmp1,llm,nqtot))
225  ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
226  ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
227  END IF
228 c$OMP END MASTER
229 c$OMP BARRIER
230 
231 c-----------------------------------------------------------------------
232 c On initialise la pression et la fonction d'Exner :
233 c --------------------------------------------------
234 
235 c$OMP MASTER
236  dq(:,:,:)=0.
237  CALL pression( ip1jmp1, ap, bp, ps, p )
238  if (pressure_exner) then
239  CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
240  else
241  CALL exner_milieu( ip1jmp1, ps, p, beta, pks, pk, pkf )
242  endif
243 c$OMP END MASTER
244 c-----------------------------------------------------------------------
245 c Debut de l'integration temporelle:
246 c ----------------------------------
247 c et du parallelisme !!
248 
249  1 CONTINUE ! Matsuno Forward step begins here
250 
251  jd_cur = jd_ref + day_ini - day_ref + &
252  & itau/day_step
253  jh_cur = jh_ref + start_time + &
254  & mod(itau,day_step)/float(day_step)
255  if (jh_cur > 1.0 ) then
256  jd_cur = jd_cur +1.
257  jh_cur = jh_cur -1.
258  endif
259 
260 
261 #ifdef CPP_IOIPSL
262  if (ok_guide) then
263 !$OMP MASTER
264  call guide_main(itau,ucov,vcov,teta,q,masse,ps)
265 !$OMP END MASTER
266 !$OMP BARRIER
267  endif
268 #endif
269 
270 c
271 c IF( MOD( itau, 10* day_step ).EQ.0 ) THEN
272 c CALL test_period ( ucov,vcov,teta,q,p,phis )
273 c PRINT *,' ---- Test_period apres continue OK ! -----', itau
274 c ENDIF
275 c
276 cym CALL SCOPY( ijmllm ,vcov , 1, vcovm1 , 1 )
277 cym CALL SCOPY( ijp1llm,ucov , 1, ucovm1 , 1 )
278 cym CALL SCOPY( ijp1llm,teta , 1, tetam1 , 1 )
279 cym CALL SCOPY( ijp1llm,masse, 1, massem1, 1 )
280 cym CALL SCOPY( ip1jmp1, ps , 1, psm1 , 1 )
281 
282  if (firstcaldyn) then
283 c$OMP MASTER
284  ucovm1=ucov
285  vcovm1=vcov
286  tetam1= teta
287  massem1= masse
288  psm1= ps
289 
290 ! Ehouarn: finvmaold is actually not used
291 ! finvmaold = masse
292 ! CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
293 c$OMP END MASTER
294 c$OMP BARRIER
295  else
296 ! Save fields obtained at previous time step as '...m1'
297  ijb=ij_begin
298  ije=ij_end
299 
300 c$OMP MASTER
301  psm1(ijb:ije) = ps(ijb:ije)
302 c$OMP END MASTER
303 
304 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
305  DO l=1,llm
306  ije=ij_end
307  ucovm1(ijb:ije,l) = ucov(ijb:ije,l)
308  tetam1(ijb:ije,l) = teta(ijb:ije,l)
309  massem1(ijb:ije,l) = masse(ijb:ije,l)
310 ! finvmaold(ijb:ije,l)=masse(ijb:ije,l)
311 
312  if (pole_sud) ije=ij_end-iip1
313  vcovm1(ijb:ije,l) = vcov(ijb:ije,l)
314 
315 
316  ENDDO
317 c$OMP ENDDO
318 
319 ! Ehouarn: finvmaold not used
320 ! CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
321 ! . llm, -2,2, .TRUE., 1 )
322 
323  endif ! of if (FirstCaldyn)
324 
325  forward = .true.
326  leapf = .false.
327  dt = dtvr
328 
329 c ... P.Le Van .26/04/94 ....
330 
331 cym CALL SCOPY ( ijp1llm, masse, 1, finvmaold, 1 )
332 cym CALL filtreg ( finvmaold ,jjp1, llm, -2,2, .TRUE., 1 )
333 
334 cym ne sert a rien
335 cym call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
336 
337  2 CONTINUE ! Matsuno backward or leapfrog step begins here
338 
339 c$OMP MASTER
340  itcount=itcount+1
341  if (mod(itcount,1)==1) then
342  debug=.true.
343  else
344  debug=.false.
345  endif
346 c$OMP END MASTER
347 c-----------------------------------------------------------------------
348 
349 c date:
350 c -----
351 
352 
353 c gestion des appels de la physique et des dissipations:
354 c ------------------------------------------------------
355 c
356 c ... P.Le Van ( 6/02/95 ) ....
357 
358  apphys = .false.
359  statcl = .false.
360  conser = .false.
361  apdiss = .false.
362 
363  IF( purmats ) THEN
364  ! Purely Matsuno time stepping
365  IF( mod(itau,iconser) .EQ.0.AND. forward ) conser = .true.
366  IF( mod(itau,dissip_period ).EQ.0.AND..NOT.forward )
367  s apdiss = .true.
368  IF( mod(itau,iphysiq ).EQ.0.AND..NOT.forward
369  s .and. physic ) apphys = .true.
370  ELSE
371  ! Leapfrog/Matsuno time stepping
372  IF( mod(itau ,iconser) .EQ. 0 ) conser = .true.
373  IF( mod(itau+1,dissip_period).EQ.0 .AND. .NOT. forward )
374  s apdiss = .true.
375  IF( mod(itau+1,iphysiq).EQ.0.AND.physic) apphys=.true.
376  END IF
377 
378 ! Ehouarn: for Shallow Water case (ie: 1 vertical layer),
379 ! supress dissipation step
380  if (llm.eq.1) then
381  apdiss=.false.
382  endif
383 
384 cym ---> Pour le moment
385 cym apphys = .FALSE.
386  statcl = .false.
387  conser = .false. ! ie: no output of control variables to stdout in //
388 
389  if (firstcaldyn) then
390 c$OMP MASTER
391  call setdistrib(jj_nb_caldyn)
392 c$OMP END MASTER
393 c$OMP BARRIER
394  firstcaldyn=.false.
395 cym call InitTime
396 c$OMP MASTER
397  call init_timer
398 c$OMP END MASTER
399  endif
400 
401 c$OMP MASTER
402  IF (ok_start_timer) THEN
403  CALL inittime
404  ok_start_timer=.false.
405  ENDIF
406 c$OMP END MASTER
407 
408  if (adjust) then
409 c$OMP MASTER
410  adjustcount=adjustcount+1
411  if (iapptrac==iapp_tracvl .and. (forward. or . leapf)
412  & .and. itau/iphysiq>2 .and. adjustcount>30) then
413  adjustcount=0
414  call allgather_timer_average
415 
416  if (prt_level > 9) then
417 
418  print *,'*********************************'
419  print *,'****** TIMER CALDYN ******'
420  do i=0,mpi_size-1
421  print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
422  & ' : temps moyen :',
423  & timer_average(jj_nb_caldyn(i),timer_caldyn,i),
424  & '+-',timer_delta(jj_nb_caldyn(i),timer_caldyn,i)
425  enddo
426 
427  print *,'*********************************'
428  print *,'****** TIMER VANLEER ******'
429  do i=0,mpi_size-1
430  print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
431  & ' : temps moyen :',
432  & timer_average(jj_nb_vanleer(i),timer_vanleer,i),
433  & '+-',timer_delta(jj_nb_vanleer(i),timer_vanleer,i)
434  enddo
435 
436  print *,'*********************************'
437  print *,'****** TIMER DISSIP ******'
438  do i=0,mpi_size-1
439  print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
440  & ' : temps moyen :',
441  & timer_average(jj_nb_dissip(i),timer_dissip,i),
442  & '+-',timer_delta(jj_nb_dissip(i),timer_dissip,i)
443  enddo
444 
445  if (mpi_rank==0) call writebands
446 
447  endif
448 
449  call adjustbands_caldyn
450  if (mpi_rank==0) call writebands
451 
452  call register_swapfieldhallo(ucov,ucov,ip1jmp1,llm,
453  & jj_nb_caldyn,0,0,testrequest)
454  call register_swapfieldhallo(ucovm1,ucovm1,ip1jmp1,llm,
455  & jj_nb_caldyn,0,0,testrequest)
456  call register_swapfieldhallo(vcov,vcov,ip1jm,llm,
457  & jj_nb_caldyn,0,0,testrequest)
458  call register_swapfieldhallo(vcovm1,vcovm1,ip1jm,llm,
459  & jj_nb_caldyn,0,0,testrequest)
460  call register_swapfieldhallo(teta,teta,ip1jmp1,llm,
461  & jj_nb_caldyn,0,0,testrequest)
462  call register_swapfieldhallo(tetam1,tetam1,ip1jmp1,llm,
463  & jj_nb_caldyn,0,0,testrequest)
464  call register_swapfieldhallo(masse,masse,ip1jmp1,llm,
465  & jj_nb_caldyn,0,0,testrequest)
466  call register_swapfieldhallo(massem1,massem1,ip1jmp1,llm,
467  & jj_nb_caldyn,0,0,testrequest)
468  call register_swapfieldhallo(ps,ps,ip1jmp1,1,
469  & jj_nb_caldyn,0,0,testrequest)
470  call register_swapfieldhallo(psm1,psm1,ip1jmp1,1,
471  & jj_nb_caldyn,0,0,testrequest)
472  call register_swapfieldhallo(pkf,pkf,ip1jmp1,llm,
473  & jj_nb_caldyn,0,0,testrequest)
474  call register_swapfieldhallo(pk,pk,ip1jmp1,llm,
475  & jj_nb_caldyn,0,0,testrequest)
476  call register_swapfieldhallo(pks,pks,ip1jmp1,1,
477  & jj_nb_caldyn,0,0,testrequest)
478  call register_swapfieldhallo(phis,phis,ip1jmp1,1,
479  & jj_nb_caldyn,0,0,testrequest)
480  call register_swapfieldhallo(phi,phi,ip1jmp1,llm,
481  & jj_nb_caldyn,0,0,testrequest)
482 ! call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
483 ! & jj_Nb_caldyn,0,0,TestRequest)
484 
485  do j=1,nqtot
486  call register_swapfieldhallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
487  & jj_nb_caldyn,0,0,testrequest)
488  enddo
489 
490  call setdistrib(jj_nb_caldyn)
491  call sendrequest(testrequest)
492  call waitrequest(testrequest)
493 
494  call adjustbands_dissip
495  call adjustbands_physic
496 
497  endif
498 c$OMP END MASTER
499  endif
500 
501 
502 
503 c-----------------------------------------------------------------------
504 c calcul des tendances dynamiques:
505 c --------------------------------
506 c$OMP BARRIER
507 c$OMP MASTER
508  call vtb(vthallo)
509 c$OMP END MASTER
510 
511  call register_hallo(ucov,ip1jmp1,llm,1,1,1,1,testrequest)
512  call register_hallo(vcov,ip1jm,llm,1,1,1,1,testrequest)
513  call register_hallo(teta,ip1jmp1,llm,1,1,1,1,testrequest)
514  call register_hallo(ps,ip1jmp1,1,1,2,2,1,testrequest)
515  call register_hallo(pkf,ip1jmp1,llm,1,1,1,1,testrequest)
516  call register_hallo(pk,ip1jmp1,llm,1,1,1,1,testrequest)
517  call register_hallo(pks,ip1jmp1,1,1,1,1,1,testrequest)
518  call register_hallo(p,ip1jmp1,llmp1,1,1,1,1,testrequest)
519 
520 c do j=1,nqtot
521 c call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
522 c * TestRequest)
523 c enddo
524 
525  call sendrequest(testrequest)
526 c$OMP BARRIER
527  call waitrequest(testrequest)
528 
529 c$OMP MASTER
530  call vte(vthallo)
531 c$OMP END MASTER
532 c$OMP BARRIER
533 
534  if (debug) then
535 !$OMP BARRIER
536 !$OMP MASTER
537  call writefield_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
538  call writefield_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
539  call writefield_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
540  call writefield_p('ps',reshape(ps,(/iip1,jmp1/)))
541  call writefield_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
542  call writefield_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
543  call writefield_p('pks',reshape(pks,(/iip1,jmp1/)))
544  call writefield_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
545  call writefield_p('phis',reshape(phis,(/iip1,jmp1/)))
546  do j=1,nqtot
547  call writefield_p('q'//trim(int2str(j)),
548  . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
549  enddo
550 !$OMP END MASTER
551 c$OMP BARRIER
552  endif
553 
554 
555  true_itau=true_itau+1
556 
557 c$OMP MASTER
558  IF (prt_level>9) THEN
559  WRITE(lunout,*)"leapfrog_p: Iteration No",true_itau
560  ENDIF
561 
562 
563  call start_timer(timer_caldyn)
564 
565  ! compute geopotential phi()
566  CALL geopot_p( ip1jmp1, teta , pk , pks, phis , phi )
567 
568 
569  call vtb(vtcaldyn)
570 c$OMP END MASTER
571 ! var_time=time+iday-day_ini
572 
573 c$OMP BARRIER
574 ! CALL FTRACE_REGION_BEGIN("caldyn")
575  time = jd_cur + jh_cur
576  CALL caldyn_p
577  $ ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
578  $ phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
579 
580 ! CALL FTRACE_REGION_END("caldyn")
581 
582 c$OMP MASTER
583  call vte(vtcaldyn)
584 c$OMP END MASTER
585 
586 cc$OMP BARRIER
587 cc$OMP MASTER
588 ! call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
589 ! call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
590 ! call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
591 ! call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
592 ! call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
593 ! call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
594 ! call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
595 ! call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
596 ! call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
597 ! call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
598 cc$OMP END MASTER
599 
600 c-----------------------------------------------------------------------
601 c calcul des tendances advection des traceurs (dont l'humidite)
602 c -------------------------------------------------------------
603 
604  IF( forward. or . leapf ) THEN
605 cc$OMP PARALLEL DEFAULT(SHARED)
606 c
607  CALL caladvtrac_p(q,pbaru,pbarv,
608  * p, masse, dq, teta,
609  . flxw,pk, iapptrac)
610 
611 C Stokage du flux de masse pour traceurs OFF-LINE
612  IF (offline .AND. .NOT. adjust) THEN
613  CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
614  . dtvr, itau)
615  ENDIF
616 
617  ENDIF ! of IF( forward. OR . leapf )
618 cc$OMP END PARALLEL
619 
620 c-----------------------------------------------------------------------
621 c integrations dynamique et traceurs:
622 c ----------------------------------
623 
624 c$OMP MASTER
625  call vtb(vtintegre)
626 c$OMP END MASTER
627 c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
628 c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
629 c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
630 c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
631 c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
632 c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
633 c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
634 c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
635 cc$OMP PARALLEL DEFAULT(SHARED)
636 c$OMP BARRIER
637 ! CALL FTRACE_REGION_BEGIN("integrd")
638 
639  CALL integrd_p( 2,vcovm1,ucovm1,tetam1,psm1,massem1 ,
640  $ dv,du,dteta,dq,dp,vcov,ucov,teta,q,ps,masse,phis )
641 ! $ finvmaold )
642 
643 ! CALL FTRACE_REGION_END("integrd")
644 c$OMP BARRIER
645 cc$OMP MASTER
646 c call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
647 c call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
648 c call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
649 c call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
650 c call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
651 c call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
652 c call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
653 c call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
654 c
655 c call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
656 c do j=1,nqtot
657 c call WriteField_p('q'//trim(int2str(j)),
658 c . reshape(q(:,:,j),(/iip1,jmp1,llm/)))
659 c call WriteField_p('dq'//trim(int2str(j)),
660 c . reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
661 c enddo
662 cc$OMP END MASTER
663 
664 
665 c$OMP MASTER
666  call vte(vtintegre)
667 c$OMP END MASTER
668 c .P.Le Van (26/04/94 ajout de finvpold dans l'appel d'integrd)
669 c
670 c-----------------------------------------------------------------------
671 c calcul des tendances physiques:
672 c -------------------------------
673 c ######## P.Le Van ( Modif le 6/02/95 ) ###########
674 c
675  IF( purmats ) THEN
676  IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .true.
677  ELSE
678  IF( itau+1. eq. itaufin ) lafin = .true.
679  ENDIF
680 
681 cc$OMP END PARALLEL
682 
683 c
684 c
685  IF( apphys ) THEN
686 c
687 c ....... Ajout P.Le Van ( 17/04/96 ) ...........
688 c
689 cc$OMP PARALLEL DEFAULT(SHARED)
690 cc$OMP+ PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
691 
692 c$OMP MASTER
693  call suspend_timer(timer_caldyn)
694 
695  if (prt_level >= 10) then
696  write(lunout,*)
697  & 'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
698  endif
699 c$OMP END MASTER
700 
701  CALL pression_p( ip1jmp1, ap, bp, ps, p )
702 
703 c$OMP BARRIER
704  if (pressure_exner) then
705  CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
706  else
707  CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
708  endif
709 c$OMP BARRIER
710  jd_cur = jd_ref + day_ini - day_ref
711  $ + itau/day_step
712 
713  IF (planet_type .eq."generic") THEN
714  ! AS: we make jD_cur to be pday
715  jd_cur = int(day_ini + itau/day_step)
716  ENDIF
717 
718  jh_cur = jh_ref + start_time + &
719  & mod(itau,day_step)/float(day_step)
720 ! call ju2ymds(jD_cur+jH_cur, an, mois, jour, secondes)
721  if (jh_cur > 1.0 ) then
722  jd_cur = jd_cur +1.
723  jh_cur = jh_cur -1.
724  endif
725 
726 c rajout debug
727 c lafin = .true.
728 
729 
730 c Inbterface avec les routines de phylmd (phymars ... )
731 c -----------------------------------------------------
732 
733 c+jld
734 
735 c Diagnostique de conservation de l'energie : initialisation
736  IF (ip_ebil_dyn.ge.1 ) THEN
737  ztit='bil dyn'
738 ! Ehouarn: be careful, diagedyn is Earth-specific (includes ../phylmd/..)!
739  IF (planet_type.eq."earth") THEN
740 #ifdef CPP_EARTH
741  CALL diagedyn(ztit,2,1,1,dtphys
742  & , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
743 #endif
744  ENDIF
745  ENDIF
746 c-jld
747 c$OMP BARRIER
748 c$OMP MASTER
749  call vtb(vthallo)
750 c$OMP END MASTER
751 
752  call settag(request_physic,800)
753 
754  call register_swapfieldhallo(ucov,ucov,ip1jmp1,llm,
755  * jj_nb_physic,2,2,request_physic)
756 
757  call register_swapfieldhallo(vcov,vcov,ip1jm,llm,
758  * jj_nb_physic,2,2,request_physic)
759 
760  call register_swapfieldhallo(teta,teta,ip1jmp1,llm,
761  * jj_nb_physic,2,2,request_physic)
762 
763  call register_swapfieldhallo(masse,masse,ip1jmp1,llm,
764  * jj_nb_physic,1,2,request_physic)
765 
766  call register_swapfieldhallo(p,p,ip1jmp1,llmp1,
767  * jj_nb_physic,2,2,request_physic)
768 
769  call register_swapfieldhallo(pk,pk,ip1jmp1,llm,
770  * jj_nb_physic,2,2,request_physic)
771 
772  call register_swapfieldhallo(phis,phis,ip1jmp1,1,
773  * jj_nb_physic,2,2,request_physic)
774 
775  call register_swapfieldhallo(phi,phi,ip1jmp1,llm,
776  * jj_nb_physic,2,2,request_physic)
777 
778  call register_swapfieldhallo(w,w,ip1jmp1,llm,
779  * jj_nb_physic,2,2,request_physic)
780 
781 c call SetDistrib(jj_nb_vanleer)
782  do j=1,nqtot
783 
784  call register_swapfieldhallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
785  * jj_nb_physic,2,2,request_physic)
786  enddo
787 
788  call register_swapfieldhallo(flxw,flxw,ip1jmp1,llm,
789  * jj_nb_physic,2,2,request_physic)
790 
791  call sendrequest(request_physic)
792 c$OMP BARRIER
793  call waitrequest(request_physic)
794 
795 c$OMP BARRIER
796 c$OMP MASTER
797  call setdistrib(jj_nb_physic)
798  call vte(vthallo)
799 
800  call vtb(vtphysiq)
801 c$OMP END MASTER
802 c$OMP BARRIER
803 
804 cc$OMP MASTER
805 c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
806 c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
807 c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
808 c call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
809 c call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
810 cc$OMP END MASTER
811 cc$OMP BARRIER
812 ! CALL FTRACE_REGION_BEGIN("calfis")
813  CALL calfis_p(lafin ,jd_cur, jh_cur,
814  $ ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
815  $ du,dv,dteta,dq,
816  $ flxw,
817  $ clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi )
818 ! CALL FTRACE_REGION_END("calfis")
819  ijb=ij_begin
820  ije=ij_end
821  if ( .not. pole_nord) then
822 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
823  DO l=1,llm
824  dufi_tmp(1:iip1,l) = dufi(ijb:ijb+iim,l)
825  dvfi_tmp(1:iip1,l) = dvfi(ijb:ijb+iim,l)
826  dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l)
827  dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:)
828  ENDDO
829 c$OMP END DO NOWAIT
830 
831 c$OMP MASTER
832  dpfi_tmp(1:iip1) = dpfi(ijb:ijb+iim)
833 c$OMP END MASTER
834  endif ! of if ( .not. pole_nord)
835 
836 c$OMP BARRIER
837 c$OMP MASTER
838  call setdistrib(jj_nb_physic_bis)
839 
840  call vtb(vthallo)
841 c$OMP END MASTER
842 c$OMP BARRIER
843 
844  call register_hallo(dufi,ip1jmp1,llm,
845  * 1,0,0,1,request_physic)
846 
847  call register_hallo(dvfi,ip1jm,llm,
848  * 1,0,0,1,request_physic)
849 
850  call register_hallo(dtetafi,ip1jmp1,llm,
851  * 1,0,0,1,request_physic)
852 
853  call register_hallo(dpfi,ip1jmp1,1,
854  * 1,0,0,1,request_physic)
855 
856  do j=1,nqtot
857  call register_hallo(dqfi(1,1,j),ip1jmp1,llm,
858  * 1,0,0,1,request_physic)
859  enddo
860 
861  call sendrequest(request_physic)
862 c$OMP BARRIER
863  call waitrequest(request_physic)
864 
865 c$OMP BARRIER
866 c$OMP MASTER
867  call vte(vthallo)
868 
869  call setdistrib(jj_nb_physic)
870 c$OMP END MASTER
871 c$OMP BARRIER
872  ijb=ij_begin
873  if (.not. pole_nord) then
874 
875 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
876  DO l=1,llm
877  dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
878  dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
879  dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
880  & +dtetafi_tmp(1:iip1,l)
881  dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
882  & + dqfi_tmp(1:iip1,l,:)
883  ENDDO
884 c$OMP END DO NOWAIT
885 
886 c$OMP MASTER
887  dpfi(ijb:ijb+iim) = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
888 c$OMP END MASTER
889 
890  endif ! of if (.not. pole_nord)
891 c$OMP BARRIER
892 cc$OMP MASTER
893 c call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
894 c call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
895 c call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
896 c call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
897 cc$OMP END MASTER
898 c
899 c do j=1,nqtot
900 c call WriteField_p('dqfi'//trim(int2str(j)),
901 c . reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
902 c enddo
903 
904 c ajout des tendances physiques:
905 c ------------------------------
906  IF (ok_strato) THEN
907  CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
908  ENDIF
909 
910  CALL addfi_p( dtphys, leapf, forward ,
911  $ ucov, vcov, teta , q ,ps ,
912  $ dufi, dvfi, dtetafi , dqfi ,dpfi )
913 
914 c$OMP BARRIER
915 c$OMP MASTER
916  call vte(vtphysiq)
917 
918  call vtb(vthallo)
919 c$OMP END MASTER
920 
921  call settag(request_physic,800)
922  call register_swapfield(ucov,ucov,ip1jmp1,llm,
923  * jj_nb_caldyn,request_physic)
924 
925  call register_swapfield(vcov,vcov,ip1jm,llm,
926  * jj_nb_caldyn,request_physic)
927 
928  call register_swapfield(teta,teta,ip1jmp1,llm,
929  * jj_nb_caldyn,request_physic)
930 
931  call register_swapfield(masse,masse,ip1jmp1,llm,
932  * jj_nb_caldyn,request_physic)
933 
934  call register_swapfield(p,p,ip1jmp1,llmp1,
935  * jj_nb_caldyn,request_physic)
936 
937  call register_swapfield(pk,pk,ip1jmp1,llm,
938  * jj_nb_caldyn,request_physic)
939 
940  call register_swapfield(phis,phis,ip1jmp1,1,
941  * jj_nb_caldyn,request_physic)
942 
943  call register_swapfield(phi,phi,ip1jmp1,llm,
944  * jj_nb_caldyn,request_physic)
945 
946  call register_swapfield(w,w,ip1jmp1,llm,
947  * jj_nb_caldyn,request_physic)
948 
949  do j=1,nqtot
950 
951  call register_swapfield(q(1,1,j),q(1,1,j),ip1jmp1,llm,
952  * jj_nb_caldyn,request_physic)
953 
954  enddo
955 
956  call sendrequest(request_physic)
957 c$OMP BARRIER
958  call waitrequest(request_physic)
959 
960 c$OMP BARRIER
961 c$OMP MASTER
962  call vte(vthallo)
963  call setdistrib(jj_nb_caldyn)
964 c$OMP END MASTER
965 c$OMP BARRIER
966 c
967 c Diagnostique de conservation de l'energie : difference
968  IF (ip_ebil_dyn.ge.1 ) THEN
969  ztit='bil phys'
970  CALL diagedyn(ztit,2,1,1,dtphys
971  e , ucov , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
972  ENDIF
973 
974 cc$OMP MASTER
975 c if (debug) then
976 c call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
977 c call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
978 c call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
979 c endif
980 cc$OMP END MASTER
981 
982 
983 c-jld
984 c$OMP MASTER
985  call resume_timer(timer_caldyn)
986  if (firstphysic) then
987  ok_start_timer=.true.
988  firstphysic=.false.
989  endif
990 c$OMP END MASTER
991  ENDIF ! of IF( apphys )
992 
993  IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
994 ! Academic case : Simple friction and Newtonan relaxation
995 ! -------------------------------------------------------
996  ijb=ij_begin
997  ije=ij_end
998 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
999  do l=1,llm
1000  teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
1001  & (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
1002  & (knewt_g+knewt_t(l)*clat4(ijb:ije))
1003  enddo ! of do l=1,llm
1004 !$OMP END DO
1005 
1006 !$OMP MASTER
1007  if (planet_type.eq."giant") then
1008  ! add an intrinsic heat flux at the base of the atmosphere
1009  teta(ijb:ije,1) = teta(ijb:ije,1)
1010  & + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
1011  endif
1012 !$OMP END MASTER
1013 !$OMP BARRIER
1014 
1015  call register_hallo(ucov,ip1jmp1,llm,0,1,1,0,request_physic)
1016  call register_hallo(vcov,ip1jm,llm,1,1,1,1,request_physic)
1017  call sendrequest(request_physic)
1018 c$OMP BARRIER
1019  call waitrequest(request_physic)
1020 c$OMP BARRIER
1021  call friction_p(ucov,vcov,dtvr)
1022 !$OMP BARRIER
1023 
1024  ! Sponge layer (if any)
1025  IF (ok_strato) THEN
1026  ! set dufi,dvfi,... to zero
1027  ijb=ij_begin
1028  ije=ij_end
1029 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1030  do l=1,llm
1031  dufi(ijb:ije,l)=0
1032  dtetafi(ijb:ije,l)=0
1033  dqfi(ijb:ije,l,1:nqtot)=0
1034  enddo
1035 !$OMP END DO
1036 !$OMP MASTER
1037  dpfi(ijb:ije)=0
1038 !$OMP END MASTER
1039  ijb=ij_begin
1040  ije=ij_end
1041  if (pole_sud) ije=ije-iip1
1042 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1043  do l=1,llm
1044  dvfi(ijb:ije,l)=0
1045  enddo
1046 !$OMP END DO
1047 
1048  CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
1049  CALL addfi_p( dtvr, leapf, forward ,
1050  $ ucov, vcov, teta , q ,ps ,
1051  $ dufi, dvfi, dtetafi , dqfi ,dpfi )
1052 !$OMP BARRIER
1053  ENDIF ! of IF (ok_strato)
1054  ENDIF ! of IF(iflag_phys.EQ.2)
1055 
1056 
1057  CALL pression_p( ip1jmp1, ap, bp, ps, p )
1058 c$OMP BARRIER
1059  if (pressure_exner) then
1060  CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
1061  else
1062  CALL exner_milieu_p( ip1jmp1, ps, p, beta, pks, pk, pkf )
1063  endif
1064 c$OMP BARRIER
1065 
1066 cc$OMP END PARALLEL
1067 
1068 c-----------------------------------------------------------------------
1069 c dissipation horizontale et verticale des petites echelles:
1070 c ----------------------------------------------------------
1071 
1072  IF(apdiss) THEN
1073 cc$OMP PARALLEL DEFAULT(SHARED)
1074 cc$OMP+ PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
1075 c$OMP MASTER
1076  call suspend_timer(timer_caldyn)
1077 
1078 c print*,'Entree dans la dissipation : Iteration No ',true_itau
1079 c calcul de l'energie cinetique avant dissipation
1080 c print *,'Passage dans la dissipation'
1081 
1082  call vtb(vthallo)
1083 c$OMP END MASTER
1084 
1085 c$OMP BARRIER
1086 
1087  call register_swapfieldhallo(ucov,ucov,ip1jmp1,llm,
1088  * jj_nb_dissip,1,1,request_dissip)
1089 
1090  call register_swapfieldhallo(vcov,vcov,ip1jm,llm,
1091  * jj_nb_dissip,1,1,request_dissip)
1092 
1093  call register_swapfield(teta,teta,ip1jmp1,llm,
1094  * jj_nb_dissip,request_dissip)
1095 
1096  call register_swapfield(p,p,ip1jmp1,llmp1,
1097  * jj_nb_dissip,request_dissip)
1098 
1099  call register_swapfield(pk,pk,ip1jmp1,llm,
1100  * jj_nb_dissip,request_dissip)
1101 
1102  call sendrequest(request_dissip)
1103 c$OMP BARRIER
1104  call waitrequest(request_dissip)
1105 
1106 c$OMP BARRIER
1107 c$OMP MASTER
1108  call setdistrib(jj_nb_dissip)
1109  call vte(vthallo)
1110  call vtb(vtdissipation)
1111  call start_timer(timer_dissip)
1112 c$OMP END MASTER
1113 c$OMP BARRIER
1114 
1115  call covcont_p(llm,ucov,vcov,ucont,vcont)
1116  call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1117 
1118 c dissipation
1119 
1120 ! CALL FTRACE_REGION_BEGIN("dissip")
1121  CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1122 ! CALL FTRACE_REGION_END("dissip")
1123 
1124  ijb=ij_begin
1125  ije=ij_end
1126 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1127  DO l=1,llm
1128  ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1129  ENDDO
1130 c$OMP END DO NOWAIT
1131  if (pole_sud) ije=ije-iip1
1132 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1133  DO l=1,llm
1134  vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1135  ENDDO
1136 c$OMP END DO NOWAIT
1137 
1138 c teta=teta+dtetadis
1139 
1140 
1141 c------------------------------------------------------------------------
1142  if (dissip_conservative) then
1143 C On rajoute la tendance due a la transform. Ec -> E therm. cree
1144 C lors de la dissipation
1145 c$OMP BARRIER
1146 c$OMP MASTER
1147  call suspend_timer(timer_dissip)
1148  call vtb(vthallo)
1149 c$OMP END MASTER
1150  call register_hallo(ucov,ip1jmp1,llm,1,1,1,1,request_dissip)
1151  call register_hallo(vcov,ip1jm,llm,1,1,1,1,request_dissip)
1152  call sendrequest(request_dissip)
1153 c$OMP BARRIER
1154  call waitrequest(request_dissip)
1155 c$OMP MASTER
1156  call vte(vthallo)
1157  call resume_timer(timer_dissip)
1158 c$OMP END MASTER
1159 c$OMP BARRIER
1160  call covcont_p(llm,ucov,vcov,ucont,vcont)
1161  call enercin_p(vcov,ucov,vcont,ucont,ecin)
1162 
1163  ijb=ij_begin
1164  ije=ij_end
1165 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1166  do l=1,llm
1167  do ij=ijb,ije
1168  dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1169  dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1170  enddo
1171  enddo
1172 c$OMP END DO NOWAIT
1173  endif ! of if (dissip_conservative)
1174 
1175  ijb=ij_begin
1176  ije=ij_end
1177 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1178  do l=1,llm
1179  do ij=ijb,ije
1180  teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1181  enddo
1182  enddo
1183 c$OMP END DO NOWAIT
1184 c------------------------------------------------------------------------
1185 
1186 
1187 c ....... P. Le Van ( ajout le 17/04/96 ) ...........
1188 c ... Calcul de la valeur moyenne, unique de h aux poles .....
1189 c
1190 
1191  ijb=ij_begin
1192  ije=ij_end
1193 
1194  if (pole_nord) then
1195 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1196  DO l = 1, llm
1197  DO ij = 1,iim
1198  tppn(ij) = aire( ij ) * teta( ij ,l)
1199  ENDDO
1200  tpn = ssum(iim,tppn,1)/apoln
1201 
1202  DO ij = 1, iip1
1203  teta( ij ,l) = tpn
1204  ENDDO
1205  ENDDO
1206 c$OMP END DO NOWAIT
1207 
1208  if (1 == 0) then
1209 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1210 !!! 2) should probably not be here anyway
1211 !!! but are kept for those who would want to revert to previous behaviour
1212 c$OMP MASTER
1213  DO ij = 1,iim
1214  tppn(ij) = aire( ij ) * ps( ij )
1215  ENDDO
1216  tpn = ssum(iim,tppn,1)/apoln
1217 
1218  DO ij = 1, iip1
1219  ps( ij ) = tpn
1220  ENDDO
1221 c$OMP END MASTER
1222  endif ! of if (1 == 0)
1223  endif ! of of (pole_nord)
1224 
1225  if (pole_sud) then
1226 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1227  DO l = 1, llm
1228  DO ij = 1,iim
1229  tpps(ij) = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1230  ENDDO
1231  tps = ssum(iim,tpps,1)/apols
1232 
1233  DO ij = 1, iip1
1234  teta(ij+ip1jm,l) = tps
1235  ENDDO
1236  ENDDO
1237 c$OMP END DO NOWAIT
1238 
1239  if (1 == 0) then
1240 !!! Ehouarn: lines here 1) kill 1+1=2 in the dynamics
1241 !!! 2) should probably not be here anyway
1242 !!! but are kept for those who would want to revert to previous behaviour
1243 c$OMP MASTER
1244  DO ij = 1,iim
1245  tpps(ij) = aire(ij+ip1jm) * ps(ij+ip1jm)
1246  ENDDO
1247  tps = ssum(iim,tpps,1)/apols
1248 
1249  DO ij = 1, iip1
1250  ps(ij+ip1jm) = tps
1251  ENDDO
1252 c$OMP END MASTER
1253  endif ! of if (1 == 0)
1254  endif ! of if (pole_sud)
1255 
1256 
1257 c$OMP BARRIER
1258 c$OMP MASTER
1259  call vte(vtdissipation)
1260 
1261  call stop_timer(timer_dissip)
1262 
1263  call vtb(vthallo)
1264 c$OMP END MASTER
1265  call register_swapfield(ucov,ucov,ip1jmp1,llm,
1266  * jj_nb_caldyn,request_dissip)
1267 
1268  call register_swapfield(vcov,vcov,ip1jm,llm,
1269  * jj_nb_caldyn,request_dissip)
1270 
1271  call register_swapfield(teta,teta,ip1jmp1,llm,
1272  * jj_nb_caldyn,request_dissip)
1273 
1274  call register_swapfield(p,p,ip1jmp1,llmp1,
1275  * jj_nb_caldyn,request_dissip)
1276 
1277  call register_swapfield(pk,pk,ip1jmp1,llm,
1278  * jj_nb_caldyn,request_dissip)
1279 
1280  call sendrequest(request_dissip)
1281 c$OMP BARRIER
1282  call waitrequest(request_dissip)
1283 
1284 c$OMP BARRIER
1285 c$OMP MASTER
1286  call setdistrib(jj_nb_caldyn)
1287  call vte(vthallo)
1288  call resume_timer(timer_caldyn)
1289 c print *,'fin dissipation'
1290 c$OMP END MASTER
1291 c$OMP BARRIER
1292  END IF ! of IF(apdiss)
1293 
1294 cc$OMP END PARALLEL
1295 
1296 c ajout debug
1297 c IF( lafin ) then
1298 c abort_message = 'Simulation finished'
1299 c call abort_gcm(modname,abort_message,0)
1300 c ENDIF
1301 
1302 c ********************************************************************
1303 c ********************************************************************
1304 c .... fin de l'integration dynamique et physique pour le pas itau ..
1305 c ********************************************************************
1306 c ********************************************************************
1307 
1308 c preparation du pas d'integration suivant ......
1309 cym call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1310 cym call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1311 c$OMP MASTER
1312  call stop_timer(timer_caldyn)
1313 c$OMP END MASTER
1314  IF (itau==itaumax) then
1315 c$OMP MASTER
1316  call allgather_timer_average
1317 
1318  if (mpi_rank==0) then
1319 
1320  print *,'*********************************'
1321  print *,'****** TIMER CALDYN ******'
1322  do i=0,mpi_size-1
1323  print *,'proc',i,' : Nb Bandes :',jj_nb_caldyn(i),
1324  & ' : temps moyen :',
1325  & timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1326  enddo
1327 
1328  print *,'*********************************'
1329  print *,'****** TIMER VANLEER ******'
1330  do i=0,mpi_size-1
1331  print *,'proc',i,' : Nb Bandes :',jj_nb_vanleer(i),
1332  & ' : temps moyen :',
1333  & timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1334  enddo
1335 
1336  print *,'*********************************'
1337  print *,'****** TIMER DISSIP ******'
1338  do i=0,mpi_size-1
1339  print *,'proc',i,' : Nb Bandes :',jj_nb_dissip(i),
1340  & ' : temps moyen :',
1341  & timer_average(jj_nb_dissip(i),timer_dissip,i)
1342  enddo
1343 
1344  print *,'*********************************'
1345  print *,'****** TIMER PHYSIC ******'
1346  do i=0,mpi_size-1
1347  print *,'proc',i,' : Nb Bandes :',jj_nb_physic(i),
1348  & ' : temps moyen :',
1349  & timer_average(jj_nb_physic(i),timer_physic,i)
1350  enddo
1351 
1352  endif
1353 
1354  print *,'Taille du Buffer MPI (REAL*8)',maxbuffersize
1355  print *,'Taille du Buffer MPI utilise (REAL*8)',maxbuffersize_used
1356  print *, 'Temps total ecoule sur la parallelisation :',difftime()
1357  print *, 'Temps CPU ecoule sur la parallelisation :',diffcputime()
1358  CALL print_filtre_timer
1359  call fin_getparam
1360  call finalize_parallel
1361 c$OMP END MASTER
1362 c$OMP BARRIER
1363  RETURN
1364  ENDIF
1365 
1366  IF ( .NOT.purmats ) THEN
1367 c ........................................................
1368 c .............. schema matsuno + leapfrog ..............
1369 c ........................................................
1370 
1371  IF(forward. or. leapf) THEN
1372  itau= itau + 1
1373 ! iday= day_ini+itau/day_step
1374 ! time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1375 ! IF(time.GT.1.) THEN
1376 ! time = time-1.
1377 ! iday = iday+1
1378 ! ENDIF
1379  ENDIF
1380 
1381 
1382  IF( itau. eq. itaufinp1 ) then
1383 
1384  if (flag_verif) then
1385  write(79,*) 'ucov',ucov
1386  write(80,*) 'vcov',vcov
1387  write(81,*) 'teta',teta
1388  write(82,*) 'ps',ps
1389  write(83,*) 'q',q
1390  WRITE(85,*) 'q1 = ',q(:,:,1)
1391  WRITE(86,*) 'q3 = ',q(:,:,3)
1392  endif
1393 
1394 
1395 c$OMP MASTER
1396  call fin_getparam
1397  call finalize_parallel
1398 c$OMP END MASTER
1399  abort_message = 'Simulation finished'
1400  call abort_gcm(modname,abort_message,0)
1401  RETURN
1402  ENDIF
1403 c-----------------------------------------------------------------------
1404 c ecriture du fichier histoire moyenne:
1405 c -------------------------------------
1406 
1407  IF(mod(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1408 c$OMP BARRIER
1409  IF(itau.EQ.itaufin) THEN
1410  iav=1
1411  ELSE
1412  iav=0
1413  ENDIF
1414 #ifdef CPP_IOIPSL
1415  IF (ok_dynzon) THEN
1416  call register_hallo(vcov,ip1jm,llm,1,0,0,1,testrequest)
1417  call sendrequest(testrequest)
1418 c$OMP BARRIER
1419  call waitrequest(testrequest)
1420 c$OMP BARRIER
1421 c$OMP MASTER
1422 ! CALL writedynav_p(histaveid, itau,vcov ,
1423 ! , ucov,teta,pk,phi,q,masse,ps,phis)
1424 
1425 c ATTENTION!!! bilan_dyn_p ne marche probablement pas avec OpenMP
1426  CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1427  , ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1428 c$OMP END MASTER
1429  ENDIF !ok_dynzon
1430 #endif
1431  IF (ok_dyn_ave) THEN
1432 !$OMP MASTER
1433 #ifdef CPP_IOIPSL
1434 ! Ehouarn: Gather fields and make master send to output
1435  call gather_field(vcov,ip1jm,llm,0)
1436  call gather_field(ucov,ip1jmp1,llm,0)
1437  call gather_field(teta,ip1jmp1,llm,0)
1438  call gather_field(pk,ip1jmp1,llm,0)
1439  call gather_field(phi,ip1jmp1,llm,0)
1440  do iq=1,nqtot
1441  call gather_field(q(1,1,iq),ip1jmp1,llm,0)
1442  enddo
1443  call gather_field(masse,ip1jmp1,llm,0)
1444  call gather_field(ps,ip1jmp1,1,0)
1445  call gather_field(phis,ip1jmp1,1,0)
1446  if (mpi_rank==0) then
1447  CALL writedynav(itau,vcov,
1448  & ucov,teta,pk,phi,q,masse,ps,phis)
1449  endif
1450 #endif
1451 !$OMP END MASTER
1452  ENDIF ! of IF (ok_dyn_ave)
1453  ENDIF ! of IF((MOD(itau,iperiod).EQ.0).OR.(itau.EQ.itaufin))
1454 
1455 c-----------------------------------------------------------------------
1456 c ecriture de la bande histoire:
1457 c ------------------------------
1458 
1459  IF( mod(itau,iecri).EQ.0) THEN
1460  ! Ehouarn: output only during LF or Backward Matsuno
1461  if (leapf.or.(.not.leapf.and.(.not.forward))) then
1462 c$OMP BARRIER
1463 c$OMP MASTER
1464  CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1465 
1466 cym unat=0.
1467 
1468  ijb=ij_begin
1469  ije=ij_end
1470 
1471  if (pole_nord) then
1472  ijb=ij_begin+iip1
1473  unat(1:iip1,:)=0.
1474  endif
1475 
1476  if (pole_sud) then
1477  ije=ij_end-iip1
1478  unat(ij_end-iip1+1:ij_end,:)=0.
1479  endif
1480 
1481  do l=1,llm
1482  unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1483  enddo
1484 
1485  ijb=ij_begin
1486  ije=ij_end
1487  if (pole_sud) ije=ij_end-iip1
1488 
1489  do l=1,llm
1490  vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1491  enddo
1492 
1493 #ifdef CPP_IOIPSL
1494  if (ok_dyn_ins) then
1495 ! Ehouarn: Gather fields and make master write to output
1496  call gather_field(vcov,ip1jm,llm,0)
1497  call gather_field(ucov,ip1jmp1,llm,0)
1498  call gather_field(teta,ip1jmp1,llm,0)
1499  call gather_field(phi,ip1jmp1,llm,0)
1500  do iq=1,nqtot
1501  call gather_field(q(1,1,iq),ip1jmp1,llm,0)
1502  enddo
1503  call gather_field(masse,ip1jmp1,llm,0)
1504  call gather_field(ps,ip1jmp1,1,0)
1505  call gather_field(phis,ip1jmp1,1,0)
1506  if (mpi_rank==0) then
1507  CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1508  endif
1509 ! CALL writehist_p(histid,histvid, itau,vcov,
1510 ! & ucov,teta,phi,q,masse,ps,phis)
1511 ! or use writefield_p
1512 ! call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1513 ! call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1514 ! call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
1515 ! call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
1516  endif ! of if (ok_dyn_ins)
1517 #endif
1518 ! For some Grads outputs of fields
1519  if (output_grads_dyn) then
1520 ! Ehouarn: hope this works the way I think it does:
1521  call gather_field(unat,ip1jmp1,llm,0)
1522  call gather_field(vnat,ip1jm,llm,0)
1523  call gather_field(teta,ip1jmp1,llm,0)
1524  call gather_field(ps,ip1jmp1,1,0)
1525  do iq=1,nqtot
1526  call gather_field(q(1,1,iq),ip1jmp1,llm,0)
1527  enddo
1528  if (mpi_rank==0) then
1529 #include "write_grads_dyn.h"
1530  endif
1531  endif ! of if (output_grads_dyn)
1532 c$OMP END MASTER
1533  endif ! of if (leapf.or.(.not.leapf.and.(.not.forward)))
1534  ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1535 
1536  IF(itau.EQ.itaufin) THEN
1537 
1538 c$OMP BARRIER
1539 c$OMP MASTER
1540 
1541 ! if (planet_type.eq."earth") then
1542 ! Write an Earth-format restart file
1543  CALL dynredem1_p("restart.nc",0.0,
1544  & vcov,ucov,teta,q,masse,ps)
1545 ! endif ! of if (planet_type.eq."earth")
1546 
1547 ! CLOSE(99)
1548 c$OMP END MASTER
1549  ENDIF ! of IF (itau.EQ.itaufin)
1550 
1551 c-----------------------------------------------------------------------
1552 c gestion de l'integration temporelle:
1553 c ------------------------------------
1554 
1555  IF( mod(itau,iperiod).EQ.0 ) THEN
1556  go to 1
1557  ELSE IF ( mod(itau-1,iperiod). eq. 0 ) THEN
1558 
1559  IF( forward ) THEN
1560 c fin du pas forward et debut du pas backward
1561 
1562  forward = .false.
1563  leapf = .false.
1564  go to 2
1565 
1566  ELSE
1567 c fin du pas backward et debut du premier pas leapfrog
1568 
1569  leapf = .true.
1570  dt = 2.*dtvr
1571  go to 2
1572  END IF
1573  ELSE
1574 
1575 c ...... pas leapfrog .....
1576 
1577  leapf = .true.
1578  dt = 2.*dtvr
1579  go to 2
1580  END IF ! of IF (MOD(itau,iperiod).EQ.0)
1581  ! ELSEIF (MOD(itau-1,iperiod).EQ.0)
1582 
1583 
1584  ELSE ! of IF (.not.purmats)
1585 
1586 c ........................................................
1587 c .............. schema matsuno ...............
1588 c ........................................................
1589  IF( forward ) THEN
1590 
1591  itau = itau + 1
1592 ! iday = day_ini+itau/day_step
1593 ! time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1594 !
1595 ! IF(time.GT.1.) THEN
1596 ! time = time-1.
1597 ! iday = iday+1
1598 ! ENDIF
1599 
1600  forward = .false.
1601  IF( itau. eq. itaufinp1 ) then
1602 c$OMP MASTER
1603  call fin_getparam
1604  call finalize_parallel
1605 c$OMP END MASTER
1606  abort_message = 'Simulation finished'
1607  call abort_gcm(modname,abort_message,0)
1608  RETURN
1609  ENDIF
1610  go to 2
1611 
1612  ELSE ! of IF(forward) i.e. backward step
1613 
1614  IF(mod(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1615  IF(itau.EQ.itaufin) THEN
1616  iav=1
1617  ELSE
1618  iav=0
1619  ENDIF
1620 #ifdef CPP_IOIPSL
1621  IF (ok_dynzon) THEN
1622 c$OMP BARRIER
1623  call register_hallo(vcov,ip1jm,llm,1,0,0,1,testrequest)
1624  call sendrequest(testrequest)
1625 c$OMP BARRIER
1626  call waitrequest(testrequest)
1627 c$OMP BARRIER
1628 c$OMP MASTER
1629 ! CALL writedynav_p(histaveid, itau,vcov ,
1630 ! , ucov,teta,pk,phi,q,masse,ps,phis)
1631  CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1632  , ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1633 c$OMP END MASTER
1634  END IF !ok_dynzon
1635 #endif
1636  IF (ok_dyn_ave) THEN
1637 !$OMP MASTER
1638 #ifdef CPP_IOIPSL
1639 ! Ehouarn: Gather fields and make master send to output
1640  call gather_field(vcov,ip1jm,llm,0)
1641  call gather_field(ucov,ip1jmp1,llm,0)
1642  call gather_field(teta,ip1jmp1,llm,0)
1643  call gather_field(pk,ip1jmp1,llm,0)
1644  call gather_field(phi,ip1jmp1,llm,0)
1645  do iq=1,nqtot
1646  call gather_field(q(1,1,iq),ip1jmp1,llm,0)
1647  enddo
1648  call gather_field(masse,ip1jmp1,llm,0)
1649  call gather_field(ps,ip1jmp1,1,0)
1650  call gather_field(phis,ip1jmp1,1,0)
1651  if (mpi_rank==0) then
1652  CALL writedynav(itau,vcov,
1653  & ucov,teta,pk,phi,q,masse,ps,phis)
1654  endif
1655 #endif
1656 !$OMP END MASTER
1657  ENDIF ! of IF (ok_dyn_ave)
1658 
1659  ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1660 
1661 
1662  IF(mod(itau,iecri ).EQ.0) THEN
1663 c IF(MOD(itau,iecri*day_step).EQ.0) THEN
1664 c$OMP BARRIER
1665 c$OMP MASTER
1666  CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1667 
1668 cym unat=0.
1669  ijb=ij_begin
1670  ije=ij_end
1671 
1672  if (pole_nord) then
1673  ijb=ij_begin+iip1
1674  unat(1:iip1,:)=0.
1675  endif
1676 
1677  if (pole_sud) then
1678  ije=ij_end-iip1
1679  unat(ij_end-iip1+1:ij_end,:)=0.
1680  endif
1681 
1682  do l=1,llm
1683  unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1684  enddo
1685 
1686  ijb=ij_begin
1687  ije=ij_end
1688  if (pole_sud) ije=ij_end-iip1
1689 
1690  do l=1,llm
1691  vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1692  enddo
1693 
1694 #ifdef CPP_IOIPSL
1695  if (ok_dyn_ins) then
1696 ! Ehouarn: Gather fields and make master send to output
1697  call gather_field(vcov,ip1jm,llm,0)
1698  call gather_field(ucov,ip1jmp1,llm,0)
1699  call gather_field(teta,ip1jmp1,llm,0)
1700  call gather_field(phi,ip1jmp1,llm,0)
1701  do iq=1,nqtot
1702  call gather_field(q(1,1,iq),ip1jmp1,llm,0)
1703  enddo
1704  call gather_field(masse,ip1jmp1,llm,0)
1705  call gather_field(ps,ip1jmp1,1,0)
1706  call gather_field(phis,ip1jmp1,1,0)
1707  if (mpi_rank==0) then
1708  CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
1709  endif
1710 ! CALL writehist_p(histid, histvid, itau,vcov ,
1711 ! & ucov,teta,phi,q,masse,ps,phis)
1712  endif ! of if (ok_dyn_ins)
1713 #endif
1714 ! For some Grads output (but does it work?)
1715  if (output_grads_dyn) then
1716  call gather_field(unat,ip1jmp1,llm,0)
1717  call gather_field(vnat,ip1jm,llm,0)
1718  call gather_field(teta,ip1jmp1,llm,0)
1719  call gather_field(ps,ip1jmp1,1,0)
1720  do iq=1,nqtot
1721  call gather_field(q(1,1,iq),ip1jmp1,llm,0)
1722  enddo
1723 c
1724  if (mpi_rank==0) then
1725 #include "write_grads_dyn.h"
1726  endif
1727  endif ! of if (output_grads_dyn)
1728 
1729 c$OMP END MASTER
1730  ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1731 
1732  IF(itau.EQ.itaufin) THEN
1733 ! if (planet_type.eq."earth") then
1734 c$OMP MASTER
1735  CALL dynredem1_p("restart.nc",0.0,
1736  . vcov,ucov,teta,q,masse,ps)
1737 c$OMP END MASTER
1738 ! endif ! of if (planet_type.eq."earth")
1739  ENDIF ! of IF(itau.EQ.itaufin)
1740 
1741  forward = .true.
1742  go to 1
1743 
1744  ENDIF ! of IF (forward)
1745 
1746  END IF ! of IF(.not.purmats)
1747 c$OMP MASTER
1748  call fin_getparam
1749  call finalize_parallel
1750 c$OMP END MASTER
1751  RETURN
1752  END