source: LMDZ5/branches/IPSLCM5A2.1/libf/dyn3dpar/leapfrog_p.F @ 3128

Last change on this file since 3128 was 2603, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn logic.h into module logic_mod.F90
EM

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