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

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

Some code tidying: turn ener.h into ener_mod.F90
EM

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