source: LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_loc.F @ 1646

Last change on this file since 1646 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT


Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

File size: 23.6 KB
Line 
1!
2! $Id: bilan_dyn_p.F 1299 2010-01-20 14:27:21Z fairhead $
3!
4      SUBROUTINE bilan_dyn_loc (ntrac,dt_app,dt_cum,
5     s  ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
6
7c   AFAIRE
8c   Prevoir en champ nq+1 le diagnostique de l'energie
9c   en faisant Qzon=Cv T + L * ...
10c             vQ..A=Cp T + L * ...
11
12#ifdef CPP_IOIPSL
13      USE IOIPSL
14#endif
15      USE parallel
16      USE mod_hallo
17      use misc_mod
18      use write_field
19      IMPLICIT NONE
20
21#include "dimensions.h"
22#include "paramet.h"
23#include "comconst.h"
24#include "comvert.h"
25#include "comgeom2.h"
26#include "temps.h"
27#include "iniprint.h"
28
29c====================================================================
30c
31c   Sous-programme consacre à des diagnostics dynamiques de base
32c
33c
34c   De facon generale, les moyennes des scalaires Q sont ponderees par
35c   la masse.
36c
37c   Les flux de masse sont eux simplement moyennes.
38c
39c====================================================================
40
41c   Arguments :
42c   ===========
43
44      integer ntrac
45      real dt_app,dt_cum
46      real ps(iip1,jjb_u:jje_u)
47      real masse(iip1,jjb_u:jje_u,llm),pk(iip1,jjb_u:jje_u,llm)
48      real flux_u(iip1,jjb_u:jje_u,llm)
49      real flux_v(iip1,jjb_v:jje_v,llm)
50      real teta(iip1,jjb_u:jje_u,llm)
51      real phi(iip1,jjb_u:jje_u,llm)
52      real ucov(iip1,jjb_u:jje_u,llm)
53      real vcov(iip1,jjb_v:jje_v,llm)
54      real trac(iip1,jjb_u:jje_u,llm,ntrac)
55
56c   Local :
57c   =======
58
59      integer,SAVE :: icum,ncum
60!$OMP THREADPRIVATE(icum,ncum)
61      LOGICAL,SAVE :: first=.TRUE.
62!$OMP THREADPRIVATE(first)     
63     
64      real zz,zqy
65      REAl,SAVE,ALLOCATABLE :: zfactv(:,:)
66
67      INTEGER,PARAMETER :: nQ=7
68
69
70cym      character*6 nom(nQ)
71cym      character*6 unites(nQ)
72      character(len=6),save :: nom(nQ)
73      character(len=6),save :: unites(nQ)
74
75      character(len=10) file
76      integer ifile
77      parameter (ifile=4)
78
79      integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
80      INTEGER,PARAMETER :: iovap=6,iun=7
81      integer,PARAMETER :: i_sortie=1
82
83      real,SAVE :: time=0.
84      integer,SAVE :: itau=0.
85!$OMP THREADPRIVATE(time,itau)
86
87      real ww
88
89c   variables dynamiques intermédiaires
90      REAL,SAVE,ALLOCATABLE :: vcont(:,:,:),ucont(:,:,:)
91      REAL,SAVE,ALLOCATABLE :: ang(:,:,:),unat(:,:,:)
92      REAL,SAVE,ALLOCATABLE :: massebx(:,:,:),masseby(:,:,:)
93      REAL,SAVE,ALLOCATABLE :: vorpot(:,:,:)
94      REAL,SAVE,ALLOCATABLE :: w(:,:,:),ecin(:,:,:),convm(:,:,:)
95      REAL,SAVE,ALLOCATABLE :: bern(:,:,:)
96
97c   champ contenant les scalaires advectés.
98      real,SAVE,ALLOCATABLE :: Q(:,:,:,:)
99   
100c   champs cumulés
101      real,SAVE,ALLOCATABLE ::  ps_cum(:,:)
102      real,SAVE,ALLOCATABLE ::  masse_cum(:,:,:)
103      real,SAVE,ALLOCATABLE ::  flux_u_cum(:,:,:)
104      real,SAVE,ALLOCATABLE ::  flux_v_cum(:,:,:)
105      real,SAVE,ALLOCATABLE ::  Q_cum(:,:,:,:)
106      real,SAVE,ALLOCATABLE ::  flux_uQ_cum(:,:,:,:)
107      real,SAVE,ALLOCATABLE ::  flux_vQ_cum(:,:,:,:)
108      real,SAVE,ALLOCATABLE ::  flux_wQ_cum(:,:,:,:)
109      real,SAVE,ALLOCATABLE ::  dQ(:,:,:,:)
110
111 
112c   champs de tansport en moyenne zonale
113      integer ntr,itr
114      parameter (ntr=5)
115
116cym      character*10 znom(ntr,nQ)
117cym      character*20 znoml(ntr,nQ)
118cym      character*10 zunites(ntr,nQ)
119      character*10,save :: znom(ntr,nQ)
120      character*20,save :: znoml(ntr,nQ)
121      character*10,save :: zunites(ntr,nQ)
122
123      INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5
124
125      character*3 ctrs(ntr)
126      data ctrs/'  ','TOT','MMC','TRS','STN'/
127
128      real,SAVE,ALLOCATABLE ::  zvQ(:,:,:,:),zvQtmp(:,:)
129      real,SAVE,ALLOCATABLE ::  zavQ(:,:,:),psiQ(:,:,:)
130      real,SAVE,ALLOCATABLE ::  zmasse(:,:),zamasse(:)
131
132      real,SAVE,ALLOCATABLE ::  zv(:,:),psi(:,:)
133
134      integer i,j,l,iQ
135
136
137c   Initialisation du fichier contenant les moyennes zonales.
138c   ---------------------------------------------------------
139
140      character*10 infile
141
142      integer fileid
143      integer thoriid, zvertiid
144      save fileid
145
146      INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)
147
148C   Variables locales
149C
150      integer tau0
151      real zjulian
152      character*3 str
153      character*10 ctrac
154      integer ii,jj
155      integer zan, dayref
156C
157      real,SAVE,ALLOCATABLE :: rlong(:),rlatg(:)
158      integer :: jjb,jje,jjn,ijb,ije
159      type(Request) :: Req
160
161! definition du domaine d'ecriture pour le rebuild
162
163      INTEGER,DIMENSION(1) :: ddid
164      INTEGER,DIMENSION(1) :: dsg
165      INTEGER,DIMENSION(1) :: dsl
166      INTEGER,DIMENSION(1) :: dpf
167      INTEGER,DIMENSION(1) :: dpl
168      INTEGER,DIMENSION(1) :: dhs
169      INTEGER,DIMENSION(1) :: dhe
170     
171      INTEGER :: bilan_dyn_domain_id
172
173
174c=====================================================================
175c   Initialisation
176c=====================================================================
177      if (adjust) return
178     
179      time=time+dt_app
180      itau=itau+1
181
182      if (first) then
183!$OMP BARRIER
184!$OMP MASTER
185      ALLOCATE(zfactv(jjb_v:jje_v,llm))
186      ALLOCATE(vcont(iip1,jjb_v:jje_v,llm))
187      ALLOCATE(ucont(iip1,jjb_u:jje_u,llm))
188      ALLOCATE(ang(iip1,jjb_u:jje_u,llm))
189      ALLOCATE(unat(iip1,jjb_u:jje_u,llm))
190      ALLOCATE(massebx(iip1,jjb_u:jje_u,llm))
191      ALLOCATE(masseby(iip1,jjb_v:jje_v,llm))
192      ALLOCATE(vorpot(iip1,jjb_v:jje_v,llm))
193      ALLOCATE(w(iip1,jjb_u:jje_u,llm))
194      ALLOCATE(ecin(iip1,jjb_u:jje_u,llm))
195      ALLOCATE(convm(iip1,jjb_u:jje_u,llm))
196      ALLOCATE(bern(iip1,jjb_u:jje_u,llm))     
197      ALLOCATE(Q(iip1,jjb_u:jje_u,llm,nQ))     
198      ALLOCATE(ps_cum(iip1,jjb_u:jje_u))
199      ALLOCATE(masse_cum(iip1,jjb_u:jje_u,llm))
200      ALLOCATE(flux_u_cum(iip1,jjb_u:jje_u,llm))
201      ALLOCATE(flux_v_cum(iip1,jjb_v:jje_v,llm))
202      ALLOCATE(Q_cum(iip1,jjb_u:jje_u,llm,nQ))
203      ALLOCATE(flux_uQ_cum(iip1,jjb_u:jje_u,llm,nQ))
204      ALLOCATE(flux_vQ_cum(iip1,jjb_v:jje_v,llm,nQ))
205      ALLOCATE(flux_wQ_cum(iip1,jjb_u:jje_u,llm,nQ))
206      ALLOCATE(dQ(iip1,jjb_u:jje_u,llm,nQ))
207      ALLOCATE(zvQ(jjb_v:jje_v,llm,ntr,nQ))
208      ALLOCATE(zvQtmp(jjb_v:jje_v,llm))
209      ALLOCATE(zavQ(jjb_v:jje_v,ntr,nQ))
210      ALLOCATE(psiQ(jjb_v:jje_v,llm+1,nQ))
211      ALLOCATE(zmasse(jjb_v:jje_v,llm))
212      ALLOCATE(zamasse(jjb_v:jje_v))
213      ALLOCATE(zv(jjb_v:jje_v,llm))
214      ALLOCATE(psi(jjb_v:jje_v,llm+1))
215      ALLOCATE(ndex3d(jjb_v:jje_v*llm))
216      ndex3d=0
217      ALLOCATE(rlong(jjb_v:jje_v))
218      ALLOCATE(rlatg(jjb_v:jje_v))
219     
220!$OMP END MASTER
221!$OMP BARRIER
222        icum=0
223c       initialisation des fichiers
224        first=.false.
225c   ncum est la frequence de stokage en pas de temps
226        ncum=dt_cum/dt_app
227        if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
228           WRITE(lunout,*)
229     .            'Pb : le pas de cumule doit etre multiple du pas'
230           WRITE(lunout,*)'dt_app=',dt_app
231           WRITE(lunout,*)'dt_cum=',dt_cum
232           stop
233        endif
234
235!$OMP MASTER
236        nom(itemp)='T'
237        nom(igeop)='gz'
238        nom(iecin)='K'
239        nom(iang)='ang'
240        nom(iu)='u'
241        nom(iovap)='ovap'
242        nom(iun)='un'
243
244        unites(itemp)='K'
245        unites(igeop)='m2/s2'
246        unites(iecin)='m2/s2'
247        unites(iang)='ang'
248        unites(iu)='m/s'
249        unites(iovap)='kg/kg'
250        unites(iun)='un'
251
252
253c   Initialisation du fichier contenant les moyennes zonales.
254c   ---------------------------------------------------------
255
256      infile='dynzon'
257
258      zan = annee_ref
259      dayref = day_ref
260      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
261      tau0 = itau_dyn
262     
263      rlong=0.
264      rlatg=rlatv*180./pi
265
266      jjb=jj_begin
267      jje=jj_end
268      jjn=jj_nb
269      IF (pole_sud) THEN
270        jjn=jj_nb-1
271        jje=jj_end-1
272      ENDIF
273
274      ddid=(/ 2 /)
275      dsg=(/ jjm /)
276      dsl=(/ jjn /)
277      dpf=(/ jjb /)
278      dpl=(/ jje /)
279      dhs=(/ 0 /)
280      dhe=(/ 0 /)
281
282      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
283     .                 'box',bilan_dyn_domain_id)
284       
285      call histbeg(trim(infile),
286     .             1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
287     .             1, 1, 1, jjn,
288     .             tau0, zjulian, dt_cum, thoriid, fileid,
289     .             bilan_dyn_domain_id)
290
291C
292C  Appel a histvert pour la grille verticale
293C
294      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
295     .              llm, presnivs, zvertiid)
296C
297C  Appels a histdef pour la definition des variables a sauvegarder
298      do iQ=1,nQ
299         do itr=1,ntr
300            if(itr.eq.1) then
301               znom(itr,iQ)=nom(iQ)
302               znoml(itr,iQ)=nom(iQ)
303               zunites(itr,iQ)=unites(iQ)
304            else
305               znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
306               znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
307               zunites(itr,iQ)='m/s * '//unites(iQ)
308            endif
309         enddo
310      enddo
311
312c   Declarations des champs avec dimension verticale
313c      print*,'1HISTDEF'
314      do iQ=1,nQ
315         do itr=1,ntr
316      IF (prt_level > 5)
317     . WRITE(lunout,*)'var ',itr,iQ
318     .      ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
319            call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
320     .        zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
321     .        32,'ave(X)',dt_cum,dt_cum)
322         enddo
323c   Declarations pour les fonctions de courant
324c      print*,'2HISTDEF'
325          call histdef(fileid,'psi'//nom(iQ)
326     .      ,'stream fn. '//znoml(itot,iQ),
327     .      zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
328     .      32,'ave(X)',dt_cum,dt_cum)
329      enddo
330
331
332c   Declarations pour les champs de transport d'air
333c      print*,'3HISTDEF'
334      call histdef(fileid, 'masse', 'masse',
335     .             'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
336     .             32, 'ave(X)', dt_cum, dt_cum)
337      call histdef(fileid, 'v', 'v',
338     .             'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
339     .             32, 'ave(X)', dt_cum, dt_cum)
340c   Declarations pour les fonctions de courant
341c      print*,'4HISTDEF'
342          call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
343     .      1,jjn,thoriid,llm,1,llm,zvertiid,
344     .      32,'ave(X)',dt_cum,dt_cum)
345
346
347c   Declaration des champs 1D de transport en latitude
348c      print*,'5HISTDEF'
349      do iQ=1,nQ
350         do itr=2,ntr
351            call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
352     .        zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
353     .        32,'ave(X)',dt_cum,dt_cum)
354         enddo
355      enddo
356
357
358c      print*,'8HISTDEF'
359               CALL histend(fileid)
360
361!$OMP END MASTER
362      endif
363
364
365c=====================================================================
366c   Calcul des champs dynamiques
367c   ----------------------------
368
369      jjb=jj_begin
370      jje=jj_end
371   
372c   énergie cinétique
373!      ucont(:,jjb:jje,:)=0
374
375      call Register_Hallo_u(ucov,llm,1,1,1,1,Req)
376      call Register_Hallo_v(vcov,llm,1,1,1,1,Req)
377      call SendRequest(Req)
378c$OMP BARRIER
379      call WaitRequest(Req)
380
381      CALL covcont_loc(llm,ucov,vcov,ucont,vcont)
382      CALL enercin_loc(vcov,ucov,vcont,ucont,ecin)
383
384c   moment cinétique
385!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
386      do l=1,llm
387         ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
388         unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
389      enddo
390!$OMP END DO NOWAIT
391
392!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
393      DO l=1,llm
394        Q(:,jjb:jje,l,itemp)=teta(:,jjb:jje,l)*pk(:,jjb:jje,l)/cpp
395        Q(:,jjb:jje,l,igeop)=phi(:,jjb:jje,l)
396        Q(:,jjb:jje,l,iecin)=ecin(:,jjb:jje,l)
397        Q(:,jjb:jje,l,iang)=ang(:,jjb:jje,l)
398        Q(:,jjb:jje,l,iu)=unat(:,jjb:jje,l)
399        Q(:,jjb:jje,l,iovap)=trac(:,jjb:jje,l,1)
400        Q(:,jjb:jje,l,iun)=1.
401      ENDDO
402!$OMP END DO NOWAIT
403
404c=====================================================================
405c   Cumul
406c=====================================================================
407c
408      if(icum.EQ.0) then
409         jjb=jj_begin
410         jje=jj_end
411
412!$OMP MASTER
413         ps_cum(:,jjb:jje)=0.
414!$OMP END MASTER
415
416
417!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
418        DO l=1,llm
419          masse_cum(:,jjb:jje,l)=0.
420          flux_u_cum(:,jjb:jje,l)=0.
421          Q_cum(:,jjb:jje,:,l)=0.
422          flux_uQ_cum(:,jjb:jje,l,:)=0.
423          flux_v_cum(:,jjb:jje,l)=0.
424          if (pole_sud) jje=jj_end-1
425          flux_v_cum(:,jjb:jje,l)=0.
426          flux_vQ_cum(:,jjb:jje,l,:)=0.
427        ENDDO
428!$OMP END DO NOWAIT
429      endif
430
431      IF (prt_level > 5)
432     . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
433      icum=icum+1
434
435c   accumulation des flux de masse horizontaux
436      jjb=jj_begin
437      jje=jj_end
438
439!$OMP MASTER
440      ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
441!$OMP END MASTER
442
443
444!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
445      DO l=1,llm
446        masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
447        flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
448     .                         +flux_u(:,jjb:jje,:)
449      ENDDO
450!$OMP END DO NOWAIT
451     
452      if (pole_sud) jje=jj_end-1
453
454!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
455      DO l=1,llm
456       flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)
457     .                          +flux_v(:,jjb:jje,:)
458      ENDDO
459!$OMP END DO NOWAIT
460     
461      jjb=jj_begin
462      jje=jj_end
463
464      do iQ=1,nQ
465!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
466        DO l=1,llm
467          Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
468     .                       +Q(:,jjb:jje,:,iQ)*masse(:,jjb:jje,:)
469        ENDDO
470!$OMP END DO NOWAIT
471      enddo
472
473c=====================================================================
474c  FLUX ET TENDANCES
475c=====================================================================
476
477c   Flux longitudinal
478c   -----------------
479      do iQ=1,nQ
480!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
481         do l=1,llm
482            do j=jjb,jje
483               do i=1,iim
484                  flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
485     s            +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
486               enddo
487               flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
488            enddo
489         enddo
490!$OMP END DO NOWAIT
491      enddo
492
493c    flux méridien
494c    -------------
495      do iQ=1,nQ
496        call Register_Hallo_u(Q(1,1,1,iQ),llm,0,1,1,0,Req)
497      enddo
498      call SendRequest(Req)
499!$OMP BARRIER     
500      call WaitRequest(Req)
501     
502      jjb=jj_begin
503      jje=jj_end
504      if (pole_sud) jje=jj_end-1
505     
506      do iQ=1,nQ
507!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
508         do l=1,llm
509            do j=jjb,jje
510               do i=1,iip1
511                  flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
512     s            +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
513               enddo
514            enddo
515         enddo
516!$OMP END DO NOWAIT
517      enddo
518
519
520c    tendances
521c    ---------
522
523c   convergence horizontale
524      call Register_Hallo_u(flux_uQ_cum,llm,2,2,2,2,Req)
525      call Register_Hallo_v(flux_vQ_cum,llm,2,2,2,2,Req)
526      call SendRequest(Req)
527!$OMP BARRIER     
528      call WaitRequest(Req)
529
530      call  convflu_loc(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
531
532c   calcul de la vitesse verticale
533      call Register_Hallo_u(flux_u_cum,llm,2,2,2,2,Req)
534      call Register_Hallo_v(flux_v_cum,llm,2,2,2,2,Req)
535      call SendRequest(Req)
536!$OMP BARRIER     
537      call WaitRequest(Req)
538
539      call convmas_loc(flux_u_cum,flux_v_cum,convm)
540      CALL vitvert_loc(convm,w)
541!$OMP BARRIER
542
543      jjb=jj_begin
544      jje=jj_end
545
546!      do iQ=1,nQ
547!         do l=1,llm-1
548!            do j=jjb,jje
549!               do i=1,iip1
550!                  ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
551!                  dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
552!                  dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
553!               enddo
554!            enddo
555!          enddo
556!       enddo
557       
558      do iQ=1,nQ
559!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
560         do l=1,llm
561            IF (l<llm) THEN
562              do j=jjb,jje
563                 do i=1,iip1
564                    ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
565                    dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
566                    dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
567                 enddo
568              enddo
569            ENDIF
570            IF (l>2) THEN
571              do j=jjb,jje
572                do i=1,iip1
573                  ww=-0.5*w(i,j,l)*(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
574                  dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
575                enddo
576              enddo
577            ENDIF
578         enddo
579!$OMP ENDDO NOWAIT
580      enddo
581      IF (prt_level > 5)
582     . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
583c=====================================================================
584c   PAS DE TEMPS D'ECRITURE
585c=====================================================================
586      if (icum.eq.ncum) then
587c=====================================================================
588
589      IF (prt_level > 5)
590     . WRITE(lunout,*)'Pas d ecriture'
591
592      jjb=jj_begin
593      jje=jj_end
594
595c   Normalisation
596      do iQ=1,nQ
597!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
598        do l=1,llm
599          Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
600     .                          /masse_cum(:,jjb:jje,l)
601        enddo
602!$OMP ENDDO NOWAIT
603      enddo   
604
605      zz=1./REAL(ncum)
606
607!$OMP MASTER
608        ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
609!$OMP END MASTER
610
611!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
612      DO l=1,llm
613        masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
614        flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
615        flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
616        dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
617      ENDDO
618!$OMP ENDDO NOWAIT
619         
620     
621      IF (pole_sud) jje=jj_end-1
622!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
623      DO l=1,llm
624        flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)*zz
625        flux_vQ_cum(:,jjb:jje,l,:)=flux_vQ_cum(:,jjb:jje,l,:)*zz
626      ENDDO
627!$OMP ENDDO NOWAIT
628         
629      jjb=jj_begin
630      jje=jj_end
631
632
633c   A retravailler eventuellement
634c   division de dQ par la masse pour revenir aux bonnes grandeurs
635      do iQ=1,nQ
636!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
637        DO l=1,llm
638           dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
639        ENDDO
640!$OMP ENDDO NOWAIT
641      enddo
642 
643c=====================================================================
644c   Transport méridien
645c=====================================================================
646
647c   cumul zonal des masses des mailles
648c   ----------------------------------
649      jjb=jj_begin
650      jje=jj_end
651      if (pole_sud) jje=jj_end-1
652
653!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
654        DO l=1,llm
655          zv(jjb:jje,l)=0.
656          zmasse(jjb:jje,l)=0.
657        ENDDO
658!$OMP ENDDO NOWAIT
659
660      call Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
661      do iQ=1,nQ
662        call Register_Hallo_u(Q_cum(1,1,1,iQ),llm,0,1,1,0,Req)
663      enddo
664
665      call SendRequest(Req)
666!$OMP BARRIER
667      call WaitRequest(Req)
668
669      call massbar_loc(masse_cum,massebx,masseby)
670     
671      jjb=jj_begin
672      jje=jj_end
673      if (pole_sud) jje=jj_end-1
674     
675!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
676      do l=1,llm
677         do j=jjb,jje
678            do i=1,iim
679               zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
680               zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
681            enddo
682            zfactv(j,l)=cv(1,j)/zmasse(j,l)
683         enddo
684      enddo
685!$OMP ENDDO NOWAIT
686
687c     print*,'3OK'
688c   --------------------------------------------------------------
689c   calcul de la moyenne zonale du transport :
690c   ------------------------------------------
691c
692c                                     --
693c TOT : la circulation totale       [ vq ]
694c
695c                                      -     -
696c MMC : mean meridional circulation [ v ] [ q ]
697c
698c                                     ----      --       - -
699c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
700c
701c                                     - * - *       - -       -     -
702c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
703c
704c                                              - -
705c    on utilise aussi l'intermediaire TMP :  [ v q ]
706c
707c    la variable zfactv transforme un transport meridien cumule
708c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
709c
710c   --------------------------------------------------------------
711
712
713c   ----------------------------------------
714c   Transport dans le plan latitude-altitude
715c   ----------------------------------------
716
717      jjb=jj_begin
718      jje=jj_end
719      if (pole_sud) jje=jj_end-1
720     
721      zvQ=0.
722      psiQ=0.
723      do iQ=1,nQ
724!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
725         do l=1,llm
726            zvQtmp(:,l)=0.
727            do j=jjb,jje
728c              print*,'j,l,iQ=',j,l,iQ
729c   Calcul des moyennes zonales du transort total et de zvQtmp
730               do i=1,iim
731                  zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
732     s                            +flux_vQ_cum(i,j,l,iQ)
733                  zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
734     s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
735                  zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
736     s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
737                  zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
738               enddo
739c              print*,'aOK'
740c   Decomposition
741               zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
742               zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
743               zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
744               zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
745               zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
746               zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
747            enddo
748         enddo
749!$OMP ENDDO NOWAIT
750c   fonction de courant meridienne pour la quantite Q
751!$OMP BARRIER
752!$OMP MASTER
753         do l=llm,1,-1
754            do j=jjb,jje
755               psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
756            enddo
757         enddo
758!$OMP END MASTER
759!$OMP BARRIER
760      enddo
761
762c   fonction de courant pour la circulation meridienne moyenne
763!$OMP BARRIER
764!$OMP MASTER
765      psi(jjb:jje,:)=0.
766      do l=llm,1,-1
767         do j=jjb,jje
768            psi(j,l)=psi(j,l+1)+zv(j,l)
769            zv(j,l)=zv(j,l)*zfactv(j,l)
770         enddo
771      enddo
772!$OMP END MASTER
773!$OMP BARRIER
774
775c     print*,'4OK'
776c   sorties proprement dites
777!$OMP MASTER     
778      if (i_sortie.eq.1) then
779      jjb=jj_begin
780      jje=jj_end
781      jjn=jj_nb
782      if (pole_sud) jje=jj_end-1
783      if (pole_sud) jjn=jj_nb-1
784      do iQ=1,nQ
785         do itr=1,ntr
786            call histwrite(fileid,znom(itr,iQ),itau,
787     s                     zvQ(jjb:jje,:,itr,iQ)
788     s                     ,jjn*llm,ndex3d)
789         enddo
790         call histwrite(fileid,'psi'//nom(iQ),
791     s                  itau,psiQ(jjb:jje,1:llm,iQ)
792     s                  ,jjn*llm,ndex3d)
793      enddo
794
795      call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
796     s   ,jjn*llm,ndex3d)
797      call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
798     s   ,jjn*llm,ndex3d)
799      psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
800      call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
801     s               jjn*llm,ndex3d)
802
803      endif
804
805 
806c   -----------------
807c   Moyenne verticale
808c   -----------------
809
810      zamasse(jjb:jje)=0.
811      do l=1,llm
812         zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
813      enddo
814     
815      zavQ(jjb:jje,:,:)=0.
816      do iQ=1,nQ
817         do itr=2,ntr
818            do l=1,llm
819               zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
820     s                             +zvQ(jjb:jje,l,itr,iQ)
821     s                             *zmasse(jjb:jje,l)
822            enddo
823            zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
824            call histwrite(fileid,'a'//znom(itr,iQ),itau,
825     s                     zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
826         enddo
827      enddo
828!$OMP END MASTER
829c     on doit pouvoir tracer systematiquement la fonction de courant.
830
831c=====================================================================
832c/////////////////////////////////////////////////////////////////////
833      icum=0                  !///////////////////////////////////////
834      endif ! icum.eq.ncum    !///////////////////////////////////////
835c/////////////////////////////////////////////////////////////////////
836c=====================================================================
837
838      return
839      end
Note: See TracBrowser for help on using the repository browser.