source: LMDZ6/branches/LMDZ-COSP/libf/dyn3dpar/leapfrog_p.F @ 5409

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

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