source: LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F @ 1897

Last change on this file since 1897 was 1403, checked in by Laurent Fairhead, 14 years ago

Merged LMDZ4V5.0-dev branch changes r1292:r1399 to trunk.

Validation:
Validation consisted in compiling the HEAD revision of the trunk,
LMDZ4V5.0-dev branch and the merged sources and running different
configurations on local and SX8 machines comparing results.

Local machine: bench configuration, 32x24x11, gfortran

  • IPSLCM5A configuration (comparison between trunk and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent
  • MH07 configuration, new physics package (comparison between LMDZ4V5.0-dev branch and merged sources):
    • numerical convergence on dynamical fields over 3 days
    • start files are equivalent (except for RN and PB fields)
    • daily history files equivalent

SX8 machine (brodie), 96x95x39 on 4 processors:

  • IPSLCM5A configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent
  • MH07 configuration:
    • start files are equivalent (except for RN and PB fields)
    • monthly history files equivalent

Changes to the makegcm and create_make_gcm scripts to take into account
main programs in F90 files


Fusion de la branche LMDZ4V5.0-dev (r1292:r1399) au tronc principal

Validation:
La validation a consisté à compiler la HEAD de le trunk et de la banche
LMDZ4V5.0-dev et les sources fusionnées et de faire tourner le modéle selon
différentes configurations en local et sur SX8 et de comparer les résultats

En local: 32x24x11, config bench/gfortran

  • pour une config IPSLCM5A (comparaison tronc/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux (à part sur RN et Pb)
    • fichiers histoire égaux
  • pour une config nlle physique (MH07) (comparaison LMDZ4v5.0-dev/fusion):
    • convergence numérique sur les champs dynamiques après 3 jours
    • restart et restartphy égaux
    • fichiers histoire équivalents

Sur brodie, 96x95x39 sur 4 proc:

  • pour une config IPSLCM5A:
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc
  • pour une config MH07
    • restart et restartphy égaux (à part sur RN et PB)
    • pas de différence dans les fichiers histmth.nc

Changement sur makegcm et create_make-gcm pour pouvoir prendre en compte des
programmes principaux en *F90

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