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

Last change on this file since 3405 was 2622, checked in by Ehouarn Millour, 8 years ago

Some code tidying: turn ener.h into ener_mod.F90
EM

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