source: LMDZ6/branches/contrails/libf/dyn3dmem/bilan_dyn_loc.f90 @ 5445

Last change on this file since 5445 was 5285, checked in by abarral, 8 weeks ago

As discussed internally, remove generic ONLY: ... for new _mod_h 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: 22.0 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 iniprint_mod_h
13  USE comgeom2_mod_h
14  USE IOIPSL
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  USE dimensions_mod, ONLY: iim, jjm, llm, ndm
24USE paramet_mod_h
25IMPLICIT NONE
26
27
28
29
30  !====================================================================
31  !
32  !   Sous-programme consacre � des diagnostics dynamiques de base
33  !
34  !
35  !   De facon generale, les moyennes des scalaires Q sont ponderees par
36  !   la masse.
37  !
38  !   Les flux de masse sont eux simplement moyennes.
39  !
40  !====================================================================
41
42  !   Arguments :
43  !   ===========
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
57  !   Local :
58  !   =======
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
71  !ym      character*6 nom(nQ)
72  !ym      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
90  !   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
98  !   champ contenant les scalaires advect�s.
99  real,SAVE,ALLOCATABLE :: Q(:,:,:,:)
100
101  !   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
113  !   champs de tansport en moyenne zonale
114  integer :: ntr,itr
115  parameter (ntr=5)
116
117  !ym      character*10 znom(ntr,nQ)
118  !ym      character*20 znoml(ntr,nQ)
119  !ym      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(len=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
138  !   Initialisation du fichier contenant les moyennes zonales.
139  !   ---------------------------------------------------------
140
141  character(len=10) :: infile
142
143  integer, save :: fileid
144  integer :: thoriid, zvertiid
145
146  INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)
147
148  !   Variables locales
149  !
150  integer :: tau0
151  real :: zjulian
152  character(len=3) :: str
153  character(len=10) :: ctrac
154  integer :: ii,jj
155  integer :: zan, dayref
156  !
157  real,SAVE,ALLOCATABLE :: rlong(:),rlatg(:)
158  integer :: jjb,jje,jjn,ijb,ije
159  type(Request),SAVE :: Req
160!$OMP THREADPRIVATE(Req)
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
174  !=====================================================================
175  !   Initialisation
176  !=====================================================================
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
217  ALLOCATE(rlong(1))
218  ALLOCATE(rlatg(jjm))
219
220!$OMP END MASTER
221!$OMP BARRIER
222    icum=0
223    ! initialisation des fichiers
224    first=.false.
225  !   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
232       CALL abort_gcm("conf_gcmbilan_dyn_loc","stopped",1)
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
253  !   Initialisation du fichier contenant les moyennes zonales.
254  !   ---------------------------------------------------------
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), &
286        1, rlong, jjn, rlatg(jjb:jje), &
287        1, 1, 1, jjn, &
288        tau0, zjulian, dt_cum, thoriid, fileid, &
289        bilan_dyn_domain_id)
290
291  !
292  !  Appel a histvert pour la grille verticale
293  !
294  call histvert(fileid, 'presnivs', 'Niveaux sigma','mb', &
295        llm, presnivs, zvertiid)
296  !
297  !  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
312  !   Declarations des champs avec dimension verticale
313   ! 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
323  !   Declarations pour les fonctions de courant
324   ! 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
332  !   Declarations pour les champs de transport d'air
333   ! 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)
340  !   Declarations pour les fonctions de courant
341   ! 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
347  !   Declaration des champs 1D de transport en latitude
348   ! 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
358   ! print*,'8HISTDEF'
359           CALL histend(fileid)
360
361!$OMP END MASTER
362  endif
363
364
365  !=====================================================================
366  !   Calcul des champs dynamiques
367  !   ----------------------------
368
369  jjb=jj_begin
370  jje=jj_end
371
372  !   �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)
378!$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
384  !   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
404  !=====================================================================
405  !   Cumul
406  !=====================================================================
407  !
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.
421      Q_cum(:,jjb:jje,l,:)=0.
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
434  !   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
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)
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
455   flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l) &
456         +flux_v(:,jjb:jje,l)
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
466      Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ) &
467            +Q(:,jjb:jje,l,iQ)*masse(:,jjb:jje,l)
468    ENDDO
469!$OMP END DO NOWAIT
470  enddo
471
472  !=====================================================================
473  !  FLUX ET TENDANCES
474  !=====================================================================
475
476  !   Flux longitudinal
477  !   -----------------
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                    +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
492  !    flux m�ridien
493  !    -------------
494  do iQ=1,nQ
495    call Register_Hallo_u(Q(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
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                    +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
512           enddo
513        enddo
514     enddo
515!$OMP ENDDO NOWAIT
516!$OMP BARRIER
517  enddo
518
519  !    tendances
520  !    ---------
521
522  !   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
531  !   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
542
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'
583  !=====================================================================
584  !   PAS DE TEMPS D'ECRITURE
585  !=====================================================================
586  if (icum.eq.ncum) then
587  !=====================================================================
588
589  IF (prt_level > 5) &
590        WRITE(lunout,*)'Pas d ecriture'
591
592  jjb=jj_begin
593  jje=jj_end
594
595  !   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) &
600            /masse_cum(:,jjb:jje,l)
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
627!$OMP BARRIER
628
629  jjb=jj_begin
630  jje=jj_end
631
632
633  !   A retravailler eventuellement
634  !   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
642
643  !=====================================================================
644  !   Transport m�ridien
645  !=====================================================================
646
647  !   cumul zonal des masses des mailles
648  !   ----------------------------------
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
659!$OMP BARRIER
660
661  call Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
662  do iQ=1,nQ
663    call Register_Hallo_u(Q_cum(1,jjb_u,1,iQ),llm,0,1,1,0,Req)
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
687!$OMP BARRIER
688
689  ! print*,'3OK'
690  !   --------------------------------------------------------------
691  !   calcul de la moyenne zonale du transport :
692  !   ------------------------------------------
693  !
694  !                                 --
695  ! TOT : la circulation totale       [ vq ]
696  !
697  !                                  -     -
698  ! MMC : mean meridional circulation [ v ] [ q ]
699  !
700  !                                 ----      --       - -
701  ! TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
702  !
703  !                                 - * - *       - -       -     -
704  ! STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
705  !
706  !                                          - -
707  !    on utilise aussi l'intermediaire TMP :  [ v q ]
708  !
709  !    la variable zfactv transforme un transport meridien cumule
710  !    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
711  !
712  !   --------------------------------------------------------------
713
714
715  !   ----------------------------------------
716  !   Transport dans le plan latitude-altitude
717  !   ----------------------------------------
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
730           ! print*,'j,l,iQ=',j,l,iQ
731  !   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                    +flux_vQ_cum(i,j,l,iQ)
735              zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+ &
736                    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                    /(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
741           ! print*,'aOK'
742  !   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
752  !   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
764  !   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
777  ! print*,'4OK'
778  !   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              zvQ(jjb:jje,:,itr,iQ) &
790              ,jjn*llm,ndex3d)
791     enddo
792     call histwrite(fileid,'psi'//nom(iQ), &
793           itau,psiQ(jjb:jje,1:llm,iQ) &
794           ,jjn*llm,ndex3d)
795  enddo
796
797  call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm) &
798        ,jjn*llm,ndex3d)
799  call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm) &
800        ,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        jjn*llm,ndex3d)
804
805  endif
806
807
808  !   -----------------
809  !   Moyenne verticale
810  !   -----------------
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                 +zvQ(jjb:jje,l,itr,iQ) &
823                 *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              zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
828     enddo
829  enddo
830!$OMP END MASTER
831  ! on doit pouvoir tracer systematiquement la fonction de courant.
832
833  !=====================================================================
834  !/////////////////////////////////////////////////////////////////////
835  icum=0                  !///////////////////////////////////////
836  endif ! icum.eq.ncum    !///////////////////////////////////////
837  !/////////////////////////////////////////////////////////////////////
838  !=====================================================================
839!$OMP MASTER
840  call histsync(fileid)
841!$OMP END MASTER
842
843
844  return
845end subroutine bilan_dyn_loc
Note: See TracBrowser for help on using the repository browser.