source: LMDZ6/trunk/libf/dyn3dmem/bilan_dyn_loc.f90 @ 5272

Last change on this file since 5272 was 5272, checked in by abarral, 3 months ago

Turn paramet.h into a module

  • 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: 22.1 KB
Line 
1!
2! $Id: bilan_dyn_p.F 1299 2010-01-20 14:27:21Z fairhead $
3!
4SUBROUTINE bilan_dyn_loc (ntrac,dt_app,dt_cum, &
5        ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
6
7  !   AFAIRE
8  !   Prevoir en champ nq+1 le diagnostique de l'energie
9  !   en faisant Qzon=Cv T + L * ...
10  !             vQ..A=Cp T + L * ...
11
12  USE IOIPSL
13  USE parallel_lmdz
14  USE mod_hallo
15  use misc_mod
16  USE write_field_loc
17  USE comconst_mod, ONLY: cpp, pi
18  USE comvert_mod, ONLY: presnivs
19  USE temps_mod, ONLY: annee_ref, day_ref, itau_dyn
20
21  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
22USE paramet_mod_h, ONLY: iip1, iip2, iip3, jjp1, llmp1, llmp2, llmm1, kftd, ip1jm, ip1jmp1, &
23          ip1jmi1, ijp1llm, ijmllm, mvar, jcfil, jcfllm
24IMPLICIT NONE
25
26
27
28  include "comgeom2.h"
29  include "iniprint.h"
30
31  !====================================================================
32  !
33  !   Sous-programme consacre � des diagnostics dynamiques de base
34  !
35  !
36  !   De facon generale, les moyennes des scalaires Q sont ponderees par
37  !   la masse.
38  !
39  !   Les flux de masse sont eux simplement moyennes.
40  !
41  !====================================================================
42
43  !   Arguments :
44  !   ===========
45
46  integer :: ntrac
47  real :: dt_app,dt_cum
48  real :: ps(iip1,jjb_u:jje_u)
49  real :: masse(iip1,jjb_u:jje_u,llm),pk(iip1,jjb_u:jje_u,llm)
50  real :: flux_u(iip1,jjb_u:jje_u,llm)
51  real :: flux_v(iip1,jjb_v:jje_v,llm)
52  real :: teta(iip1,jjb_u:jje_u,llm)
53  real :: phi(iip1,jjb_u:jje_u,llm)
54  real :: ucov(iip1,jjb_u:jje_u,llm)
55  real :: vcov(iip1,jjb_v:jje_v,llm)
56  real :: trac(iip1,jjb_u:jje_u,llm,ntrac)
57
58  !   Local :
59  !   =======
60
61  integer,SAVE :: icum,ncum
62!$OMP THREADPRIVATE(icum,ncum)
63  LOGICAL,SAVE :: first=.TRUE.
64!$OMP THREADPRIVATE(first)
65
66  real :: zz,zqy
67  REAl,SAVE,ALLOCATABLE :: zfactv(:,:)
68
69  INTEGER,PARAMETER :: nQ=7
70
71
72  !ym      character*6 nom(nQ)
73  !ym      character*6 unites(nQ)
74  character(len=6),save :: nom(nQ)
75  character(len=6),save :: unites(nQ)
76
77  character(len=10) file
78  integer :: ifile
79  parameter (ifile=4)
80
81  integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
82  INTEGER,PARAMETER :: iovap=6,iun=7
83  integer,PARAMETER :: i_sortie=1
84
85  real,SAVE :: time=0.
86  integer,SAVE :: itau=0.
87!$OMP THREADPRIVATE(time,itau)
88
89  real :: ww
90
91  !   variables dynamiques interm�diaires
92  REAL,SAVE,ALLOCATABLE :: vcont(:,:,:),ucont(:,:,:)
93  REAL,SAVE,ALLOCATABLE :: ang(:,:,:),unat(:,:,:)
94  REAL,SAVE,ALLOCATABLE :: massebx(:,:,:),masseby(:,:,:)
95  REAL,SAVE,ALLOCATABLE :: vorpot(:,:,:)
96  REAL,SAVE,ALLOCATABLE :: w(:,:,:),ecin(:,:,:),convm(:,:,:)
97  REAL,SAVE,ALLOCATABLE :: bern(:,:,:)
98
99  !   champ contenant les scalaires advect�s.
100  real,SAVE,ALLOCATABLE :: Q(:,:,:,:)
101
102  !   champs cumul�s
103  real,SAVE,ALLOCATABLE ::  ps_cum(:,:)
104  real,SAVE,ALLOCATABLE ::  masse_cum(:,:,:)
105  real,SAVE,ALLOCATABLE ::  flux_u_cum(:,:,:)
106  real,SAVE,ALLOCATABLE ::  flux_v_cum(:,:,:)
107  real,SAVE,ALLOCATABLE ::  Q_cum(:,:,:,:)
108  real,SAVE,ALLOCATABLE ::  flux_uQ_cum(:,:,:,:)
109  real,SAVE,ALLOCATABLE ::  flux_vQ_cum(:,:,:,:)
110  real,SAVE,ALLOCATABLE ::  flux_wQ_cum(:,:,:,:)
111  real,SAVE,ALLOCATABLE ::  dQ(:,:,:,:)
112
113
114  !   champs de tansport en moyenne zonale
115  integer :: ntr,itr
116  parameter (ntr=5)
117
118  !ym      character*10 znom(ntr,nQ)
119  !ym      character*20 znoml(ntr,nQ)
120  !ym      character*10 zunites(ntr,nQ)
121  character*10,save :: znom(ntr,nQ)
122  character*20,save :: znoml(ntr,nQ)
123  character*10,save :: zunites(ntr,nQ)
124
125  INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5
126
127  character(len=3) :: ctrs(ntr)
128  data ctrs/'  ','TOT','MMC','TRS','STN'/
129
130  real,SAVE,ALLOCATABLE ::  zvQ(:,:,:,:),zvQtmp(:,:)
131  real,SAVE,ALLOCATABLE ::  zavQ(:,:,:),psiQ(:,:,:)
132  real,SAVE,ALLOCATABLE ::  zmasse(:,:),zamasse(:)
133
134  real,SAVE,ALLOCATABLE ::  zv(:,:),psi(:,:)
135
136  integer :: i,j,l,iQ
137
138
139  !   Initialisation du fichier contenant les moyennes zonales.
140  !   ---------------------------------------------------------
141
142  character(len=10) :: infile
143
144  integer, save :: fileid
145  integer :: thoriid, zvertiid
146
147  INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)
148
149  !   Variables locales
150  !
151  integer :: tau0
152  real :: zjulian
153  character(len=3) :: str
154  character(len=10) :: ctrac
155  integer :: ii,jj
156  integer :: zan, dayref
157  !
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
175  !=====================================================================
176  !   Initialisation
177  !=====================================================================
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
224    ! initialisation des fichiers
225    first=.false.
226  !   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       CALL abort_gcm("conf_gcmbilan_dyn_loc","stopped",1)
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
254  !   Initialisation du fichier contenant les moyennes zonales.
255  !   ---------------------------------------------------------
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
292  !
293  !  Appel a histvert pour la grille verticale
294  !
295  call histvert(fileid, 'presnivs', 'Niveaux sigma','mb', &
296        llm, presnivs, zvertiid)
297  !
298  !  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
313  !   Declarations des champs avec dimension verticale
314   ! 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
324  !   Declarations pour les fonctions de courant
325   ! 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
333  !   Declarations pour les champs de transport d'air
334   ! 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)
341  !   Declarations pour les fonctions de courant
342   ! 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
348  !   Declaration des champs 1D de transport en latitude
349   ! 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
359   ! print*,'8HISTDEF'
360           CALL histend(fileid)
361
362!$OMP END MASTER
363  endif
364
365
366  !=====================================================================
367  !   Calcul des champs dynamiques
368  !   ----------------------------
369
370  jjb=jj_begin
371  jje=jj_end
372
373  !   �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)
379!$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
385  !   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
405  !=====================================================================
406  !   Cumul
407  !=====================================================================
408  !
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
435  !   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
473  !=====================================================================
474  !  FLUX ET TENDANCES
475  !=====================================================================
476
477  !   Flux longitudinal
478  !   -----------------
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                    +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
493  !    flux m�ridien
494  !    -------------
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                    +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
520  !    tendances
521  !    ---------
522
523  !   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
532  !   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'
584  !=====================================================================
585  !   PAS DE TEMPS D'ECRITURE
586  !=====================================================================
587  if (icum.eq.ncum) then
588  !=====================================================================
589
590  IF (prt_level > 5) &
591        WRITE(lunout,*)'Pas d ecriture'
592
593  jjb=jj_begin
594  jje=jj_end
595
596  !   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
634  !   A retravailler eventuellement
635  !   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
644  !=====================================================================
645  !   Transport m�ridien
646  !=====================================================================
647
648  !   cumul zonal des masses des mailles
649  !   ----------------------------------
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
690  ! print*,'3OK'
691  !   --------------------------------------------------------------
692  !   calcul de la moyenne zonale du transport :
693  !   ------------------------------------------
694  !
695  !                                 --
696  ! TOT : la circulation totale       [ vq ]
697  !
698  !                                  -     -
699  ! MMC : mean meridional circulation [ v ] [ q ]
700  !
701  !                                 ----      --       - -
702  ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
703  !
704  !                                 - * - *       - -       -     -
705  ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
706  !
707  !                                          - -
708  !    on utilise aussi l'intermediaire TMP :  [ v q ]
709  !
710  !    la variable zfactv transforme un transport meridien cumule
711  !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
712  !
713  !   --------------------------------------------------------------
714
715
716  !   ----------------------------------------
717  !   Transport dans le plan latitude-altitude
718  !   ----------------------------------------
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
731           ! print*,'j,l,iQ=',j,l,iQ
732  !   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                    +flux_vQ_cum(i,j,l,iQ)
736              zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ &
737                    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                    /(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
742           ! print*,'aOK'
743  !   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
753  !   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
765  !   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
778  ! print*,'4OK'
779  !   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              zvQ(jjb:jje,:,itr,iQ) &
791              ,jjn*llm,ndex3d)
792     enddo
793     call histwrite(fileid,'psi'//nom(iQ), &
794           itau,psiQ(jjb:jje,1:llm,iQ) &
795           ,jjn*llm,ndex3d)
796  enddo
797
798  call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm) &
799        ,jjn*llm,ndex3d)
800  call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm) &
801        ,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        jjn*llm,ndex3d)
805
806  endif
807
808
809  !   -----------------
810  !   Moyenne verticale
811  !   -----------------
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                 +zvQ(jjb:jje,l,itr,iQ) &
824                 *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              zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
829     enddo
830  enddo
831!$OMP END MASTER
832  ! on doit pouvoir tracer systematiquement la fonction de courant.
833
834  !=====================================================================
835  !/////////////////////////////////////////////////////////////////////
836  icum=0                  !///////////////////////////////////////
837  endif ! icum.eq.ncum    !///////////////////////////////////////
838  !/////////////////////////////////////////////////////////////////////
839  !=====================================================================
840!$OMP MASTER
841  call histsync(fileid)
842!$OMP END MASTER
843
844
845  return
846end subroutine bilan_dyn_loc
Note: See TracBrowser for help on using the repository browser.