source: LMDZ4/branches/LMDZ4_V3_patches/libf/dyn3dpar/leapfrog_p.F @ 853

Last change on this file since 853 was 774, checked in by Laurent Fairhead, 17 years ago

Suite du merge entre la version et la HEAD: quelques modifications de
Yann sur le

LF

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