source: LMDZ6/trunk/libf/dyn3dmem/bilan_dyn_loc.F @ 5069

Last change on this file since 5069 was 4605, checked in by fhourdin, 17 months ago

Synchronisation des dynzon.nc a chaque ecriture

  • 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
File size: 23.9 KB
RevLine 
[1632]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
[1823]15      USE parallel_lmdz
[1632]16      USE mod_hallo
17      use misc_mod
[2475]18      USE write_field_loc
[2597]19      USE comconst_mod, ONLY: cpp, pi
[2600]20      USE comvert_mod, ONLY: presnivs
[2601]21      USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
[2600]22     
[1632]23      IMPLICIT NONE
24
[2597]25      include "dimensions.h"
26      include "paramet.h"
27      include "comgeom2.h"
28      include "iniprint.h"
[1632]29
30c====================================================================
31c
32c   Sous-programme consacre à des diagnostics dynamiques de base
33c
34c
35c   De facon generale, les moyennes des scalaires Q sont ponderees par
36c   la masse.
37c
38c   Les flux de masse sont eux simplement moyennes.
39c
40c====================================================================
41
42c   Arguments :
43c   ===========
44
45      integer ntrac
46      real dt_app,dt_cum
47      real ps(iip1,jjb_u:jje_u)
48      real masse(iip1,jjb_u:jje_u,llm),pk(iip1,jjb_u:jje_u,llm)
49      real flux_u(iip1,jjb_u:jje_u,llm)
50      real flux_v(iip1,jjb_v:jje_v,llm)
51      real teta(iip1,jjb_u:jje_u,llm)
52      real phi(iip1,jjb_u:jje_u,llm)
53      real ucov(iip1,jjb_u:jje_u,llm)
54      real vcov(iip1,jjb_v:jje_v,llm)
55      real trac(iip1,jjb_u:jje_u,llm,ntrac)
56
57c   Local :
58c   =======
59
60      integer,SAVE :: icum,ncum
61!$OMP THREADPRIVATE(icum,ncum)
62      LOGICAL,SAVE :: first=.TRUE.
63!$OMP THREADPRIVATE(first)     
64     
65      real zz,zqy
66      REAl,SAVE,ALLOCATABLE :: zfactv(:,:)
67
68      INTEGER,PARAMETER :: nQ=7
69
70
71cym      character*6 nom(nQ)
72cym      character*6 unites(nQ)
73      character(len=6),save :: nom(nQ)
74      character(len=6),save :: unites(nQ)
75
76      character(len=10) file
77      integer ifile
78      parameter (ifile=4)
79
80      integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
81      INTEGER,PARAMETER :: iovap=6,iun=7
82      integer,PARAMETER :: i_sortie=1
83
84      real,SAVE :: time=0.
85      integer,SAVE :: itau=0.
86!$OMP THREADPRIVATE(time,itau)
87
88      real ww
89
90c   variables dynamiques intermédiaires
91      REAL,SAVE,ALLOCATABLE :: vcont(:,:,:),ucont(:,:,:)
92      REAL,SAVE,ALLOCATABLE :: ang(:,:,:),unat(:,:,:)
93      REAL,SAVE,ALLOCATABLE :: massebx(:,:,:),masseby(:,:,:)
94      REAL,SAVE,ALLOCATABLE :: vorpot(:,:,:)
95      REAL,SAVE,ALLOCATABLE :: w(:,:,:),ecin(:,:,:),convm(:,:,:)
96      REAL,SAVE,ALLOCATABLE :: bern(:,:,:)
97
98c   champ contenant les scalaires advectés.
99      real,SAVE,ALLOCATABLE :: Q(:,:,:,:)
100   
101c   champs cumulés
102      real,SAVE,ALLOCATABLE ::  ps_cum(:,:)
103      real,SAVE,ALLOCATABLE ::  masse_cum(:,:,:)
104      real,SAVE,ALLOCATABLE ::  flux_u_cum(:,:,:)
105      real,SAVE,ALLOCATABLE ::  flux_v_cum(:,:,:)
106      real,SAVE,ALLOCATABLE ::  Q_cum(:,:,:,:)
107      real,SAVE,ALLOCATABLE ::  flux_uQ_cum(:,:,:,:)
108      real,SAVE,ALLOCATABLE ::  flux_vQ_cum(:,:,:,:)
109      real,SAVE,ALLOCATABLE ::  flux_wQ_cum(:,:,:,:)
110      real,SAVE,ALLOCATABLE ::  dQ(:,:,:,:)
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,ALLOCATABLE ::  zvQ(:,:,:,:),zvQtmp(:,:)
130      real,SAVE,ALLOCATABLE ::  zavQ(:,:,:),psiQ(:,:,:)
131      real,SAVE,ALLOCATABLE ::  zmasse(:,:),zamasse(:)
132
133      real,SAVE,ALLOCATABLE ::  zv(:,:),psi(:,:)
134
135      integer i,j,l,iQ
136
137
138c   Initialisation du fichier contenant les moyennes zonales.
139c   ---------------------------------------------------------
140
141      character*10 infile
142
[4605]143      integer, save :: fileid
[1632]144      integer thoriid, zvertiid
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
[1848]159      type(Request),SAVE :: Req
160!$OMP THREADPRIVATE(Req)
[1632]161
162! definition du domaine d'ecriture pour le rebuild
163
164      INTEGER,DIMENSION(1) :: ddid
165      INTEGER,DIMENSION(1) :: dsg
166      INTEGER,DIMENSION(1) :: dsl
167      INTEGER,DIMENSION(1) :: dpf
168      INTEGER,DIMENSION(1) :: dpl
169      INTEGER,DIMENSION(1) :: dhs
170      INTEGER,DIMENSION(1) :: dhe
171     
172      INTEGER :: bilan_dyn_domain_id
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
[2475]217      ALLOCATE(rlong(1))
218      ALLOCATE(rlatg(jjm))
[1632]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
[4469]232           CALL abort_gcm("conf_gcmbilan_dyn_loc","stopped",1)
[1632]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),
[2475]286     .             1, rlong, jjn, rlatg(jjb:jje),
[1632]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.
[1874]421          Q_cum(:,jjb:jje,l,:)=0.
[1632]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
434c   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
[1885]445        masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)+masse(:,jjb:jje,l)
446        flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)
447     .                         +flux_u(:,jjb:jje,l)
[1632]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
[1885]455       flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)
456     .                          +flux_v(:,jjb:jje,l)
[1632]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
[1885]466          Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
467     .                       +Q(:,jjb:jje,l,iQ)*masse(:,jjb:jje,l)
[1632]468        ENDDO
469!$OMP END DO NOWAIT
470      enddo
471
472c=====================================================================
473c  FLUX ET TENDANCES
474c=====================================================================
475
476c   Flux longitudinal
477c   -----------------
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
492c    flux méridien
493c    -------------
494      do iQ=1,nQ
[1874]495        call Register_Hallo_u(Q(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
[1632]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
[2475]515!$OMP ENDDO NOWAIT
516!$OMP BARRIER
[1632]517      enddo
518
519c    tendances
520c    ---------
521
522c   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
531c   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
[2475]542
[1632]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)
[2600]600     .                                /masse_cum(:,jjb:jje,l)
[1632]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      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
[2475]627!$OMP BARRIER
[1632]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
[2475]642
[1632]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
[2475]659!$OMP BARRIER
[1632]660
661      call Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
662      do iQ=1,nQ
[1874]663        call Register_Hallo_u(Q_cum(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
[1632]664      enddo
665
666      call SendRequest(Req)
667!$OMP BARRIER
668      call WaitRequest(Req)
669
670      call massbar_loc(masse_cum,massebx,masseby)
671     
672      jjb=jj_begin
673      jje=jj_end
674      if (pole_sud) jje=jj_end-1
675     
676!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
677      do l=1,llm
678         do j=jjb,jje
679            do i=1,iim
680               zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
681               zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
682            enddo
683            zfactv(j,l)=cv(1,j)/zmasse(j,l)
684         enddo
685      enddo
686!$OMP ENDDO NOWAIT
[2475]687!$OMP BARRIER
[1632]688
689c     print*,'3OK'
690c   --------------------------------------------------------------
691c   calcul de la moyenne zonale du transport :
692c   ------------------------------------------
693c
694c                                     --
695c TOT : la circulation totale       [ vq ]
696c
697c                                      -     -
698c MMC : mean meridional circulation [ v ] [ q ]
699c
700c                                     ----      --       - -
701c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
702c
703c                                     - * - *       - -       -     -
704c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
705c
706c                                              - -
707c    on utilise aussi l'intermediaire TMP :  [ v q ]
708c
709c    la variable zfactv transforme un transport meridien cumule
710c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
711c
712c   --------------------------------------------------------------
713
714
715c   ----------------------------------------
716c   Transport dans le plan latitude-altitude
717c   ----------------------------------------
718
719      jjb=jj_begin
720      jje=jj_end
721      if (pole_sud) jje=jj_end-1
722     
723      zvQ=0.
724      psiQ=0.
725      do iQ=1,nQ
726!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
727         do l=1,llm
728            zvQtmp(:,l)=0.
729            do j=jjb,jje
730c              print*,'j,l,iQ=',j,l,iQ
731c   Calcul des moyennes zonales du transort total et de zvQtmp
732               do i=1,iim
733                  zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
734     s                            +flux_vQ_cum(i,j,l,iQ)
735                  zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
736     s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
737                  zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
738     s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
739                  zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
740               enddo
741c              print*,'aOK'
742c   Decomposition
743               zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
744               zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
745               zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
746               zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
747               zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
748               zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
749            enddo
750         enddo
751!$OMP ENDDO NOWAIT
752c   fonction de courant meridienne pour la quantite Q
753!$OMP BARRIER
754!$OMP MASTER
755         do l=llm,1,-1
756            do j=jjb,jje
757               psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
758            enddo
759         enddo
760!$OMP END MASTER
761!$OMP BARRIER
762      enddo
763
764c   fonction de courant pour la circulation meridienne moyenne
765!$OMP BARRIER
766!$OMP MASTER
767      psi(jjb:jje,:)=0.
768      do l=llm,1,-1
769         do j=jjb,jje
770            psi(j,l)=psi(j,l+1)+zv(j,l)
771            zv(j,l)=zv(j,l)*zfactv(j,l)
772         enddo
773      enddo
774!$OMP END MASTER
775!$OMP BARRIER
776
777c     print*,'4OK'
778c   sorties proprement dites
779!$OMP MASTER     
780      if (i_sortie.eq.1) then
781      jjb=jj_begin
782      jje=jj_end
783      jjn=jj_nb
784      if (pole_sud) jje=jj_end-1
785      if (pole_sud) jjn=jj_nb-1
786      do iQ=1,nQ
787         do itr=1,ntr
788            call histwrite(fileid,znom(itr,iQ),itau,
789     s                     zvQ(jjb:jje,:,itr,iQ)
790     s                     ,jjn*llm,ndex3d)
791         enddo
792         call histwrite(fileid,'psi'//nom(iQ),
793     s                  itau,psiQ(jjb:jje,1:llm,iQ)
794     s                  ,jjn*llm,ndex3d)
795      enddo
796
797      call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
798     s   ,jjn*llm,ndex3d)
799      call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
800     s   ,jjn*llm,ndex3d)
801      psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
802      call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
803     s               jjn*llm,ndex3d)
804
805      endif
806
807 
808c   -----------------
809c   Moyenne verticale
810c   -----------------
811
812      zamasse(jjb:jje)=0.
813      do l=1,llm
814         zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
815      enddo
816     
817      zavQ(jjb:jje,:,:)=0.
818      do iQ=1,nQ
819         do itr=2,ntr
820            do l=1,llm
821               zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
822     s                             +zvQ(jjb:jje,l,itr,iQ)
823     s                             *zmasse(jjb:jje,l)
824            enddo
825            zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
826            call histwrite(fileid,'a'//znom(itr,iQ),itau,
827     s                     zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
828         enddo
829      enddo
830!$OMP END MASTER
831c     on doit pouvoir tracer systematiquement la fonction de courant.
832
833c=====================================================================
834c/////////////////////////////////////////////////////////////////////
835      icum=0                  !///////////////////////////////////////
836      endif ! icum.eq.ncum    !///////////////////////////////////////
837c/////////////////////////////////////////////////////////////////////
838c=====================================================================
[4605]839!$OMP MASTER
840      call histsync(fileid)
841!$OMP END MASTER
[1632]842
[4605]843
[1632]844      return
845      end
Note: See TracBrowser for help on using the repository browser.