source: trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F @ 832

Last change on this file since 832 was 776, checked in by emillour, 12 years ago

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