source: LMDZ6/branches/LMDZ-tracers/libf/dyn3dmem/leapfrog_loc.F @ 3852

Last change on this file since 3852 was 3852, checked in by dcugnet, 3 years ago

Extension of the tracers management.

The tracers files can be:

1) "traceur.def": old format, with:

  • the number of tracers on the first line
  • one line for each tracer: <tracer name> <hadv> <vadv> [<parent name>]

2) "tracer.def": new format with one section each model component.
3) "tracer_<name>.def": new format with a single section.

The formats 2 and 3 reading is driven by the "type_trac" key, which can be a

coma-separated list of components.

  • Format 2: read the sections from the "tracer.def" file.
  • format 3: read one section each "tracer_<section name>.def" file.
  • the first line of a section is "&<section name>
  • the other lines start with a tracer name followed by <key>=<val> pairs.
  • the "default" tracer name is reserved ; the other tracers of the section inherit its <key>=<val>, except for the keys that are redefined locally.

This format helps keeping the tracers files compact, thanks to the "default"
special tracer and the three levels of factorization:

  • on the tracers names: a tracer name can be a coma-separated list of tracers => all the tracers of the list have the same <key>=<val> properties
  • on the parents names: the value of the "parent" property can be a coma-separated list of tracers => only possible for geographic tagging tracers
  • on the phases: the property "phases" is [g](l][s] (gas/liquid/solid)

Read information is stored in the vector "tracers(:)", of derived type "tra".

"isotopes_params.def" is a similar file, with one section each isotopes family.
It contains a database of isotopes properties ; if there are second generation
tracers (isotopes), the corresponding sections are read.

Read information is stored in the vector "isotopes(:)", of derived type "iso".

The "getKey" function helps to get the values of the parameters stored in
"tracers" or "isotopes".

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