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

Last change on this file since 1453 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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