source: LMDZ4/tags/LMDZ4_V3_3/libf/dyn3dpar/leapfrog_p.F @ 919

Last change on this file since 919 was 919, checked in by (none), 16 years ago

This commit was manufactured by cvs2svn to create tag 'LMDZ4_V3_3'.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 44.8 KB
Line 
1!
2! $Header$
3!
4c
5c
6#define IO_DEBUG
7
8c#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)
483c$OMP END MASTER
484
485       call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,TestRequest)
486       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,TestRequest)
487       call Register_Hallo(teta,ip1jmp1,llm,1,1,1,1,TestRequest)
488       call Register_Hallo(ps,ip1jmp1,1,1,2,2,1,TestRequest)
489       call Register_Hallo(pkf,ip1jmp1,llm,1,1,1,1,TestRequest)
490       call Register_Hallo(pk,ip1jmp1,llm,1,1,1,1,TestRequest)
491       call Register_Hallo(pks,ip1jmp1,1,1,1,1,1,TestRequest)
492       call Register_Hallo(p,ip1jmp1,llmp1,1,1,1,1,TestRequest)
493       
494c       do j=1,nqmx
495c         call Register_Hallo(q(1,1,j),ip1jmp1,llm,1,1,1,1,
496c     *                       TestRequest)
497c        enddo
498
499       call SendRequest(TestRequest)
500c$OMP BARRIER
501       call WaitRequest(TestRequest)
502
503c$OMP MASTER
504       call VTe(VThallo)
505c$OMP END MASTER
506c$OMP BARRIER
507     
508      if (debug) then   
509        call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
510        call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
511        call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
512        call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
513        call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
514        call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
515        call WriteField_p('pks',reshape(pks,(/iip1,jmp1/)))
516        call WriteField_p('pkf',reshape(pkf,(/iip1,jmp1,llm/)))
517        call WriteField_p('phis',reshape(phis,(/iip1,jmp1/)))
518        do j=1,nqmx
519          call WriteField_p('q'//trim(int2str(j)),
520     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
521        enddo       
522c$OMP BARRIER
523      endif
524
525     
526      True_itau=True_itau+1
527
528c$OMP MASTER
529      print*,"Iteration No",True_itau
530
531
532      call start_timer(timer_caldyn)
533
534      CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi   )
535
536     
537      call VTb(VTcaldyn)
538c$OMP END MASTER
539      var_time=time+iday-day_ini
540
541c$OMP BARRIER
542!      CALL FTRACE_REGION_BEGIN("caldyn")
543      CALL caldyn_p
544     $  ( itau,ucov,vcov,teta,ps,masse,pk,pkf,phis ,
545     $    phi,conser,du,dv,dteta,dp,w, pbaru,pbarv, time+iday-day_ini )
546
547!      CALL FTRACE_REGION_END("caldyn")
548c$OMP MASTER
549      call VTe(VTcaldyn)
550c$OMP END MASTER     
551
552cc$OMP BARRIER
553cc$OMP MASTER
554c      call WriteField_p('du',reshape(du,(/iip1,jmp1,llm/)))
555c      call WriteField_p('dv',reshape(dv,(/iip1,jjm,llm/)))
556c      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
557c      call WriteField_p('dp',reshape(dp,(/iip1,jmp1/)))
558c      call WriteField_p('w',reshape(w,(/iip1,jmp1,llm/)))
559c      call WriteField_p('pbaru',reshape(pbaru,(/iip1,jmp1,llm/)))
560c      call WriteField_p('pbarv',reshape(pbarv,(/iip1,jjm,llm/)))
561c      call WriteField_p('p',reshape(p,(/iip1,jmp1,llmp1/)))
562c      call WriteField_p('masse',reshape(masse,(/iip1,jmp1,llm/)))
563c      call WriteField_p('pk',reshape(pk,(/iip1,jmp1,llm/)))
564cc$OMP END MASTER
565
566c-----------------------------------------------------------------------
567c   calcul des tendances advection des traceurs (dont l'humidite)
568c   -------------------------------------------------------------
569
570      IF( forward. OR . leapf )  THEN
571cc$OMP PARALLEL DEFAULT(SHARED)
572c
573#ifdef INCA
574             CALL caladvtrac_p(q,pbaru,pbarv,
575     *                      p, masse, dq,  teta,
576     .             flxw,
577     .             pk,
578     .             iapptrac)
579#else
580             CALL caladvtrac_p(q,pbaru,pbarv,
581     *                      p, masse, dq,  teta,
582     .             pk,iapptrac)
583#endif
584
585
586
587c      do j=1,nqmx
588c        call WriteField_p('q'//trim(int2str(j)),
589c     .                reshape(q(:,:,j),(/iip1,jmp1,llm/)))
590c        call WriteField_p('dq'//trim(int2str(j)),
591c     .                reshape(dq(:,:,j),(/iip1,jmp1,llm/)))
592c      enddo
593       IF (offline) THEN
594Cmaf stokage du flux de masse pour traceurs OFF-LINE
595
596#ifdef CPP_IOIPSL
597           CALL fluxstokenc_p(pbaru,pbarv,masse,teta,phi,phis,
598     .   dtvr, itau)
599#endif
600
601
602         ENDIF
603c
604      ENDIF
605cc$OMP END PARALLEL
606
607c-----------------------------------------------------------------------
608c   integrations dynamique et traceurs:
609c   ----------------------------------
610
611c$OMP MASTER
612       call VTb(VTintegre)
613c$OMP END MASTER
614c      call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
615c      call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
616c      call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
617c      call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
618c      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
619c      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
620c      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
621c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
622cc$OMP PARALLEL DEFAULT(SHARED)
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
628c      call WriteField_p('ucovm1',reshape(ucovm1,(/iip1,jmp1,llm/)))
629c      call WriteField_p('vcovm1',reshape(vcovm1,(/iip1,jjm,llm/)))
630c      call WriteField_p('tetam1',reshape(tetam1,(/iip1,jmp1,llm/)))
631c      call WriteField_p('psm1',reshape(psm1,(/iip1,jmp1/)))
632c      call WriteField_p('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
633c      call WriteField_p('vcov',reshape(vcov,(/iip1,jjm,llm/)))
634c      call WriteField_p('teta',reshape(teta,(/iip1,jmp1,llm/)))
635c      call WriteField_p('dteta',reshape(dteta,(/iip1,jmp1,llm/)))
636
637c      call WriteField_p('ps',reshape(ps,(/iip1,jmp1/)))
638c$OMP MASTER
639       call VTe(VTintegre)
640c$OMP END MASTER
641c .P.Le Van (26/04/94  ajout de  finvpold dans l'appel d'integrd)
642c
643c-----------------------------------------------------------------------
644c   calcul des tendances physiques:
645c   -------------------------------
646c    ########   P.Le Van ( Modif le  6/02/95 )   ###########
647c
648       IF( purmats )  THEN
649          IF( itau.EQ.itaufin.AND..NOT.forward ) lafin = .TRUE.
650       ELSE
651          IF( itau+1. EQ. itaufin )              lafin = .TRUE.
652       ENDIF
653
654cc$OMP END PARALLEL
655
656c
657c
658       IF( apphys )  THEN
659c
660c     .......   Ajout   P.Le Van ( 17/04/96 )   ...........
661c
662cc$OMP PARALLEL DEFAULT(SHARED)
663cc$OMP+         PRIVATE(rdaym_ini,rdayvrai,ijb,ije)
664
665c$OMP MASTER
666         call suspend_timer(timer_caldyn)
667         print*,'Entree dans la physique : Iteration No ',true_itau
668c$OMP END MASTER
669
670         CALL pression_p (  ip1jmp1, ap, bp, ps,  p      )
671
672c$OMP BARRIER
673
674         CALL exner_hyb_p(  ip1jmp1, ps, p,alpha,beta,pks, pk, pkf )
675c$OMP BARRIER
676           rdaym_ini  = itau * dtvr / daysec
677           rdayvrai   = rdaym_ini  + day_ini
678
679
680c rajout debug
681c       lafin = .true.
682
683
684c   Inbterface avec les routines de phylmd (phymars ... )
685c   -----------------------------------------------------
686
687#ifdef CPP_PHYS
688c+jld
689
690c  Diagnostique de conservation de l'energie : initialisation
691      IF (ip_ebil_dyn.ge.1 ) THEN
692          ztit='bil dyn'
693          CALL diagedyn(ztit,2,1,1,dtphys
694     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
695      ENDIF
696c-jld
697c$OMP BARRIER
698c$OMP MASTER
699        call VTb(VThallo)
700c$OMP END MASTER
701
702        call SetTag(Request_physic,800)
703       
704        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
705     *                               jj_Nb_physic,2,2,Request_physic)
706       
707        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
708     *                               jj_Nb_physic,2,2,Request_physic)
709       
710        call Register_SwapFieldHallo(teta,teta,ip1jmp1,llm,
711     *                               jj_Nb_physic,2,2,Request_physic)
712       
713        call Register_SwapFieldHallo(p,p,ip1jmp1,llmp1,
714     *                               jj_Nb_physic,2,2,Request_physic)
715       
716        call Register_SwapFieldHallo(pk,pk,ip1jmp1,llm,
717     *                               jj_Nb_physic,2,2,Request_physic)
718       
719        call Register_SwapFieldHallo(phis,phis,ip1jmp1,1,
720     *                               jj_Nb_physic,2,2,Request_physic)
721       
722        call Register_SwapFieldHallo(phi,phi,ip1jmp1,llm,
723     *                               jj_Nb_physic,2,2,Request_physic)
724       
725        call Register_SwapFieldHallo(w,w,ip1jmp1,llm,
726     *                               jj_Nb_physic,2,2,Request_physic)
727       
728c       call SetDistrib(jj_nb_vanleer)
729        do j=1,nqmx
730 
731          call Register_SwapFieldHallo(q(1,1,j),q(1,1,j),ip1jmp1,llm,
732     *                               jj_Nb_physic,2,2,Request_physic)
733        enddo
734#ifdef INCA
735        call Register_SwapFieldHallo(flxw,flxw,ip1jmp1,llm,
736     *                               jj_Nb_physic,2,2,Request_physic)
737#endif
738       
739        call SendRequest(Request_Physic)
740c$OMP BARRIER
741        call WaitRequest(Request_Physic)       
742
743c$OMP BARRIER
744c$OMP MASTER
745        call SetDistrib(jj_nb_Physic)
746        call VTe(VThallo)
747        call VTb(VTphysiq)
748c$OMP END MASTER
749c$OMP BARRIER
750
751cc$OMP MASTER   
752c      call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
753c      call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
754c      call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
755c      call WriteField_p('pfi',reshape(p,(/iip1,jmp1,llmp1/)))
756c      call WriteField_p('pkfi',reshape(pk,(/iip1,jmp1,llm/)))
757cc$OMP END MASTER
758cc$OMP BARRIER
759       
760        CALL calfis_p( nq, lafin ,rdayvrai,time  ,
761     $               ucov,vcov,teta,q,masse,ps,p,pk,phis,phi ,
762     $               du,dv,dteta,dq,w,
763#ifdef INCA
764     $               flxw,
765#endif
766     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
767        ijb=ij_begin
768        ije=ij_end 
769        if ( .not. pole_nord) then
770c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
771          DO l=1,llm
772          dufi_tmp(1:iip1,l)   = dufi(ijb:ijb+iim,l)
773          dvfi_tmp(1:iip1,l)   = dvfi(ijb:ijb+iim,l) 
774          dtetafi_tmp(1:iip1,l)= dtetafi(ijb:ijb+iim,l) 
775          dqfi_tmp(1:iip1,l,:) = dqfi(ijb:ijb+iim,l,:) 
776          ENDDO
777c$OMP END DO NOWAIT
778
779c$OMP MASTER
780          dpfi_tmp(1:iip1)     = dpfi(ijb:ijb+iim) 
781c$OMP END MASTER
782        endif
783
784c$OMP BARRIER
785c$OMP MASTER
786        call SetDistrib(jj_nb_Physic_bis)
787        call VTb(VThallo)
788c$OMP END MASTER
789c$OMP BARRIER
790 
791        call Register_Hallo(dufi,ip1jmp1,llm,
792     *                      1,0,0,1,Request_physic)
793       
794        call Register_Hallo(dvfi,ip1jm,llm,
795     *                      1,0,0,1,Request_physic)
796       
797        call Register_Hallo(dtetafi,ip1jmp1,llm,
798     *                      1,0,0,1,Request_physic)
799
800        call Register_Hallo(dpfi,ip1jmp1,1,
801     *                      1,0,0,1,Request_physic)
802
803        do j=1,nqmx
804          call Register_Hallo(dqfi(1,1,j),ip1jmp1,llm,
805     *                        1,0,0,1,Request_physic)
806        enddo
807       
808        call SendRequest(Request_Physic)
809c$OMP BARRIER
810        call WaitRequest(Request_Physic)
811             
812c$OMP BARRIER
813c$OMP MASTER
814        call VTe(VThallo)
815        call SetDistrib(jj_nb_Physic)
816c$OMP END MASTER
817c$OMP BARRIER   
818       
819        ijb=ij_begin
820        if (.not. pole_nord) then
821       
822c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
823          DO l=1,llm
824            dufi(ijb:ijb+iim,l) = dufi(ijb:ijb+iim,l)+dufi_tmp(1:iip1,l)
825            dvfi(ijb:ijb+iim,l) = dvfi(ijb:ijb+iim,l)+dvfi_tmp(1:iip1,l)
826            dtetafi(ijb:ijb+iim,l) = dtetafi(ijb:ijb+iim,l)
827     &                              +dtetafi_tmp(1:iip1,l)
828            dqfi(ijb:ijb+iim,l,:) = dqfi(ijb:ijb+iim,l,:)
829     &                              + dqfi_tmp(1:iip1,l,:)
830          ENDDO
831c$OMP END DO NOWAIT
832
833c$OMP MASTER
834          dpfi(ijb:ijb+iim)   = dpfi(ijb:ijb+iim)+ dpfi_tmp(1:iip1)
835c$OMP END MASTER
836         
837        endif
838c$OMP BARRIER
839cc$OMP MASTER   
840c      call WriteField_p('dufi',reshape(dufi,(/iip1,jmp1,llm/)))
841c      call WriteField_p('dvfi',reshape(dvfi,(/iip1,jjm,llm/)))
842c      call WriteField_p('dtetafi',reshape(dtetafi,(/iip1,jmp1,llm/)))
843c      call WriteField_p('dpfi',reshape(dpfi,(/iip1,jmp1/)))
844cc$OMP END MASTER
845c     
846c      do j=1,nqmx
847c        call WriteField_p('dqfi'//trim(int2str(j)),
848c     .                reshape(dqfi(:,:,j),(/iip1,jmp1,llm/)))
849c      enddo
850
851c      ajout des tendances physiques:
852c      ------------------------------
853          CALL addfi_p( nqmx, dtphys, leapf, forward   ,
854     $                  ucov, vcov, teta , q   ,ps ,
855     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
856
857c$OMP BARRIER
858c$OMP MASTER
859        call VTe(VTphysiq)
860        call VTb(VThallo)
861c$OMP END MASTER
862
863        call SetTag(Request_physic,800)
864        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
865     *                               jj_Nb_caldyn,Request_physic)
866       
867        call Register_SwapField(vcov,vcov,ip1jm,llm,
868     *                               jj_Nb_caldyn,Request_physic)
869       
870        call Register_SwapField(teta,teta,ip1jmp1,llm,
871     *                               jj_Nb_caldyn,Request_physic)
872       
873        call Register_SwapField(p,p,ip1jmp1,llmp1,
874     *                               jj_Nb_caldyn,Request_physic)
875       
876        call Register_SwapField(pk,pk,ip1jmp1,llm,
877     *                               jj_Nb_caldyn,Request_physic)
878       
879        call Register_SwapField(phis,phis,ip1jmp1,1,
880     *                               jj_Nb_caldyn,Request_physic)
881       
882        call Register_SwapField(phi,phi,ip1jmp1,llm,
883     *                               jj_Nb_caldyn,Request_physic)
884       
885        call Register_SwapField(w,w,ip1jmp1,llm,
886     *                               jj_Nb_caldyn,Request_physic)
887
888        do j=1,nqmx
889       
890          call Register_SwapField(q(1,1,j),q(1,1,j),ip1jmp1,llm,
891     *                               jj_Nb_caldyn,Request_physic)
892       
893        enddo
894
895        call SendRequest(Request_Physic)
896c$OMP BARRIER
897        call WaitRequest(Request_Physic)     
898
899c$OMP BARRIER
900c$OMP MASTER
901       call VTe(VThallo)
902       call SetDistrib(jj_Nb_caldyn)
903c$OMP END MASTER
904c$OMP BARRIER
905c
906c  Diagnostique de conservation de l'energie : difference
907      IF (ip_ebil_dyn.ge.1 ) THEN
908          ztit='bil phys'
909          CALL diagedyn(ztit,2,1,1,dtphys
910     e  , ucov    , vcov , ps, p ,pk , teta , q(:,:,1), q(:,:,2))
911      ENDIF
912
913cc$OMP MASTER     
914c      if (debug) then
915c       call WriteField_p('ucovfi',reshape(ucov,(/iip1,jmp1,llm/)))
916c       call WriteField_p('vcovfi',reshape(vcov,(/iip1,jjm,llm/)))
917c       call WriteField_p('tetafi',reshape(teta,(/iip1,jmp1,llm/)))
918c      endif
919cc$OMP END MASTER
920
921#else
922
923c   Calcul academique de la physique = Rappel Newtonien + fritcion
924c   --------------------------------------------------------------
925cym       teta(:,:)=teta(:,:)
926cym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
927       ijb=ij_begin
928       ije=ij_end
929       teta(ijb:ije,:)=teta(ijb:ije,:)
930     s  -iphysiq*dtvr*(teta(ijb:ije,:)-tetarappel(ijb:ije,:))/taurappel
931
932       call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
933       call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Physic)
934       call SendRequest(Request_Physic)
935       call WaitRequest(Request_Physic)     
936
937       call friction_p(ucov,vcov,iphysiq*dtvr)
938
939#endif
940
941c-jld
942c$OMP MASTER
943         call resume_timer(timer_caldyn)
944         if (FirstPhysic) then
945           ok_start_timer=.TRUE.
946           FirstPhysic=.false.
947         endif
948c$OMP END MASTER
949       ENDIF
950
951        CALL pression_p ( ip1jmp1, ap, bp, ps, p                  )
952c$OMP BARRIER
953
954
955        CALL exner_hyb_p( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
956c$OMP BARRIER
957
958cc$OMP END PARALLEL
959
960c-----------------------------------------------------------------------
961c   dissipation horizontale et verticale  des petites echelles:
962c   ----------------------------------------------------------
963
964      IF(apdiss) THEN
965cc$OMP  PARALLEL DEFAULT(SHARED)
966cc$OMP+          PRIVATE(ijb,ije,tppn,tpn,tpps,tps)
967c$OMP MASTER
968        call suspend_timer(timer_caldyn)
969       
970        print*,'Entree dans la dissipation : Iteration No ',true_itau
971c   calcul de l'energie cinetique avant dissipation
972        print *,'Passage dans la dissipation'
973
974        call VTb(VThallo)
975c$OMP END MASTER
976
977c$OMP BARRIER
978
979        call Register_SwapFieldHallo(ucov,ucov,ip1jmp1,llm,
980     *                          jj_Nb_dissip,1,1,Request_dissip)
981
982        call Register_SwapFieldHallo(vcov,vcov,ip1jm,llm,
983     *                          jj_Nb_dissip,1,1,Request_dissip)
984
985        call Register_SwapField(teta,teta,ip1jmp1,llm,
986     *                          jj_Nb_dissip,Request_dissip)
987
988        call Register_SwapField(p,p,ip1jmp1,llmp1,
989     *                          jj_Nb_dissip,Request_dissip)
990
991        call Register_SwapField(pk,pk,ip1jmp1,llm,
992     *                          jj_Nb_dissip,Request_dissip)
993
994        call SendRequest(Request_dissip)       
995c$OMP BARRIER
996        call WaitRequest(Request_dissip)       
997
998c$OMP BARRIER
999c$OMP MASTER
1000        call SetDistrib(jj_Nb_dissip)
1001        call VTe(VThallo)
1002        call VTb(VTdissipation)
1003        call start_timer(timer_dissip)
1004c$OMP END MASTER
1005c$OMP BARRIER
1006
1007        call covcont_p(llm,ucov,vcov,ucont,vcont)
1008        call enercin_p(vcov,ucov,vcont,ucont,ecin0)
1009
1010c   dissipation
1011
1012        CALL dissip_p(vcov,ucov,teta,p,dvdis,dudis,dtetadis)
1013         
1014        ijb=ij_begin
1015        ije=ij_end
1016c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1017        DO l=1,llm
1018          ucov(ijb:ije,l)=ucov(ijb:ije,l)+dudis(ijb:ije,l)
1019        ENDDO
1020c$OMP END DO NOWAIT     
1021        if (pole_sud) ije=ije-iip1
1022c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
1023        DO l=1,llm
1024          vcov(ijb:ije,l)=vcov(ijb:ije,l)+dvdis(ijb:ije,l)
1025        ENDDO
1026c$OMP END DO NOWAIT     
1027
1028c       teta=teta+dtetadis
1029
1030
1031c------------------------------------------------------------------------
1032        if (dissip_conservative) then
1033C       On rajoute la tendance due a la transform. Ec -> E therm. cree
1034C       lors de la dissipation
1035c$OMP BARRIER
1036c$OMP MASTER
1037            call suspend_timer(timer_dissip)
1038            call VTb(VThallo)
1039c$OMP END MASTER
1040            call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Request_Dissip)
1041            call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Request_Dissip)
1042            call SendRequest(Request_Dissip)
1043c$OMP BARRIER
1044            call WaitRequest(Request_Dissip)
1045c$OMP MASTER
1046            call VTe(VThallo)
1047            call resume_timer(timer_dissip)
1048c$OMP END MASTER
1049c$OMP BARRIER       
1050            call covcont_p(llm,ucov,vcov,ucont,vcont)
1051            call enercin_p(vcov,ucov,vcont,ucont,ecin)
1052           
1053            ijb=ij_begin
1054            ije=ij_end
1055c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
1056            do l=1,llm
1057              do ij=ijb,ije
1058                dtetaecdt(ij,l)= (ecin0(ij,l)-ecin(ij,l))/ pk(ij,l)
1059                dtetadis(ij,l)=dtetadis(ij,l)+dtetaecdt(ij,l)
1060              enddo
1061            enddo
1062c$OMP END DO NOWAIT         
1063       endif
1064
1065       ijb=ij_begin
1066       ije=ij_end
1067c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)         
1068         do l=1,llm
1069           do ij=ijb,ije
1070              teta(ij,l)=teta(ij,l)+dtetadis(ij,l)
1071           enddo
1072         enddo
1073c$OMP END DO NOWAIT     
1074c------------------------------------------------------------------------
1075
1076
1077c    .......        P. Le Van (  ajout  le 17/04/96  )   ...........
1078c   ...      Calcul de la valeur moyenne, unique de h aux poles  .....
1079c
1080
1081        ijb=ij_begin
1082        ije=ij_end
1083         
1084        if (pole_nord) then
1085c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1086          DO l  =  1, llm
1087            DO ij =  1,iim
1088             tppn(ij)  = aire(  ij    ) * teta(  ij    ,l)
1089            ENDDO
1090             tpn  = SSUM(iim,tppn,1)/apoln
1091
1092            DO ij = 1, iip1
1093             teta(  ij    ,l) = tpn
1094            ENDDO
1095          ENDDO
1096c$OMP END DO NOWAIT
1097
1098c$OMP MASTER               
1099          DO ij =  1,iim
1100            tppn(ij)  = aire(  ij    ) * ps (  ij    )
1101          ENDDO
1102            tpn  = SSUM(iim,tppn,1)/apoln
1103 
1104          DO ij = 1, iip1
1105            ps(  ij    ) = tpn
1106          ENDDO
1107c$OMP END MASTER
1108        endif
1109       
1110        if (pole_sud) then
1111c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
1112          DO l  =  1, llm
1113            DO ij =  1,iim
1114             tpps(ij)  = aire(ij+ip1jm) * teta(ij+ip1jm,l)
1115            ENDDO
1116             tps  = SSUM(iim,tpps,1)/apols
1117
1118            DO ij = 1, iip1
1119             teta(ij+ip1jm,l) = tps
1120            ENDDO
1121          ENDDO
1122c$OMP END DO NOWAIT
1123
1124c$OMP MASTER               
1125          DO ij =  1,iim
1126            tpps(ij)  = aire(ij+ip1jm) * ps (ij+ip1jm)
1127          ENDDO
1128            tps  = SSUM(iim,tpps,1)/apols
1129 
1130          DO ij = 1, iip1
1131            ps(ij+ip1jm) = tps
1132          ENDDO
1133c$OMP END MASTER
1134        endif
1135
1136
1137c$OMP BARRIER
1138c$OMP MASTER
1139        call VTe(VTdissipation)
1140        call stop_timer(timer_dissip)
1141        call VTb(VThallo)
1142c$OMP END MASTER
1143        call Register_SwapField(ucov,ucov,ip1jmp1,llm,
1144     *                          jj_Nb_caldyn,Request_dissip)
1145
1146        call Register_SwapField(vcov,vcov,ip1jm,llm,
1147     *                          jj_Nb_caldyn,Request_dissip)
1148
1149        call Register_SwapField(teta,teta,ip1jmp1,llm,
1150     *                          jj_Nb_caldyn,Request_dissip)
1151
1152        call Register_SwapField(p,p,ip1jmp1,llmp1,
1153     *                          jj_Nb_caldyn,Request_dissip)
1154
1155        call Register_SwapField(pk,pk,ip1jmp1,llm,
1156     *                          jj_Nb_caldyn,Request_dissip)
1157
1158        call SendRequest(Request_dissip)       
1159c$OMP BARRIER
1160        call WaitRequest(Request_dissip)       
1161
1162c$OMP BARRIER
1163c$OMP MASTER
1164        call SetDistrib(jj_Nb_caldyn)
1165        call VTe(VThallo)
1166        call resume_timer(timer_caldyn)
1167        print *,'fin dissipation'
1168c$OMP END MASTER
1169c$OMP BARRIER
1170      END IF
1171
1172cc$OMP END PARALLEL
1173
1174c ajout debug
1175c              IF( lafin ) then 
1176c                abort_message = 'Simulation finished'
1177c                call abort_gcm(modname,abort_message,0)
1178c              ENDIF
1179       
1180c   ********************************************************************
1181c   ********************************************************************
1182c   .... fin de l'integration dynamique  et physique pour le pas itau ..
1183c   ********************************************************************
1184c   ********************************************************************
1185
1186c   preparation du pas d'integration suivant  ......
1187cym      call WriteField('ucov',reshape(ucov,(/iip1,jmp1,llm/)))
1188cym      call WriteField('vcov',reshape(vcov,(/iip1,jjm,llm/)))
1189c$OMP MASTER     
1190      call stop_timer(timer_caldyn)
1191c$OMP END MASTER
1192      IF (itau==itaumax) then
1193c$OMP MASTER
1194            call allgather_timer_average
1195
1196      if (mpi_rank==0) then
1197       
1198        print *,'*********************************'
1199        print *,'******    TIMER CALDYN     ******'
1200        do i=0,mpi_size-1
1201          print *,'proc',i,' :   Nb Bandes  :',jj_nb_caldyn(i),
1202     &            '  : temps moyen :',
1203     &             timer_average(jj_nb_caldyn(i),timer_caldyn,i)
1204        enddo
1205     
1206        print *,'*********************************'
1207        print *,'******    TIMER VANLEER    ******'
1208        do i=0,mpi_size-1
1209          print *,'proc',i,' :   Nb Bandes  :',jj_nb_vanleer(i),
1210     &            '  : temps moyen :',
1211     &             timer_average(jj_nb_vanleer(i),timer_vanleer,i)
1212        enddo
1213     
1214        print *,'*********************************'
1215        print *,'******    TIMER DISSIP    ******'
1216        do i=0,mpi_size-1
1217          print *,'proc',i,' :   Nb Bandes  :',jj_nb_dissip(i),
1218     &            '  : temps moyen :',
1219     &             timer_average(jj_nb_dissip(i),timer_dissip,i)
1220        enddo
1221       
1222        print *,'*********************************'
1223        print *,'******    TIMER PHYSIC    ******'
1224        do i=0,mpi_size-1
1225          print *,'proc',i,' :   Nb Bandes  :',jj_nb_physic(i),
1226     &            '  : temps moyen :',
1227     &             timer_average(jj_nb_physic(i),timer_physic,i)
1228        enddo
1229       
1230      endif 
1231     
1232      print *,'Taille du Buffer MPI (REAL*8)',MaxBufferSize
1233      print *,'Taille du Buffer MPI utilise (REAL*8)',MaxBufferSize_Used
1234      print *, 'Temps total ecoule sur la parallelisation :',DiffTime()
1235      print *, 'Temps CPU ecoule sur la parallelisation :',DiffCpuTime()
1236
1237        call finalize_parallel
1238c$OMP END MASTER
1239c$OMP BARRIER
1240        RETURN
1241      ENDIF
1242     
1243      IF ( .NOT.purmats ) THEN
1244c       ........................................................
1245c       ..............  schema matsuno + leapfrog  ..............
1246c       ........................................................
1247
1248            IF(forward. OR. leapf) THEN
1249              itau= itau + 1
1250              iday= day_ini+itau/day_step
1251              time= FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
1252                IF(time.GT.1.) THEN
1253                  time = time-1.
1254                  iday = iday+1
1255                ENDIF
1256            ENDIF
1257
1258
1259            IF( itau. EQ. itaufinp1 ) then 
1260
1261c$OMP MASTER
1262              call finalize_parallel
1263c$OMP END MASTER
1264              abort_message = 'Simulation finished'
1265              call abort_gcm(modname,abort_message,0)
1266            ENDIF
1267c-----------------------------------------------------------------------
1268c   ecriture du fichier histoire moyenne:
1269c   -------------------------------------
1270
1271            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1272c$OMP BARRIER
1273               IF(itau.EQ.itaufin) THEN
1274                  iav=1
1275               ELSE
1276                  iav=0
1277               ENDIF
1278c#ifdef CPP_IOIPSL
1279c             call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1280c             call SendRequest(TestRequest)
1281cc$OMP BARRIER
1282c              call WaitRequest(TestRequest)
1283cc$OMP MASTER
1284c              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
1285c     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1286cc$OMP END MASTER
1287c               call bilan_dyn_p(2,dtvr*iperiod,dtvr*day_step*periodav,
1288c     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1289c#endif
1290            ENDIF
1291
1292c-----------------------------------------------------------------------
1293c   ecriture de la bande histoire:
1294c   ------------------------------
1295
1296c      IF( MOD(itau,iecri         ).EQ.0) THEN
1297
1298           IF( MOD(itau,iecri*day_step).EQ.0) THEN
1299c$OMP BARRIER
1300c$OMP MASTER
1301               nbetat = nbetatdem
1302        CALL geopot_p  ( ip1jmp1, teta  , pk , pks,  phis  , phi )
1303       
1304cym     unat=0.
1305       
1306        ijb=ij_begin
1307        ije=ij_end
1308       
1309        if (pole_nord) then
1310          ijb=ij_begin+iip1
1311          unat(1:iip1,:)=0.
1312        endif
1313       
1314        if (pole_sud) then
1315          ije=ij_end-iip1
1316          unat(ij_end-iip1+1:ij_end,:)=0.
1317        endif
1318           
1319        do l=1,llm
1320           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1321        enddo
1322
1323        ijb=ij_begin
1324        ije=ij_end
1325        if (pole_sud) ije=ij_end-iip1
1326       
1327        do l=1,llm
1328           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1329        enddo
1330       
1331c#ifdef CPP_IOIPSL
1332 
1333c        CALL writehist_p(histid,histvid, nqmx,itau,vcov,
1334c     s                       ucov,teta,phi,q,masse,ps,phis)
1335c#else
1336c       call Gather_Field(unat,ip1jmp1,llm,0)
1337c       call Gather_Field(vnat,ip1jm,llm,0)
1338c       call Gather_Field(teta,ip1jmp1,llm,0)
1339c       call Gather_Field(ps,ip1jmp1,1,0)
1340c       do iq=1,nqmx
1341c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1342c       enddo
1343c     
1344c       if (mpi_rank==0) then
1345c#include "write_grads_dyn.h"
1346c       endif
1347c#endif
1348c$OMP END MASTER
1349           ENDIF
1350
1351            IF(itau.EQ.itaufin) THEN
1352
1353c$OMP BARRIER
1354c$OMP MASTER
1355
1356c#ifdef CPP_IOIPSL
1357
1358       CALL dynredem1_p("restart.nc",0.0,
1359     ,                     vcov,ucov,teta,q,nqmx,masse,ps)
1360c#endif
1361
1362              CLOSE(99)
1363c$OMP END MASTER
1364            ENDIF
1365
1366c-----------------------------------------------------------------------
1367c   gestion de l'integration temporelle:
1368c   ------------------------------------
1369
1370            IF( MOD(itau,iperiod).EQ.0 )    THEN
1371                    GO TO 1
1372            ELSE IF ( MOD(itau-1,iperiod). EQ. 0 ) THEN
1373
1374                   IF( forward )  THEN
1375c      fin du pas forward et debut du pas backward
1376
1377                      forward = .FALSE.
1378                        leapf = .FALSE.
1379                           GO TO 2
1380
1381                   ELSE
1382c      fin du pas backward et debut du premier pas leapfrog
1383
1384                        leapf =  .TRUE.
1385                        dt  =  2.*dtvr
1386                        GO TO 2
1387                   END IF
1388            ELSE
1389
1390c      ......   pas leapfrog  .....
1391
1392                 leapf = .TRUE.
1393                 dt  = 2.*dtvr
1394                 GO TO 2
1395            END IF
1396
1397      ELSE
1398
1399c       ........................................................
1400c       ..............       schema  matsuno        ...............
1401c       ........................................................
1402            IF( forward )  THEN
1403
1404             itau =  itau + 1
1405             iday = day_ini+itau/day_step
1406             time = FLOAT(itau-(iday-day_ini)*day_step)/day_step+time_0
1407
1408                  IF(time.GT.1.) THEN
1409                   time = time-1.
1410                   iday = iday+1
1411                  ENDIF
1412
1413               forward =  .FALSE.
1414               IF( itau. EQ. itaufinp1 ) then 
1415c$OMP MASTER
1416                 call finalize_parallel
1417c$OMP END MASTER
1418                 abort_message = 'Simulation finished'
1419                 call abort_gcm(modname,abort_message,0)
1420               ENDIF
1421               GO TO 2
1422
1423            ELSE
1424
1425            IF(MOD(itau,iperiod).EQ.0 .OR. itau.EQ.itaufin) THEN
1426               IF(itau.EQ.itaufin) THEN
1427                  iav=1
1428               ELSE
1429                  iav=0
1430               ENDIF
1431c#ifdef CPP_IOIPSL
1432cc$OMP BARRIER
1433
1434c              call Register_Hallo(vcov,ip1jm,llm,1,0,0,1,TestRequest)
1435c              call SendRequest(TestRequest)
1436cc$OMP BARRIER
1437c              call WaitRequest(TestRequest)
1438
1439cc$OMP MASTER
1440c              CALL writedynav_p(histaveid, nqmx, itau,vcov ,
1441c     ,                          ucov,teta,pk,phi,q,masse,ps,phis)
1442c               call bilan_dyn_p (2,dtvr*iperiod,dtvr*day_step*periodav,
1443c     ,           ps,masse,pk,pbaru,pbarv,teta,phi,ucov,vcov,q)
1444cc$OMP END MASTER
1445c#endif
1446            ENDIF
1447
1448c               IF(MOD(itau,iecri         ).EQ.0) THEN
1449              IF(MOD(itau,iecri*day_step).EQ.0) THEN
1450c$OMP BARRIER
1451c$OMP MASTER
1452                  nbetat = nbetatdem
1453       CALL geopot_p( ip1jmp1, teta  , pk , pks,  phis  , phi   )
1454
1455cym        unat=0.
1456        ijb=ij_begin
1457        ije=ij_end
1458       
1459        if (pole_nord) then
1460          ijb=ij_begin+iip1
1461          unat(1:iip1,:)=0.
1462        endif
1463       
1464        if (pole_sud) then
1465          ije=ij_end-iip1
1466          unat(ij_end-iip1+1:ij_end,:)=0.
1467        endif
1468           
1469        do l=1,llm
1470           unat(ijb:ije,l)=ucov(ijb:ije,l)/cu(ijb:ije)
1471        enddo
1472
1473        ijb=ij_begin
1474        ije=ij_end
1475        if (pole_sud) ije=ij_end-iip1
1476       
1477        do l=1,llm
1478           vnat(ijb:ije,l)=vcov(ijb:ije,l)/cv(ijb:ije)
1479        enddo
1480
1481c#ifdef CPP_IOIPSL
1482
1483c       CALL writehist_p( histid, histvid, nqmx, itau,vcov ,
1484c     ,                           ucov,teta,phi,q,masse,ps,phis)
1485c#else
1486c      call Gather_Field(unat,ip1jmp1,llm,0)
1487c      call Gather_Field(vnat,ip1jm,llm,0)
1488c      call Gather_Field(teta,ip1jmp1,llm,0)
1489c      call Gather_Field(ps,ip1jmp1,1,0)
1490c      do iq=1,nqmx
1491c        call Gather_Field(q(1,1,iq),ip1jmp1,llm,0)
1492c      enddo
1493c     
1494c      if (mpi_rank==0) then
1495c#include "write_grads_dyn.h"
1496c      endif
1497c#endif
1498
1499c$OMP END MASTER
1500               ENDIF
1501
1502                 IF(itau.EQ.itaufin) THEN
1503c$OMP MASTER
1504                   CALL dynredem1_p("restart.nc",0.0,
1505     .                               vcov,ucov,teta,q,nqmx,masse,ps)
1506c$OMP END MASTER
1507                 ENDIF
1508                 forward = .TRUE.
1509                 GO TO  1
1510
1511            ENDIF
1512
1513      END IF
1514c$OMP MASTER
1515        call finalize_parallel
1516c$OMP END MASTER
1517        RETURN
1518      END
Note: See TracBrowser for help on using the repository browser.