source: LMDZ5/branches/LMDZ_tree_FC/libf/dyn3dmem/bilan_dyn_loc.F

Last change on this file was 2601, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: turn temps.h into module temps_mod.F90
EM

  • 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.8 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_lmdz
16      USE mod_hallo
17      use misc_mod
18      USE write_field_loc
19      USE comconst_mod, ONLY: cpp, pi
20      USE comvert_mod, ONLY: presnivs
21      USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
22     
23      IMPLICIT NONE
24
25      include "dimensions.h"
26      include "paramet.h"
27      include "comgeom2.h"
28      include "iniprint.h"
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
143      integer fileid
144      integer thoriid, zvertiid
145      save fileid
146
147      INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)
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,ALLOCATABLE :: rlong(:),rlatg(:)
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
175c=====================================================================
176c   Initialisation
177c=====================================================================
178      if (adjust) return
179     
180      time=time+dt_app
181      itau=itau+1
182
183      if (first) then
184!$OMP BARRIER
185!$OMP MASTER
186      ALLOCATE(zfactv(jjb_v:jje_v,llm))
187      ALLOCATE(vcont(iip1,jjb_v:jje_v,llm))
188      ALLOCATE(ucont(iip1,jjb_u:jje_u,llm))
189      ALLOCATE(ang(iip1,jjb_u:jje_u,llm))
190      ALLOCATE(unat(iip1,jjb_u:jje_u,llm))
191      ALLOCATE(massebx(iip1,jjb_u:jje_u,llm))
192      ALLOCATE(masseby(iip1,jjb_v:jje_v,llm))
193      ALLOCATE(vorpot(iip1,jjb_v:jje_v,llm))
194      ALLOCATE(w(iip1,jjb_u:jje_u,llm))
195      ALLOCATE(ecin(iip1,jjb_u:jje_u,llm))
196      ALLOCATE(convm(iip1,jjb_u:jje_u,llm))
197      ALLOCATE(bern(iip1,jjb_u:jje_u,llm))     
198      ALLOCATE(Q(iip1,jjb_u:jje_u,llm,nQ))     
199      ALLOCATE(ps_cum(iip1,jjb_u:jje_u))
200      ALLOCATE(masse_cum(iip1,jjb_u:jje_u,llm))
201      ALLOCATE(flux_u_cum(iip1,jjb_u:jje_u,llm))
202      ALLOCATE(flux_v_cum(iip1,jjb_v:jje_v,llm))
203      ALLOCATE(Q_cum(iip1,jjb_u:jje_u,llm,nQ))
204      ALLOCATE(flux_uQ_cum(iip1,jjb_u:jje_u,llm,nQ))
205      ALLOCATE(flux_vQ_cum(iip1,jjb_v:jje_v,llm,nQ))
206      ALLOCATE(flux_wQ_cum(iip1,jjb_u:jje_u,llm,nQ))
207      ALLOCATE(dQ(iip1,jjb_u:jje_u,llm,nQ))
208      ALLOCATE(zvQ(jjb_v:jje_v,llm,ntr,nQ))
209      ALLOCATE(zvQtmp(jjb_v:jje_v,llm))
210      ALLOCATE(zavQ(jjb_v:jje_v,ntr,nQ))
211      ALLOCATE(psiQ(jjb_v:jje_v,llm+1,nQ))
212      ALLOCATE(zmasse(jjb_v:jje_v,llm))
213      ALLOCATE(zamasse(jjb_v:jje_v))
214      ALLOCATE(zv(jjb_v:jje_v,llm))
215      ALLOCATE(psi(jjb_v:jje_v,llm+1))
216      ALLOCATE(ndex3d(jjb_v:jje_v*llm))
217      ndex3d=0
218      ALLOCATE(rlong(1))
219      ALLOCATE(rlatg(jjm))
220     
221!$OMP END MASTER
222!$OMP BARRIER
223        icum=0
224c       initialisation des fichiers
225        first=.false.
226c   ncum est la frequence de stokage en pas de temps
227        ncum=dt_cum/dt_app
228        if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
229           WRITE(lunout,*)
230     .            'Pb : le pas de cumule doit etre multiple du pas'
231           WRITE(lunout,*)'dt_app=',dt_app
232           WRITE(lunout,*)'dt_cum=',dt_cum
233           stop
234        endif
235
236!$OMP MASTER
237        nom(itemp)='T'
238        nom(igeop)='gz'
239        nom(iecin)='K'
240        nom(iang)='ang'
241        nom(iu)='u'
242        nom(iovap)='ovap'
243        nom(iun)='un'
244
245        unites(itemp)='K'
246        unites(igeop)='m2/s2'
247        unites(iecin)='m2/s2'
248        unites(iang)='ang'
249        unites(iu)='m/s'
250        unites(iovap)='kg/kg'
251        unites(iun)='un'
252
253
254c   Initialisation du fichier contenant les moyennes zonales.
255c   ---------------------------------------------------------
256
257      infile='dynzon'
258
259      zan = annee_ref
260      dayref = day_ref
261      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
262      tau0 = itau_dyn
263     
264      rlong=0.
265      rlatg=rlatv*180./pi
266
267      jjb=jj_begin
268      jje=jj_end
269      jjn=jj_nb
270      IF (pole_sud) THEN
271        jjn=jj_nb-1
272        jje=jj_end-1
273      ENDIF
274
275      ddid=(/ 2 /)
276      dsg=(/ jjm /)
277      dsl=(/ jjn /)
278      dpf=(/ jjb /)
279      dpl=(/ jje /)
280      dhs=(/ 0 /)
281      dhe=(/ 0 /)
282
283      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
284     .                 'box',bilan_dyn_domain_id)
285       
286      call histbeg(trim(infile),
287     .             1, rlong, jjn, rlatg(jjb:jje),
288     .             1, 1, 1, jjn,
289     .             tau0, zjulian, dt_cum, thoriid, fileid,
290     .             bilan_dyn_domain_id)
291
292C
293C  Appel a histvert pour la grille verticale
294C
295      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
296     .              llm, presnivs, zvertiid)
297C
298C  Appels a histdef pour la definition des variables a sauvegarder
299      do iQ=1,nQ
300         do itr=1,ntr
301            if(itr.eq.1) then
302               znom(itr,iQ)=nom(iQ)
303               znoml(itr,iQ)=nom(iQ)
304               zunites(itr,iQ)=unites(iQ)
305            else
306               znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
307               znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
308               zunites(itr,iQ)='m/s * '//unites(iQ)
309            endif
310         enddo
311      enddo
312
313c   Declarations des champs avec dimension verticale
314c      print*,'1HISTDEF'
315      do iQ=1,nQ
316         do itr=1,ntr
317      IF (prt_level > 5)
318     . WRITE(lunout,*)'var ',itr,iQ
319     .      ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
320            call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
321     .        zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
322     .        32,'ave(X)',dt_cum,dt_cum)
323         enddo
324c   Declarations pour les fonctions de courant
325c      print*,'2HISTDEF'
326          call histdef(fileid,'psi'//nom(iQ)
327     .      ,'stream fn. '//znoml(itot,iQ),
328     .      zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
329     .      32,'ave(X)',dt_cum,dt_cum)
330      enddo
331
332
333c   Declarations pour les champs de transport d'air
334c      print*,'3HISTDEF'
335      call histdef(fileid, 'masse', 'masse',
336     .             'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
337     .             32, 'ave(X)', dt_cum, dt_cum)
338      call histdef(fileid, 'v', 'v',
339     .             'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
340     .             32, 'ave(X)', dt_cum, dt_cum)
341c   Declarations pour les fonctions de courant
342c      print*,'4HISTDEF'
343          call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
344     .      1,jjn,thoriid,llm,1,llm,zvertiid,
345     .      32,'ave(X)',dt_cum,dt_cum)
346
347
348c   Declaration des champs 1D de transport en latitude
349c      print*,'5HISTDEF'
350      do iQ=1,nQ
351         do itr=2,ntr
352            call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
353     .        zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
354     .        32,'ave(X)',dt_cum,dt_cum)
355         enddo
356      enddo
357
358
359c      print*,'8HISTDEF'
360               CALL histend(fileid)
361
362!$OMP END MASTER
363      endif
364
365
366c=====================================================================
367c   Calcul des champs dynamiques
368c   ----------------------------
369
370      jjb=jj_begin
371      jje=jj_end
372   
373c   énergie cinétique
374!      ucont(:,jjb:jje,:)=0
375
376      call Register_Hallo_u(ucov,llm,1,1,1,1,Req)
377      call Register_Hallo_v(vcov,llm,1,1,1,1,Req)
378      call SendRequest(Req)
379c$OMP BARRIER
380      call WaitRequest(Req)
381
382      CALL covcont_loc(llm,ucov,vcov,ucont,vcont)
383      CALL enercin_loc(vcov,ucov,vcont,ucont,ecin)
384
385c   moment cinétique
386!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
387      do l=1,llm
388         ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
389         unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
390      enddo
391!$OMP END DO NOWAIT
392
393!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
394      DO l=1,llm
395        Q(:,jjb:jje,l,itemp)=teta(:,jjb:jje,l)*pk(:,jjb:jje,l)/cpp
396        Q(:,jjb:jje,l,igeop)=phi(:,jjb:jje,l)
397        Q(:,jjb:jje,l,iecin)=ecin(:,jjb:jje,l)
398        Q(:,jjb:jje,l,iang)=ang(:,jjb:jje,l)
399        Q(:,jjb:jje,l,iu)=unat(:,jjb:jje,l)
400        Q(:,jjb:jje,l,iovap)=trac(:,jjb:jje,l,1)
401        Q(:,jjb:jje,l,iun)=1.
402      ENDDO
403!$OMP END DO NOWAIT
404
405c=====================================================================
406c   Cumul
407c=====================================================================
408c
409      if(icum.EQ.0) then
410         jjb=jj_begin
411         jje=jj_end
412
413!$OMP MASTER
414         ps_cum(:,jjb:jje)=0.
415!$OMP END MASTER
416
417
418!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
419        DO l=1,llm
420          masse_cum(:,jjb:jje,l)=0.
421          flux_u_cum(:,jjb:jje,l)=0.
422          Q_cum(:,jjb:jje,l,:)=0.
423          flux_uQ_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,l)=masse_cum(:,jjb:jje,l)+masse(:,jjb:jje,l)
447        flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)
448     .                         +flux_u(:,jjb:jje,l)
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,l)=flux_v_cum(:,jjb:jje,l)
457     .                          +flux_v(:,jjb:jje,l)
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,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
468     .                       +Q(:,jjb:jje,l,iQ)*masse(:,jjb:jje,l)
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,jjb_u,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 ENDDO NOWAIT
517!$OMP BARRIER
518      enddo
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
544      jjb=jj_begin
545      jje=jj_end
546
547!      do iQ=1,nQ
548!         do l=1,llm-1
549!            do j=jjb,jje
550!               do i=1,iip1
551!                  ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
552!                  dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
553!                  dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
554!               enddo
555!            enddo
556!          enddo
557!       enddo
558       
559      do iQ=1,nQ
560!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
561         do l=1,llm
562            IF (l<llm) THEN
563              do j=jjb,jje
564                 do i=1,iip1
565                    ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
566                    dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
567                    dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
568                 enddo
569              enddo
570            ENDIF
571            IF (l>2) THEN
572              do j=jjb,jje
573                do i=1,iip1
574                  ww=-0.5*w(i,j,l)*(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
575                  dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
576                enddo
577              enddo
578            ENDIF
579         enddo
580!$OMP ENDDO NOWAIT
581      enddo
582      IF (prt_level > 5)
583     . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
584c=====================================================================
585c   PAS DE TEMPS D'ECRITURE
586c=====================================================================
587      if (icum.eq.ncum) then
588c=====================================================================
589
590      IF (prt_level > 5)
591     . WRITE(lunout,*)'Pas d ecriture'
592
593      jjb=jj_begin
594      jje=jj_end
595
596c   Normalisation
597      do iQ=1,nQ
598!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
599        do l=1,llm
600          Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
601     .                                /masse_cum(:,jjb:jje,l)
602        enddo
603!$OMP ENDDO NOWAIT
604      enddo   
605
606      zz=1./REAL(ncum)
607
608!$OMP MASTER
609        ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
610!$OMP END MASTER
611
612!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
613      DO l=1,llm
614        masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
615        flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
616        flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
617        dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
618      ENDDO
619!$OMP ENDDO NOWAIT
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!$OMP BARRIER
629         
630      jjb=jj_begin
631      jje=jj_end
632
633
634c   A retravailler eventuellement
635c   division de dQ par la masse pour revenir aux bonnes grandeurs
636      do iQ=1,nQ
637!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
638        DO l=1,llm
639           dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
640        ENDDO
641!$OMP ENDDO NOWAIT
642      enddo
643
644c=====================================================================
645c   Transport méridien
646c=====================================================================
647
648c   cumul zonal des masses des mailles
649c   ----------------------------------
650      jjb=jj_begin
651      jje=jj_end
652      if (pole_sud) jje=jj_end-1
653
654!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
655        DO l=1,llm
656          zv(jjb:jje,l)=0.
657          zmasse(jjb:jje,l)=0.
658        ENDDO
659!$OMP ENDDO NOWAIT
660!$OMP BARRIER
661
662      call Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
663      do iQ=1,nQ
664        call Register_Hallo_u(Q_cum(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
665      enddo
666
667      call SendRequest(Req)
668!$OMP BARRIER
669      call WaitRequest(Req)
670
671      call massbar_loc(masse_cum,massebx,masseby)
672     
673      jjb=jj_begin
674      jje=jj_end
675      if (pole_sud) jje=jj_end-1
676     
677!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
678      do l=1,llm
679         do j=jjb,jje
680            do i=1,iim
681               zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
682               zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
683            enddo
684            zfactv(j,l)=cv(1,j)/zmasse(j,l)
685         enddo
686      enddo
687!$OMP ENDDO NOWAIT
688!$OMP BARRIER
689
690c     print*,'3OK'
691c   --------------------------------------------------------------
692c   calcul de la moyenne zonale du transport :
693c   ------------------------------------------
694c
695c                                     --
696c TOT : la circulation totale       [ vq ]
697c
698c                                      -     -
699c MMC : mean meridional circulation [ v ] [ q ]
700c
701c                                     ----      --       - -
702c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
703c
704c                                     - * - *       - -       -     -
705c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
706c
707c                                              - -
708c    on utilise aussi l'intermediaire TMP :  [ v q ]
709c
710c    la variable zfactv transforme un transport meridien cumule
711c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
712c
713c   --------------------------------------------------------------
714
715
716c   ----------------------------------------
717c   Transport dans le plan latitude-altitude
718c   ----------------------------------------
719
720      jjb=jj_begin
721      jje=jj_end
722      if (pole_sud) jje=jj_end-1
723     
724      zvQ=0.
725      psiQ=0.
726      do iQ=1,nQ
727!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
728         do l=1,llm
729            zvQtmp(:,l)=0.
730            do j=jjb,jje
731c              print*,'j,l,iQ=',j,l,iQ
732c   Calcul des moyennes zonales du transort total et de zvQtmp
733               do i=1,iim
734                  zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
735     s                            +flux_vQ_cum(i,j,l,iQ)
736                  zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
737     s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
738                  zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
739     s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
740                  zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
741               enddo
742c              print*,'aOK'
743c   Decomposition
744               zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
745               zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
746               zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
747               zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
748               zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
749               zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
750            enddo
751         enddo
752!$OMP ENDDO NOWAIT
753c   fonction de courant meridienne pour la quantite Q
754!$OMP BARRIER
755!$OMP MASTER
756         do l=llm,1,-1
757            do j=jjb,jje
758               psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
759            enddo
760         enddo
761!$OMP END MASTER
762!$OMP BARRIER
763      enddo
764
765c   fonction de courant pour la circulation meridienne moyenne
766!$OMP BARRIER
767!$OMP MASTER
768      psi(jjb:jje,:)=0.
769      do l=llm,1,-1
770         do j=jjb,jje
771            psi(j,l)=psi(j,l+1)+zv(j,l)
772            zv(j,l)=zv(j,l)*zfactv(j,l)
773         enddo
774      enddo
775!$OMP END MASTER
776!$OMP BARRIER
777
778c     print*,'4OK'
779c   sorties proprement dites
780!$OMP MASTER     
781      if (i_sortie.eq.1) then
782      jjb=jj_begin
783      jje=jj_end
784      jjn=jj_nb
785      if (pole_sud) jje=jj_end-1
786      if (pole_sud) jjn=jj_nb-1
787      do iQ=1,nQ
788         do itr=1,ntr
789            call histwrite(fileid,znom(itr,iQ),itau,
790     s                     zvQ(jjb:jje,:,itr,iQ)
791     s                     ,jjn*llm,ndex3d)
792         enddo
793         call histwrite(fileid,'psi'//nom(iQ),
794     s                  itau,psiQ(jjb:jje,1:llm,iQ)
795     s                  ,jjn*llm,ndex3d)
796      enddo
797
798      call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
799     s   ,jjn*llm,ndex3d)
800      call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
801     s   ,jjn*llm,ndex3d)
802      psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
803      call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
804     s               jjn*llm,ndex3d)
805
806      endif
807
808 
809c   -----------------
810c   Moyenne verticale
811c   -----------------
812
813      zamasse(jjb:jje)=0.
814      do l=1,llm
815         zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
816      enddo
817     
818      zavQ(jjb:jje,:,:)=0.
819      do iQ=1,nQ
820         do itr=2,ntr
821            do l=1,llm
822               zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
823     s                             +zvQ(jjb:jje,l,itr,iQ)
824     s                             *zmasse(jjb:jje,l)
825            enddo
826            zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
827            call histwrite(fileid,'a'//znom(itr,iQ),itau,
828     s                     zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
829         enddo
830      enddo
831!$OMP END MASTER
832c     on doit pouvoir tracer systematiquement la fonction de courant.
833
834c=====================================================================
835c/////////////////////////////////////////////////////////////////////
836      icum=0                  !///////////////////////////////////////
837      endif ! icum.eq.ncum    !///////////////////////////////////////
838c/////////////////////////////////////////////////////////////////////
839c=====================================================================
840
841      return
842      end
Note: See TracBrowser for help on using the repository browser.