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

Last change on this file since 3567 was 2549, checked in by emillour, 3 years ago

Common dynamics:
Fix a weird issue that only arises in mixed MPI/OpenMP mode with
deallocation (via deallocate_buffer called by WaitRequest?) if
only the halo of surface pressure is registered. Adding the registration
of a companion 3D field seems to solve the problem.
EM

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