source: LMDZ5/trunk/libf/dyn3dmem/leapfrog_loc.F @ 2180

Last change on this file since 2180 was 2180, checked in by acozic, 10 years ago

After checking I reload commit 2169 and 2170

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