source: LMDZ5/branches/LF-private/libf/dyn3dmem/leapfrog_loc.F @ 3856

Last change on this file since 3856 was 1848, checked in by yann meurdesoif, 11 years ago

Solve performance problem comming from declarations of a derived type.

YM

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