source: LMDZ6/trunk/libf/dyn3dmem/leapfrog_loc.F @ 4139

Last change on this file since 4139 was 4139, checked in by acozic, 2 years ago

added (and call) a routine for the calculation of mass flux used in offline mode

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