source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F @ 1329

Last change on this file since 1329 was 1322, checked in by Laurent Fairhead, 15 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file


Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

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