source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/leapfrog_p.F @ 1110

Last change on this file since 1110 was 1000, checked in by Laurent Fairhead, 16 years ago
  • Modifs sur le parallelisme: masquage dans la physique
  • Inclusion strato
  • mise en coherence etat0
  • le mode offline fonctionne maintenant en parallele,
  • les fichiers de la dynamiques sont correctement sortis et peuvent etre reconstruit avec rebuild
  • la version parallele de la dynamique peut s'executer sans MPI (sur 1 proc)
  • L'OPENMP fonctionne maintenant sans la parallelisation MPI.

YM
LF

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