source: LMDZ5/trunk/libf/dyn3dpar/bilan_dyn_p.F @ 2573

Last change on this file since 2573 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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