source: LMDZ4/trunk/libf/dyn3dpar/leapfrog_p.F @ 993

Last change on this file since 993 was 985, checked in by Laurent Fairhead, 17 years ago

Mise a jour de dyn3dpar par rapport a dyn3d, inclusion OpenMP et filtre FFT YM
LF

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