My Project
 All Classes Files Functions Variables Macros
bilan_dyn_loc.F
Go to the documentation of this file.
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 
7 c AFAIRE
8 c Prevoir en champ nq+1 le diagnostique de l'energie
9 c en faisant Qzon=Cv T + L * ...
10 c 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 
29 c====================================================================
30 c
31 c Sous-programme consacre à des diagnostics dynamiques de base
32 c
33 c
34 c De facon generale, les moyennes des scalaires Q sont ponderees par
35 c la masse.
36 c
37 c Les flux de masse sont eux simplement moyennes.
38 c
39 c====================================================================
40 
41 c Arguments :
42 c ===========
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 
56 c Local :
57 c =======
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 
70 cym character*6 nom(nQ)
71 cym 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 
89 c 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 
97 c champ contenant les scalaires advectés.
98  real,SAVE,ALLOCATABLE :: q(:,:,:,:)
99 
100 c 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 
112 c champs de tansport en moyenne zonale
113  integer ntr,itr
114  parameter(ntr=5)
115 
116 cym character*10 znom(ntr,nQ)
117 cym character*20 znoml(ntr,nQ)
118 cym 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 
137 c Initialisation du fichier contenant les moyennes zonales.
138 c ---------------------------------------------------------
139 
140  character*10 infile
141 
142  integer fileid
143  integer thoriid, zvertiid
144  save fileid
145 
146  INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)
147 
148 C Variables locales
149 C
150  integer tau0
151  real zjulian
152  character*3 str
153  character*10 ctrac
154  integer ii,jj
155  integer zan, dayref
156 C
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 
174 c=====================================================================
175 c Initialisation
176 c=====================================================================
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
223 c initialisation des fichiers
224  first=.false.
225 c 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 
253 c Initialisation du fichier contenant les moyennes zonales.
254 c ---------------------------------------------------------
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 
291 C
292 C Appel a histvert pour la grille verticale
293 C
294  call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
295  . llm, presnivs, zvertiid)
296 C
297 C 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 
312 c Declarations des champs avec dimension verticale
313 c 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
323 c Declarations pour les fonctions de courant
324 c 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 
332 c Declarations pour les champs de transport d'air
333 c 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)
340 c Declarations pour les fonctions de courant
341 c 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 
347 c Declaration des champs 1D de transport en latitude
348 c 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 
358 c print*,'8HISTDEF'
359  CALL histend(fileid)
360 
361 !$OMP END MASTER
362  endif
363 
364 
365 c=====================================================================
366 c Calcul des champs dynamiques
367 c ----------------------------
368 
369  jjb=jj_begin
370  jje=jj_end
371 
372 c é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)
378 c$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 
384 c 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 
404 c=====================================================================
405 c Cumul
406 c=====================================================================
407 c
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  if (pole_sud) jje=jj_end-1
424  flux_v_cum(:,jjb:jje,l)=0.
425  flux_vq_cum(:,jjb:jje,l,:)=0.
426  ENDDO
427 !$OMP END DO NOWAIT
428  endif
429 
430  IF (prt_level > 5)
431  . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
432  icum=icum+1
433 
434 c accumulation des flux de masse horizontaux
435  jjb=jj_begin
436  jje=jj_end
437 
438 !$OMP MASTER
439  ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
440 !$OMP END MASTER
441 
442 
443 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
444  DO l=1,llm
445  masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
446  flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
447  . +flux_u(:,jjb:jje,:)
448  ENDDO
449 !$OMP END DO NOWAIT
450 
451  if (pole_sud) jje=jj_end-1
452 
453 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
454  DO l=1,llm
455  flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)
456  . +flux_v(:,jjb:jje,:)
457  ENDDO
458 !$OMP END DO NOWAIT
459 
460  jjb=jj_begin
461  jje=jj_end
462 
463  do iq=1,nq
464 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
465  DO l=1,llm
466  q_cum(:,jjb:jje,:,iq)=q_cum(:,jjb:jje,:,iq)
467  . +q(:,jjb:jje,:,iq)*masse(:,jjb:jje,:)
468  ENDDO
469 !$OMP END DO NOWAIT
470  enddo
471 
472 c=====================================================================
473 c FLUX ET TENDANCES
474 c=====================================================================
475 
476 c Flux longitudinal
477 c -----------------
478  do iq=1,nq
479 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
480  do l=1,llm
481  do j=jjb,jje
482  do i=1,iim
483  flux_uq_cum(i,j,l,iq)=flux_uq_cum(i,j,l,iq)
484  s +flux_u(i,j,l)*0.5*(q(i,j,l,iq)+q(i+1,j,l,iq))
485  enddo
486  flux_uq_cum(iip1,j,l,iq)=flux_uq_cum(1,j,l,iq)
487  enddo
488  enddo
489 !$OMP END DO NOWAIT
490  enddo
491 
492 c flux méridien
493 c -------------
494  do iq=1,nq
495  call register_hallo_u(q(1,1,1,iq),llm,0,1,1,0,req)
496  enddo
497  call sendrequest(req)
498 !$OMP BARRIER
499  call waitrequest(req)
500 
501  jjb=jj_begin
502  jje=jj_end
503  if (pole_sud) jje=jj_end-1
504 
505  do iq=1,nq
506 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
507  do l=1,llm
508  do j=jjb,jje
509  do i=1,iip1
510  flux_vq_cum(i,j,l,iq)=flux_vq_cum(i,j,l,iq)
511  s +flux_v(i,j,l)*0.5*(q(i,j,l,iq)+q(i,j+1,l,iq))
512  enddo
513  enddo
514  enddo
515 !$OMP END DO NOWAIT
516  enddo
517 
518 
519 c tendances
520 c ---------
521 
522 c convergence horizontale
523  call register_hallo_u(flux_uq_cum,llm,2,2,2,2,req)
524  call register_hallo_v(flux_vq_cum,llm,2,2,2,2,req)
525  call sendrequest(req)
526 !$OMP BARRIER
527  call waitrequest(req)
528 
529  call convflu_loc(flux_uq_cum,flux_vq_cum,llm*nq,dq)
530 
531 c calcul de la vitesse verticale
532  call register_hallo_u(flux_u_cum,llm,2,2,2,2,req)
533  call register_hallo_v(flux_v_cum,llm,2,2,2,2,req)
534  call sendrequest(req)
535 !$OMP BARRIER
536  call waitrequest(req)
537 
538  call convmas_loc(flux_u_cum,flux_v_cum,convm)
539  CALL vitvert_loc(convm,w)
540 !$OMP BARRIER
541 
542  jjb=jj_begin
543  jje=jj_end
544 
545 ! do iQ=1,nQ
546 ! do l=1,llm-1
547 ! do j=jjb,jje
548 ! do i=1,iip1
549 ! ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
550 ! dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
551 ! dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
552 ! enddo
553 ! enddo
554 ! enddo
555 ! enddo
556 
557  do iq=1,nq
558 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
559  do l=1,llm
560  IF (l<llm) THEN
561  do j=jjb,jje
562  do i=1,iip1
563  ww=-0.5*w(i,j,l+1)*(q(i,j,l,iq)+q(i,j,l+1,iq))
564  dq(i,j,l ,iq)=dq(i,j,l ,iq)-ww
565  dq(i,j,l+1,iq)=dq(i,j,l+1,iq)+ww
566  enddo
567  enddo
568  ENDIF
569  IF (l>2) THEN
570  do j=jjb,jje
571  do i=1,iip1
572  ww=-0.5*w(i,j,l)*(q(i,j,l-1,iq)+q(i,j,l,iq))
573  dq(i,j,l,iq)=dq(i,j,l,iq)+ww
574  enddo
575  enddo
576  ENDIF
577  enddo
578 !$OMP ENDDO NOWAIT
579  enddo
580  IF (prt_level > 5)
581  . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
582 c=====================================================================
583 c PAS DE TEMPS D'ECRITURE
584 c=====================================================================
585  if (icum.eq.ncum) then
586 c=====================================================================
587 
588  IF (prt_level > 5)
589  . WRITE(lunout,*)'Pas d ecriture'
590 
591  jjb=jj_begin
592  jje=jj_end
593 
594 c Normalisation
595  do iq=1,nq
596 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
597  do l=1,llm
598  q_cum(:,jjb:jje,l,iq)=q_cum(:,jjb:jje,l,iq)
599  . /masse_cum(:,jjb:jje,l)
600  enddo
601 !$OMP ENDDO NOWAIT
602  enddo
603 
604  zz=1./REAL(ncum)
605 
606 !$OMP MASTER
607  ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
608 !$OMP END MASTER
609 
610 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
611  DO l=1,llm
612  masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
613  flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
614  flux_uq_cum(:,jjb:jje,l,:)=flux_uq_cum(:,jjb:jje,l,:)*zz
615  dq(:,jjb:jje,l,:)=dq(:,jjb:jje,l,:)*zz
616  ENDDO
617 !$OMP ENDDO NOWAIT
618 
619 
620  IF (pole_sud) jje=jj_end-1
621 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
622  DO l=1,llm
623  flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)*zz
624  flux_vq_cum(:,jjb:jje,l,:)=flux_vq_cum(:,jjb:jje,l,:)*zz
625  ENDDO
626 !$OMP ENDDO NOWAIT
627 
628  jjb=jj_begin
629  jje=jj_end
630 
631 
632 c A retravailler eventuellement
633 c division de dQ par la masse pour revenir aux bonnes grandeurs
634  do iq=1,nq
635 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
636  DO l=1,llm
637  dq(:,jjb:jje,l,iq)=dq(:,jjb:jje,l,iq)/masse_cum(:,jjb:jje,l)
638  ENDDO
639 !$OMP ENDDO NOWAIT
640  enddo
641 
642 c=====================================================================
643 c Transport méridien
644 c=====================================================================
645 
646 c cumul zonal des masses des mailles
647 c ----------------------------------
648  jjb=jj_begin
649  jje=jj_end
650  if (pole_sud) jje=jj_end-1
651 
652 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
653  DO l=1,llm
654  zv(jjb:jje,l)=0.
655  zmasse(jjb:jje,l)=0.
656  ENDDO
657 !$OMP ENDDO NOWAIT
658 
659  call register_hallo_u(masse_cum,llm,1,1,1,1,req)
660  do iq=1,nq
661  call register_hallo_u(q_cum(1,1,1,iq),llm,0,1,1,0,req)
662  enddo
663 
664  call sendrequest(req)
665 !$OMP BARRIER
666  call waitrequest(req)
667 
668  call massbar_loc(masse_cum,massebx,masseby)
669 
670  jjb=jj_begin
671  jje=jj_end
672  if (pole_sud) jje=jj_end-1
673 
674 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
675  do l=1,llm
676  do j=jjb,jje
677  do i=1,iim
678  zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
679  zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
680  enddo
681  zfactv(j,l)=cv(1,j)/zmasse(j,l)
682  enddo
683  enddo
684 !$OMP ENDDO NOWAIT
685 
686 c print*,'3OK'
687 c --------------------------------------------------------------
688 c calcul de la moyenne zonale du transport :
689 c ------------------------------------------
690 c
691 c --
692 c TOT : la circulation totale [ vq ]
693 c
694 c - -
695 c MMC : mean meridional circulation [ v ] [ q ]
696 c
697 c ---- -- - -
698 c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
699 c
700 c - * - * - - - -
701 c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
702 c
703 c - -
704 c on utilise aussi l'intermediaire TMP : [ v q ]
705 c
706 c la variable zfactv transforme un transport meridien cumule
707 c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
708 c
709 c --------------------------------------------------------------
710 
711 
712 c ----------------------------------------
713 c Transport dans le plan latitude-altitude
714 c ----------------------------------------
715 
716  jjb=jj_begin
717  jje=jj_end
718  if (pole_sud) jje=jj_end-1
719 
720  zvq=0.
721  psiq=0.
722  do iq=1,nq
723 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
724  do l=1,llm
725  zvqtmp(:,l)=0.
726  do j=jjb,jje
727 c print*,'j,l,iQ=',j,l,iQ
728 c Calcul des moyennes zonales du transort total et de zvQtmp
729  do i=1,iim
730  zvq(j,l,itot,iq)=zvq(j,l,itot,iq)
731  s +flux_vq_cum(i,j,l,iq)
732  zqy= 0.5*(q_cum(i,j,l,iq)*masse_cum(i,j,l)+
733  s q_cum(i,j+1,l,iq)*masse_cum(i,j+1,l))
734  zvqtmp(j,l)=zvqtmp(j,l)+flux_v_cum(i,j,l)*zqy
735  s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
736  zvq(j,l,iave,iq)=zvq(j,l,iave,iq)+zqy
737  enddo
738 c print*,'aOK'
739 c Decomposition
740  zvq(j,l,iave,iq)=zvq(j,l,iave,iq)/zmasse(j,l)
741  zvq(j,l,itot,iq)=zvq(j,l,itot,iq)*zfactv(j,l)
742  zvqtmp(j,l)=zvqtmp(j,l)*zfactv(j,l)
743  zvq(j,l,immc,iq)=zv(j,l)*zvq(j,l,iave,iq)*zfactv(j,l)
744  zvq(j,l,itrs,iq)=zvq(j,l,itot,iq)-zvqtmp(j,l)
745  zvq(j,l,istn,iq)=zvqtmp(j,l)-zvq(j,l,immc,iq)
746  enddo
747  enddo
748 !$OMP ENDDO NOWAIT
749 c fonction de courant meridienne pour la quantite Q
750 !$OMP BARRIER
751 !$OMP MASTER
752  do l=llm,1,-1
753  do j=jjb,jje
754  psiq(j,l,iq)=psiq(j,l+1,iq)+zvq(j,l,itot,iq)
755  enddo
756  enddo
757 !$OMP END MASTER
758 !$OMP BARRIER
759  enddo
760 
761 c fonction de courant pour la circulation meridienne moyenne
762 !$OMP BARRIER
763 !$OMP MASTER
764  psi(jjb:jje,:)=0.
765  do l=llm,1,-1
766  do j=jjb,jje
767  psi(j,l)=psi(j,l+1)+zv(j,l)
768  zv(j,l)=zv(j,l)*zfactv(j,l)
769  enddo
770  enddo
771 !$OMP END MASTER
772 !$OMP BARRIER
773 
774 c print*,'4OK'
775 c sorties proprement dites
776 !$OMP MASTER
777  if (i_sortie.eq.1) then
778  jjb=jj_begin
779  jje=jj_end
780  jjn=jj_nb
781  if (pole_sud) jje=jj_end-1
782  if (pole_sud) jjn=jj_nb-1
783  do iq=1,nq
784  do itr=1,ntr
785  call histwrite(fileid,znom(itr,iq),itau,
786  s zvq(jjb:jje,:,itr,iq)
787  s ,jjn*llm,ndex3d)
788  enddo
789  call histwrite(fileid,'psi'//nom(iq),
790  s itau,psiq(jjb:jje,1:llm,iq)
791  s ,jjn*llm,ndex3d)
792  enddo
793 
794  call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
795  s ,jjn*llm,ndex3d)
796  call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
797  s ,jjn*llm,ndex3d)
798  psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
799  call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
800  s jjn*llm,ndex3d)
801 
802  endif
803 
804 
805 c -----------------
806 c Moyenne verticale
807 c -----------------
808 
809  zamasse(jjb:jje)=0.
810  do l=1,llm
811  zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
812  enddo
813 
814  zavq(jjb:jje,:,:)=0.
815  do iq=1,nq
816  do itr=2,ntr
817  do l=1,llm
818  zavq(jjb:jje,itr,iq)=zavq(jjb:jje,itr,iq)
819  s +zvq(jjb:jje,l,itr,iq)
820  s *zmasse(jjb:jje,l)
821  enddo
822  zavq(jjb:jje,itr,iq)=zavq(jjb:jje,itr,iq)/zamasse(jjb:jje)
823  call histwrite(fileid,'a'//znom(itr,iq),itau,
824  s zavq(jjb:jje,itr,iq),jjn*llm,ndex3d)
825  enddo
826  enddo
827 !$OMP END MASTER
828 c on doit pouvoir tracer systematiquement la fonction de courant.
829 
830 c=====================================================================
831 c/////////////////////////////////////////////////////////////////////
832  icum=0 !///////////////////////////////////////
833  endif ! icum.eq.ncum !///////////////////////////////////////
834 c/////////////////////////////////////////////////////////////////////
835 c=====================================================================
836 
837  return
838  end