source: LMDZ5/branches/IPSLCM5A2.1/libf/dyn3dmem/leapfrog_loc.F

Last change on this file was 2616, checked in by fhourdin, 8 years ago

Reactivation de la clé "conser" pour pouvoir faire ds sorties dans
le listing tous les iconser.

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