source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3dpar/leapfrog_p.F @ 3400

Last change on this file since 3400 was 1717, checked in by Ehouarn Millour, 12 years ago

Added "arch" files for Ada (using dynamic libraries for NetCDF, you must have
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/smplocal/pub/NetCDF/4.1.3/lib:/smplocal/pub/HDF5/1.8.9/seq/lib
in your .bashrc or .bash_login or in your job to run).
Also updated some sources so that gcm bench runs in "debug" mode (note that all these changes are minor and have already been implemented in LMDZ5 trunk).
EM

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