source: LMDZ6/branches/Amaury_dev/libf/dyn3dmem/bilan_dyn_loc.f90 @ 5136

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

Put comgeom.h, comgeom2.h into modules

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