source: LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F @ 1592

Last change on this file since 1592 was 1592, checked in by Ehouarn Millour, 13 years ago

Minor cleanup: remove unused variables nbetat, nbetatmoy and nbetatdem.
EM

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