source: LMDZ6/branches/LMDZ_ECRad/libf/dyn3dmem/leapfrog_loc.F @ 4727

Last change on this file since 4727 was 4727, checked in by idelkadi, 8 months ago

Merged trunk changes -r4488:4726 LMDZ_ECRad branch

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