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

Last change on this file since 1659 was 1659, checked in by Laurent Fairhead, 12 years ago

Modification pour etre consistant


Modification for consistency

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