source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F @ 245

Last change on this file since 245 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

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