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
Line 
1!
2! $Id: leapfrog_p.F 1322 2010-03-12 10:54:11Z fairhead $
3!
4c
5c
6
7      SUBROUTINE leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
8     &                    time_0)
9
10       USE misc_mod
11       USE parallel
12       USE times
13       USE mod_hallo
14       USE Bands
15       USE Write_Field
16       USE Write_Field_p
17       USE vampir
18       USE timer_filtre, ONLY : print_filtre_timer
19       USE infotrac
20       USE guide_p_mod, ONLY : guide_main
21       USE getparam
22       USE control_mod
23
24      IMPLICIT NONE
25
26c      ......   Version  du 10/01/98    ..........
27
28c             avec  coordonnees  verticales hybrides
29c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
30
31c=======================================================================
32c
33c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
34c   -------
35c
36c   Objet:
37c   ------
38c
39c   GCM LMD nouvelle grille
40c
41c=======================================================================
42c
43c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
44c      et possibilite d'appeler une fonction f(y)  a derivee tangente
45c      hyperbolique a la  place de la fonction a derivee sinusoidale.
46
47c  ... Possibilite de choisir le shema pour l'advection de
48c        q  , en modifiant iadv dans traceur.def  (10/02) .
49c
50c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
51c      Pour Van-Leer iadv=10
52c
53c-----------------------------------------------------------------------
54c   Declarations:
55c   -------------
56
57#include "dimensions.h"
58#include "paramet.h"
59#include "comconst.h"
60#include "comdissnew.h"
61#include "comvert.h"
62#include "comgeom.h"
63#include "logic.h"
64#include "temps.h"
65#include "ener.h"
66#include "description.h"
67#include "serre.h"
68#include "com_io_dyn.h"
69#include "iniprint.h"
70#include "academic.h"
71     
72      INTEGER         longcles
73      PARAMETER     ( longcles = 20 )
74      REAL  clesphy0( longcles )
75
76      real zqmin,zqmax
77      INTEGER nbetatmoy, nbetatdem,nbetat
78
79c   variables dynamiques
80      REAL :: vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
81      REAL :: teta(ip1jmp1,llm)                 ! temperature potentielle
82      REAL :: q(ip1jmp1,llm,nqtot)              ! champs advectes
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
92
93c variables dynamiques intermediaire pour le transport
94      REAL,SAVE :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm) !flux de masse
95
96c   variables dynamiques au pas -1
97      REAL,SAVE :: vcovm1(ip1jm,llm),ucovm1(ip1jmp1,llm)
98      REAL,SAVE :: tetam1(ip1jmp1,llm),psm1(ip1jmp1)
99      REAL,SAVE :: massem1(ip1jmp1,llm)
100
101c   tendances dynamiques
102      REAL,SAVE :: dv(ip1jm,llm),du(ip1jmp1,llm)
103      REAL,SAVE :: dteta(ip1jmp1,llm),dp(ip1jmp1)
104      REAL,DIMENSION(:,:,:), ALLOCATABLE, SAVE :: dq
105
106c   tendances de la dissipation
107      REAL,SAVE :: dvdis(ip1jm,llm),dudis(ip1jmp1,llm)
108      REAL,SAVE :: dtetadis(ip1jmp1,llm)
109
110c   tendances physiques
111      REAL,SAVE :: dvfi(ip1jm,llm),dufi(ip1jmp1,llm)
112      REAL,SAVE :: dtetafi(ip1jmp1,llm)
113      REAL,SAVE :: dpfi(ip1jmp1)
114      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
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
122!      INTEGER  iday ! jour julien
123      REAL       time
124
125      REAL  SSUM
126      REAL time_0
127      REAL,SAVE :: finvmaold(ip1jmp1,llm)
128
129cym      LOGICAL  lafin
130      LOGICAL :: lafin
131      INTEGER ij,iq,l
132      INTEGER ik
133
134      real time_step, t_wrt, t_ops
135
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
142      LOGICAL first,callinigrads
143
144      data callinigrads/.true./
145      character*10 string10
146
147      REAL,SAVE :: alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
148      REAL,SAVE :: flxw(ip1jmp1,llm) ! flux de masse verticale
149
150c+jld variables test conservation energie
151      REAL,SAVE :: ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
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
155      REAL,SAVE :: dtetaecdt(ip1jmp1,llm)
156      REAL,SAVE :: vcont(ip1jm,llm),ucont(ip1jmp1,llm)
157      REAL,SAVE :: vnat(ip1jm,llm),unat(ip1jmp1,llm)
158      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
159      CHARACTER*15 ztit
160!      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
161!      SAVE      ip_ebil_dyn
162!      DATA      ip_ebil_dyn/0/
163c-jld
164
165      character*80 dynhist_file, dynhistave_file
166      character*20 modname
167      character*80 abort_message
168
169
170      logical,PARAMETER :: dissip_conservative=.TRUE.
171 
172      INTEGER testita
173      PARAMETER (testita = 9)
174
175      logical , parameter :: flag_verif = .false.
176     
177c declaration liees au parallelisme
178      INTEGER :: ierr
179      LOGICAL :: FirstCaldyn
180      LOGICAL :: FirstPhysic
181      INTEGER :: ijb,ije,j,i
182      type(Request) :: TestRequest
183      type(Request) :: Request_Dissip
184      type(Request) :: Request_physic
185      REAL,SAVE :: dvfi_tmp(iip1,llm),dufi_tmp(iip1,llm)
186      REAL,SAVE :: dtetafi_tmp(iip1,llm)
187      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi_tmp
188      REAL,SAVE :: dpfi_tmp(iip1)
189
190      INTEGER :: true_itau
191      LOGICAL :: verbose=.true.
192      INTEGER :: iapptrac
193      INTEGER :: AdjustCount
194!      INTEGER :: var_time
195      LOGICAL :: ok_start_timer=.FALSE.
196      LOGICAL, SAVE :: firstcall=.TRUE.
197
198c$OMP MASTER
199      ItCount=0
200c$OMP END MASTER     
201      true_itau=0
202      FirstCaldyn=.TRUE.
203      FirstPhysic=.TRUE.
204      iapptrac=0
205      AdjustCount = 0
206      lafin=.false.
207     
208      itaufin   = nday*day_step
209      itaufinp1 = itaufin +1
210      modname="leapfrog_p"
211
212      itau = 0
213!      iday = day_ini+itau/day_step
214!      time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
215!         IF(time.GT.1.) THEN
216!          time = time-1.
217!          iday = iday+1
218!         ENDIF
219
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
231c-----------------------------------------------------------------------
232c   On initialise la pression et la fonction d'Exner :
233c   --------------------------------------------------
234
235c$OMP MASTER
236      dq=0.
237      CALL pression ( ip1jmp1, ap, bp, ps, p       )
238      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
239c$OMP END MASTER
240c-----------------------------------------------------------------------
241c   Debut de l'integration temporelle:
242c   ----------------------------------
243c et du parallelisme !!
244
245   1  CONTINUE
246
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))
250
251
252#ifdef CPP_IOIPSL
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
258      endif
259#endif
260
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
274c$OMP MASTER
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 )
283c$OMP END MASTER
284c$OMP BARRIER
285       else
286! Save fields obtained at previous time step as '...m1'
287         ijb=ij_begin
288         ije=ij_end
289
290c$OMP MASTER           
291         psm1     (ijb:ije) = ps    (ijb:ije)
292c$OMP END MASTER
293
294c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
295         DO l=1,llm     
296           ije=ij_end
297           ucovm1   (ijb:ije,l) = ucov  (ijb:ije,l)
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)
301                 
302           if (pole_sud) ije=ij_end-iip1
303           vcovm1(ijb:ije,l) = vcov  (ijb:ije,l)
304       
305
306         ENDDO
307c$OMP ENDDO 
308
309
310          CALL filtreg_p ( finvmaold ,jj_begin,jj_end,jjp1,
311     .                    llm, -2,2, .TRUE., 1 )
312
313       endif ! of if (FirstCaldyn)
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
324cym  ne sert a rien
325cym      call minmax(ijp1llm,q(:,:,3),zqmin,zqmax)
326
327   2  CONTINUE
328
329c$OMP MASTER
330      ItCount=ItCount+1
331      if (MOD(ItCount,1)==1) then
332        debug=.true.
333      else
334        debug=.false.
335      endif
336c$OMP END MASTER
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.
352c      idissip=1
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
357     s          .and. iflag_phys.EQ.1                 ) apphys = .TRUE.
358      ELSE
359         IF( MOD(itau   ,iconser) .EQ. 0              ) conser = .TRUE.
360         IF( MOD(itau+1,idissip)  .EQ. 0              ) apdiss = .TRUE.
361         IF( MOD(itau+1,iphysiq).EQ.0.AND.iflag_phys.EQ.1) apphys=.TRUE.
362      END IF
363
364cym    ---> Pour le moment     
365cym      apphys = .FALSE.
366      statcl = .FALSE.
367      conser = .FALSE.
368     
369      if (firstCaldyn) then
370c$OMP MASTER
371          call SetDistrib(jj_Nb_Caldyn)
372c$OMP END MASTER
373c$OMP BARRIER
374          firstCaldyn=.FALSE.
375cym          call InitTime
376c$OMP MASTER
377          call Init_timer
378c$OMP END MASTER
379      endif
380
381c$OMP MASTER     
382      IF (ok_start_timer) THEN
383        CALL InitTime
384        ok_start_timer=.FALSE.
385      ENDIF     
386c$OMP END MASTER     
387     
388      if (Adjust) then
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
393           AdjustCount=0
394           call allgather_timer_average
395
396        if (Verbose) then
397       
398        print *,'*********************************'
399        print *,'******    TIMER CALDYN     ******'
400        do i=0,mpi_size-1
401          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
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 *,'*********************************'
408        print *,'******    TIMER VANLEER    ******'
409        do i=0,mpi_size-1
410          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
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 *,'*********************************'
417        print *,'******    TIMER DISSIP    ******'
418        do i=0,mpi_size-1
419          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
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       
425        if (mpi_rank==0) call WriteBands
426       
427       endif
428       
429         call AdjustBands_caldyn
430         if (mpi_rank==0) call WriteBands
431         
432         call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
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)
460         call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
461     &                                jj_Nb_caldyn,0,0,TestRequest)
462         call Register_SwapFieldHallo(finvmaold,finvmaold,ip1jmp1,llm,
463     &                                jj_Nb_caldyn,0,0,TestRequest)
464 
465        do j=1,nqtot
466         call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
467     &                                jj_nb_caldyn,0,0,TestRequest)
468        enddo
469
470         call SetDistrib(jj_nb_caldyn)
471         call SendRequest(TestRequest)
472         call WaitRequest(TestRequest)
473         
474        call AdjustBands_dissip
475        call AdjustBands_physic
476
477      endif
478c$OMP END MASTER 
479      endif       
480     
481     
482     
483c-----------------------------------------------------------------------
484c   calcul des tendances dynamiques:
485c   --------------------------------
486c$OMP BARRIER
487c$OMP MASTER
488       call VTb(VThallo)
489c$OMP END MASTER
490
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       
500c       do j=1,nqtot
501c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
502c     *                       TestRequest)
503c        enddo
504
505       call SendRequest(TestRequest)
506c$OMP BARRIER
507       call WaitRequest(TestRequest)
508
509c$OMP MASTER
510       call VTe(VThallo)
511c$OMP END MASTER
512c$OMP BARRIER
513     
514      if (debug) then       
515!$OMP BARRIER
516!$OMP MASTER
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/)))
522        call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
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/)))
526        do j=1,nqtot
527          call WriteField_p('q'//trim(int2str(j)),
528     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
529        enddo
530!$OMP END MASTER       
531c$OMP BARRIER
532      endif
533
534     
535      True_itau=True_itau+1
536
537c$OMP MASTER
538      IF (prt_level>9) THEN
539        WRITE(lunout,*)"leapfrog_p: Iteration No",True_itau
540      ENDIF
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)
549c$OMP END MASTER
550!      var_time=time+iday-day_ini
551
552c$OMP BARRIER
553!      CALL FTRACE_REGION_BEGIN("caldyn")
554      time = jD_cur + jH_cur
555      CALL caldyn_p
556     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
557     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time )
558
559!      CALL FTRACE_REGION_END("caldyn")
560
561c$OMP MASTER
562      call VTe(VTcaldyn)
563c$OMP END MASTER     
564
565cc$OMP BARRIER
566cc$OMP MASTER
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/)))
577cc$OMP END MASTER
578
579c-----------------------------------------------------------------------
580c   calcul des tendances advection des traceurs (dont l'humidite)
581c   -------------------------------------------------------------
582
583      IF( forward. OR . leapf )  THEN
584cc$OMP PARALLEL DEFAULT(SHARED)
585c
586         CALL caladvtrac_p(q,pbaru,pbarv,
587     *        p, masse, dq,  teta,
588     .        flxw,pk, iapptrac)
589
590       IF (offline) THEN
591Cmaf stokage du flux de masse pour traceurs OFF-LINE
592
593#ifdef CPP_IOIPSL
594           CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
595     .   dtvr, itau)
596#endif
597
598
599         ENDIF ! of IF (offline)
600c
601      ENDIF ! of IF( forward. OR . leapf )
602cc$OMP END PARALLEL
603
604c-----------------------------------------------------------------------
605c   integrations dynamique et traceurs:
606c   ----------------------------------
607
608c$OMP MASTER
609       call VTb(VTintegre)
610c$OMP END MASTER
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/)))
619cc$OMP PARALLEL DEFAULT(SHARED)
620c$OMP BARRIER
621!       CALL FTRACE_REGION_BEGIN("integrd")
622
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
627!       CALL FTRACE_REGION_END("integrd")
628c$OMP BARRIER
629cc$OMP MASTER
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/)))
638c
639c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
640c      do j=1,nqtot
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
647
648
649c$OMP MASTER
650       call VTe(VTintegre)
651c$OMP END MASTER
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
664
665cc$OMP END PARALLEL
666
667c
668c
669       IF( apphys )  THEN
670c
671c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
672c
673cc$OMP PARALLEL DEFAULT(SHARED)
674cc$OMP+         PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
675
676c$OMP MASTER
677         call suspend_timer(timer_caldyn)
678
679        if (prt_level >= 10) then
680         write(lunout,*)
681     &   'leapfrog_p: Entree dans la physique : Iteration No ',true_itau
682        endif
683c$OMP END MASTER
684
685         CALL pression_p (  ip1jmp1, ap, bp, ps,  p      )
686
687c$OMP BARRIER
688         CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
689c$OMP BARRIER
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)
695
696c rajout debug
697c       lafin = .true.
698
699
700c   Inbterface avec les routines de phylmd (phymars ... )
701c   -----------------------------------------------------
702
703c+jld
704
705c  Diagnostique de conservation de l'energie : initialisation
706      IF (ip_ebil_dyn.ge.1 ) THEN
707          ztit='bil dyn'
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
713      ENDIF
714c-jld
715c$OMP BARRIER
716c$OMP MASTER
717        call VTb(VThallo)
718c$OMP END MASTER
719
720        call SetTag(Request_physic,800)
721       
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       
731        call Register_SwapFieldHallo(masse,masse,ip1jmp1,llm,
732     *                               jj_Nb_physic,1,2,Request_physic)
733
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       
749c        call SetDistrib(jj_nb_vanleer)
750        do j=1,nqtot
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
755
756        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
757     *                               jj_Nb_physic,2,2,Request_physic)
758       
759        call SendRequest(Request_Physic)
760c$OMP BARRIER
761        call WaitRequest(Request_Physic)       
762
763c$OMP BARRIER
764c$OMP MASTER
765        call SetDistrib(jj_nb_Physic)
766        call VTe(VThallo)
767       
768        call VTb(VTphysiq)
769c$OMP END MASTER
770c$OMP BARRIER
771
772cc$OMP MASTER       
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
780!        CALL FTRACE_REGION_BEGIN("calfis")
781        CALL calfis_p(lafin ,jD_cur, jH_cur,
782     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
783     $               du,dv,dteta,dq,
784     $               flxw,
785     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
786!        CALL FTRACE_REGION_END("calfis")
787        ijb=ij_begin
788        ije=ij_end 
789        if ( .not. pole_nord) then
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
800          dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim) 
801c$OMP END MASTER
802        endif ! of if ( .not. pole_nord)
803
804c$OMP BARRIER
805c$OMP MASTER
806        call SetDistrib(jj_nb_Physic_bis)
807
808        call VTb(VThallo)
809c$OMP END MASTER
810c$OMP BARRIER
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
824        do j=1,nqtot
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)
830c$OMP BARRIER
831        call WaitRequest(Request_Physic)
832             
833c$OMP BARRIER
834c$OMP MASTER
835        call VTe(VThallo)
836 
837        call SetDistrib(jj_nb_Physic)
838c$OMP END MASTER
839c$OMP BARRIER       
840                ijb=ij_begin
841        if (.not. pole_nord) then
842       
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
855          dpfi(ijb:ijb+iim)   = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
856c$OMP END MASTER
857         
858        endif ! of if (.not. pole_nord)
859c$OMP BARRIER
860cc$OMP MASTER       
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/)))
865cc$OMP END MASTER
866c     
867c      do j=1,nqtot
868c        call WriteField_p('dqfi'//trim(int2str(j)),
869c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
870c      enddo
871
872c      ajout des tendances physiques:
873c      ------------------------------
874         IF (ok_strato) THEN
875           CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
876         ENDIF
877       
878          CALL addfi_p( dtphys, leapf, forward   ,
879     $                  ucov, vcov, teta , q   ,ps ,
880     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
881
882c$OMP BARRIER
883c$OMP MASTER
884        call VTe(VTphysiq)
885
886        call VTb(VThallo)
887c$OMP END MASTER
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       
899        call Register_SwapField(masse,masse,ip1jmp1,llm,
900     *                               jj_Nb_caldyn,Request_physic)
901
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
917        do j=1,nqtot
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)
925c$OMP BARRIER
926        call WaitRequest(Request_Physic)     
927
928c$OMP BARRIER
929c$OMP MASTER
930       call VTe(VThallo)
931       call SetDistrib(jj_Nb_caldyn)
932c$OMP END MASTER
933c$OMP BARRIER
934c
935c  Diagnostique de conservation de l'energie : difference
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
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
950
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
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)
974c$OMP BARRIER
975       call WaitRequest(Request_Physic)     
976
977       call friction_p(ucov,vcov,iphysiq*dtvr)
978      ENDIF ! of IF(iflag_phys.EQ.2)
979
980
981        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
982c$OMP BARRIER
983        CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
984c$OMP BARRIER
985
986cc$OMP END PARALLEL
987
988c-----------------------------------------------------------------------
989c   dissipation horizontale et verticale  des petites echelles:
990c   ----------------------------------------------------------
991
992      IF(apdiss) THEN
993cc$OMP  PARALLEL DEFAULT(SHARED)
994cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
995c$OMP MASTER
996        call suspend_timer(timer_caldyn)
997       
998c       print*,'Entree dans la dissipation : Iteration No ',true_itau
999c   calcul de l'energie cinetique avant dissipation
1000c       print *,'Passage dans la dissipation'
1001
1002        call VTb(VThallo)
1003c$OMP END MASTER
1004
1005c$OMP BARRIER
1006
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)       
1023c$OMP BARRIER
1024        call WaitRequest(Request_dissip)       
1025
1026c$OMP BARRIER
1027c$OMP MASTER
1028        call SetDistrib(jj_Nb_dissip)
1029        call VTe(VThallo)
1030        call VTb(VTdissipation)
1031        call start_timer(timer_dissip)
1032c$OMP END MASTER
1033c$OMP BARRIER
1034
1035        call covcont_p(llm,ucov,vcov,ucont,vcont)
1036        call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1037
1038c   dissipation
1039
1040!        CALL FTRACE_REGION_BEGIN("dissip")
1041        CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1042!        CALL FTRACE_REGION_END("dissip")
1043         
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)
1049        ENDDO
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)
1055        ENDDO
1056c$OMP END DO NOWAIT       
1057
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
1065c$OMP BARRIER
1066c$OMP MASTER
1067            call suspend_timer(timer_dissip)
1068            call VTb(VThallo)
1069c$OMP END MASTER
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)
1073c$OMP BARRIER
1074            call WaitRequest(Request_Dissip)
1075c$OMP MASTER
1076            call VTe(VThallo)
1077            call resume_timer(timer_dissip)
1078c$OMP END MASTER
1079c$OMP BARRIER           
1080            call covcont_p(llm,ucov,vcov,ucont,vcont)
1081            call enercin_p(vcov,ucov,vcont,ucont,ecin)
1082           
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)
1090              enddo
1091            enddo
1092c$OMP END DO NOWAIT           
1093       endif
1094
1095       ijb=ij_begin
1096       ije=ij_end
1097c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)           
1098         do l=1,llm
1099           do ij=ijb,ije
1100              teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1101           enddo
1102         enddo
1103c$OMP END DO NOWAIT         
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
1115c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
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
1126c$OMP END DO NOWAIT
1127
1128c$OMP MASTER               
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
1137c$OMP END MASTER
1138        endif
1139       
1140        if (pole_sud) then
1141c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
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
1152c$OMP END DO NOWAIT
1153
1154c$OMP MASTER               
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
1163c$OMP END MASTER
1164        endif
1165
1166
1167c$OMP BARRIER
1168c$OMP MASTER
1169        call VTe(VTdissipation)
1170
1171        call stop_timer(timer_dissip)
1172       
1173        call VTb(VThallo)
1174c$OMP END MASTER
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)       
1191c$OMP BARRIER
1192        call WaitRequest(Request_dissip)       
1193
1194c$OMP BARRIER
1195c$OMP MASTER
1196        call SetDistrib(jj_Nb_caldyn)
1197        call VTe(VThallo)
1198        call resume_timer(timer_caldyn)
1199c        print *,'fin dissipation'
1200c$OMP END MASTER
1201c$OMP BARRIER
1202      END IF
1203
1204cc$OMP END PARALLEL
1205
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/)))
1221c$OMP MASTER     
1222      call stop_timer(timer_caldyn)
1223c$OMP END MASTER
1224      IF (itau==itaumax) then
1225c$OMP MASTER
1226            call allgather_timer_average
1227
1228      if (mpi_rank==0) then
1229       
1230        print *,'*********************************'
1231        print *,'******    TIMER CALDYN     ******'
1232        do i=0,mpi_size-1
1233          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
1234     &            '  : temps moyen :',
1235     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1236        enddo
1237     
1238        print *,'*********************************'
1239        print *,'******    TIMER VANLEER    ******'
1240        do i=0,mpi_size-1
1241          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
1242     &            '  : temps moyen :',
1243     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1244        enddo
1245     
1246        print *,'*********************************'
1247        print *,'******    TIMER DISSIP    ******'
1248        do i=0,mpi_size-1
1249          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
1250     &            '  : temps moyen :',
1251     &             timer_average(jj_nb_dissip(i),timer_dissip,i)
1252        enddo
1253       
1254        print *,'*********************************'
1255        print *,'******    TIMER PHYSIC    ******'
1256        do i=0,mpi_size-1
1257          print *,'proc',i,' :   Nb Bandes  :',jj_nb_physic(i),
1258     &            '  : temps moyen :',
1259     &             timer_average(jj_nb_physic(i),timer_physic,i)
1260        enddo
1261       
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()
1268      CALL print_filtre_timer
1269      call fin_getparam
1270        call finalize_parallel
1271c$OMP END MASTER
1272c$OMP BARRIER
1273        RETURN
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
1283!              iday= day_ini+itau/day_step
1284!              time= REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1285!                IF(time.GT.1.) THEN
1286!                  time = time-1.
1287!                  iday = iday+1
1288!                ENDIF
1289            ENDIF
1290
1291
1292            IF( itau. EQ. itaufinp1 ) then
1293
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
1305c$OMP MASTER
1306              call fin_getparam
1307              call finalize_parallel
1308c$OMP END MASTER
1309              abort_message = 'Simulation finished'
1310              call abort_gcm(modname,abort_message,0)
1311              RETURN
1312            ENDIF
1313c-----------------------------------------------------------------------
1314c   ecriture du fichier histoire moyenne:
1315c   -------------------------------------
1316
1317            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1318c$OMP BARRIER
1319               IF(itau.EQ.itaufin) THEN
1320                  iav=1
1321               ELSE
1322                  iav=0
1323               ENDIF
1324#ifdef CPP_IOIPSL
1325             IF (ok_dynzon) THEN
1326             call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1327             call SendRequest(TestRequest)
1328c$OMP BARRIER
1329              call WaitRequest(TestRequest)
1330c$OMP BARRIER
1331c$OMP MASTER
1332!              CALL writedynav_p(histaveid, itau,vcov ,
1333!     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
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)
1338c$OMP END MASTER
1339              ENDIF !ok_dynzon
1340#endif
1341            ENDIF
1342
1343c-----------------------------------------------------------------------
1344c   ecriture de la bande histoire:
1345c   ------------------------------
1346
1347c      IF( MOD(itau,iecri         ).EQ.0) THEN
1348
1349            IF( MOD(itau,iecri*day_step).EQ.0) THEN
1350c$OMP BARRIER
1351c$OMP MASTER
1352              nbetat = nbetatdem
1353              CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1354       
1355cym        unat=0.
1356       
1357              ijb=ij_begin
1358              ije=ij_end
1359       
1360              if (pole_nord) then
1361                ijb=ij_begin+iip1
1362                unat(1:iip1,:)=0.
1363              endif
1364       
1365              if (pole_sud) then
1366                ije=ij_end-iip1
1367                unat(ij_end-iip1+1:ij_end,:)=0.
1368              endif
1369           
1370              do l=1,llm
1371                unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1372              enddo
1373
1374              ijb=ij_begin
1375              ije=ij_end
1376              if (pole_sud) ije=ij_end-iip1
1377       
1378              do l=1,llm
1379                vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1380              enddo
1381       
1382#ifdef CPP_IOIPSL
1383 
1384!              CALL writehist_p(histid,histvid, itau,vcov,
1385!     &                         ucov,teta,phi,q,masse,ps,phis)
1386
1387#endif
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)
1402c$OMP END MASTER
1403            ENDIF ! of IF(MOD(itau,iecri).EQ.0)
1404
1405            IF(itau.EQ.itaufin) THEN
1406
1407c$OMP BARRIER
1408c$OMP MASTER
1409
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")
1415
1416!              CLOSE(99)
1417c$OMP END MASTER
1418            ENDIF ! of IF (itau.EQ.itaufin)
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
1449            END IF ! of IF (MOD(itau,iperiod).EQ.0)
1450                   !    ELSEIF (MOD(itau-1,iperiod).EQ.0)
1451
1452
1453      ELSE ! of IF (.not.purmats)
1454
1455c       ........................................................
1456c       ..............       schema  matsuno        ...............
1457c       ........................................................
1458            IF( forward )  THEN
1459
1460             itau =  itau + 1
1461!             iday = day_ini+itau/day_step
1462!             time = REAL(itau-(iday-day_ini)*day_step)/day_step+time_0
1463!
1464!                  IF(time.GT.1.) THEN
1465!                   time = time-1.
1466!                   iday = iday+1
1467!                  ENDIF
1468
1469               forward =  .FALSE.
1470               IF( itau. EQ. itaufinp1 ) then 
1471c$OMP MASTER
1472                 call fin_getparam
1473                 call finalize_parallel
1474c$OMP END MASTER
1475                 abort_message = 'Simulation finished'
1476                 call abort_gcm(modname,abort_message,0)
1477                 RETURN
1478               ENDIF
1479               GO TO 2
1480
1481            ELSE ! of IF(forward)
1482
1483              IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1484               IF(itau.EQ.itaufin) THEN
1485                  iav=1
1486               ELSE
1487                  iav=0
1488               ENDIF
1489#ifdef CPP_IOIPSL
1490               IF (ok_dynzon) THEN
1491c$OMP BARRIER
1492
1493               call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1494               call SendRequest(TestRequest)
1495c$OMP BARRIER
1496               call WaitRequest(TestRequest)
1497
1498c$OMP BARRIER
1499c$OMP MASTER
1500!               CALL writedynav_p(histaveid, itau,vcov ,
1501!     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1502               CALL bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1503     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1504c$OMP END MASTER
1505               END IF !ok_dynzon
1506#endif
1507              ENDIF ! of IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin)
1508
1509
1510c               IF(MOD(itau,iecri         ).EQ.0) THEN
1511              IF(MOD(itau,iecri*day_step).EQ.0) THEN
1512c$OMP BARRIER
1513c$OMP MASTER
1514                nbetat = nbetatdem
1515                CALL geopot_p(ip1jmp1,teta,pk,pks,phis,phi)
1516
1517cym        unat=0.
1518                ijb=ij_begin
1519                ije=ij_end
1520       
1521                if (pole_nord) then
1522                  ijb=ij_begin+iip1
1523                  unat(1:iip1,:)=0.
1524                endif
1525       
1526                if (pole_sud) then
1527                  ije=ij_end-iip1
1528                  unat(ij_end-iip1+1:ij_end,:)=0.
1529                endif
1530           
1531                do l=1,llm
1532                  unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1533                enddo
1534
1535                ijb=ij_begin
1536                ije=ij_end
1537                if (pole_sud) ije=ij_end-iip1
1538       
1539                do l=1,llm
1540                  vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1541                enddo
1542
1543#ifdef CPP_IOIPSL
1544
1545!                CALL writehist_p(histid, histvid, itau,vcov ,
1546!     &                           ucov,teta,phi,q,masse,ps,phis)
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
1557c     
1558                  if (mpi_rank==0) then
1559#include "write_grads_dyn.h"
1560                  endif
1561                endif ! of if (output_grads_dyn)
1562
1563c$OMP END MASTER
1564              ENDIF ! of IF(MOD(itau,iecri*day_step).EQ.0)
1565
1566              IF(itau.EQ.itaufin) THEN
1567                if (planet_type.eq."earth") then
1568c$OMP MASTER
1569                   CALL dynredem1_p("restart.nc",0.0,
1570     .                               vcov,ucov,teta,q,masse,ps)
1571c$OMP END MASTER
1572                endif ! of if (planet_type.eq."earth")
1573              ENDIF ! of IF(itau.EQ.itaufin)
1574
1575              forward = .TRUE.
1576              GO TO  1
1577
1578            ENDIF ! of IF (forward)
1579
1580      END IF ! of IF(.not.purmats)
1581c$OMP MASTER
1582      call fin_getparam
1583      call finalize_parallel
1584c$OMP END MASTER
1585      RETURN
1586      END
Note: See TracBrowser for help on using the repository browser.