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

Last change on this file since 1832 was 1825, checked in by Ehouarn Millour, 11 years ago

Première étape de l'implémentation de XIOS. Modifications isolées dans des flags CPP_XIOS. Sorties opérationnelles (sauf stations et régionalisation) en modes séquentiel et omp, pas mpi.
UG
...........................................
First step of the XIOS implementation. Modifications are confined into CPP_XIOS flags. Output is operationnal (except for stations and regionalization) in sequential and omp modes (not mpi).
UG

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