source: LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F @ 2302

Last change on this file since 2302 was 2302, checked in by Ehouarn Millour, 9 years ago

Move etat0phys_netcdf.F90 to "dynlonlat_phylonlat/phylmd" as it relies on "phylmd" routines.
Some cleanup to remove obsolete and unecessary CPP_EARTH preprocessing condition.
EM

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