Context Navigation

bilan_dyn_loc.F @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 11 years ago

Version testing basée sur la r1706

Testing release based on r1706

File size: 23.6 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
7	c AFAIRE
8	c Prevoir en champ nq+1 le diagnostique de l'energie
9	c en faisant Qzon=Cv T + L * ...
10	c vQ..A=Cp T + L * ...
11
12	#ifdef CPP_IOIPSL
13	USE IOIPSL
14	#endif
15	USE parallel
16	USE mod_hallo
17	use misc_mod
18	use write_field
19	IMPLICIT NONE
20
21	#include "dimensions.h"
22	#include "paramet.h"
23	#include "comconst.h"
24	#include "comvert.h"
25	#include "comgeom2.h"
26	#include "temps.h"
27	#include "iniprint.h"
28
29	c====================================================================
30	c
31	c Sous-programme consacre à des diagnostics dynamiques de base
32	c
33	c
34	c De facon generale, les moyennes des scalaires Q sont ponderees par
35	c la masse.
36	c
37	c Les flux de masse sont eux simplement moyennes.
38	c
39	c====================================================================
40
41	c Arguments :
42	c ===========
43
44	integer ntrac
45	real dt_app,dt_cum
46	real ps(iip1,jjb_u:jje_u)
47	real masse(iip1,jjb_u:jje_u,llm),pk(iip1,jjb_u:jje_u,llm)
48	real flux_u(iip1,jjb_u:jje_u,llm)
49	real flux_v(iip1,jjb_v:jje_v,llm)
50	real teta(iip1,jjb_u:jje_u,llm)
51	real phi(iip1,jjb_u:jje_u,llm)
52	real ucov(iip1,jjb_u:jje_u,llm)
53	real vcov(iip1,jjb_v:jje_v,llm)
54	real trac(iip1,jjb_u:jje_u,llm,ntrac)
55
56	c Local :
57	c =======
58
59	integer,SAVE :: icum,ncum
60	!$OMP THREADPRIVATE(icum,ncum)
61	LOGICAL,SAVE :: first=.TRUE.
62	!$OMP THREADPRIVATE(first)
63
64	real zz,zqy
65	REAl,SAVE,ALLOCATABLE :: zfactv(:,:)
66
67	INTEGER,PARAMETER :: nQ=7
68
69
70	cym character*6 nom(nQ)
71	cym character*6 unites(nQ)
72	character(len=6),save :: nom(nQ)
73	character(len=6),save :: unites(nQ)
74
75	character(len=10) file
76	integer ifile
77	parameter (ifile=4)
78
79	integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
80	INTEGER,PARAMETER :: iovap=6,iun=7
81	integer,PARAMETER :: i_sortie=1
82
83	real,SAVE :: time=0.
84	integer,SAVE :: itau=0.
85	!$OMP THREADPRIVATE(time,itau)
86
87	real ww
88
89	c variables dynamiques intermédiaires
90	REAL,SAVE,ALLOCATABLE :: vcont(:,:,:),ucont(:,:,:)
91	REAL,SAVE,ALLOCATABLE :: ang(:,:,:),unat(:,:,:)
92	REAL,SAVE,ALLOCATABLE :: massebx(:,:,:),masseby(:,:,:)
93	REAL,SAVE,ALLOCATABLE :: vorpot(:,:,:)
94	REAL,SAVE,ALLOCATABLE :: w(:,:,:),ecin(:,:,:),convm(:,:,:)
95	REAL,SAVE,ALLOCATABLE :: bern(:,:,:)
96
97	c champ contenant les scalaires advectés.
98	real,SAVE,ALLOCATABLE :: Q(:,:,:,:)
99
100	c champs cumulés
101	real,SAVE,ALLOCATABLE :: ps_cum(:,:)
102	real,SAVE,ALLOCATABLE :: masse_cum(:,:,:)
103	real,SAVE,ALLOCATABLE :: flux_u_cum(:,:,:)
104	real,SAVE,ALLOCATABLE :: flux_v_cum(:,:,:)
105	real,SAVE,ALLOCATABLE :: Q_cum(:,:,:,:)
106	real,SAVE,ALLOCATABLE :: flux_uQ_cum(:,:,:,:)
107	real,SAVE,ALLOCATABLE :: flux_vQ_cum(:,:,:,:)
108	real,SAVE,ALLOCATABLE :: flux_wQ_cum(:,:,:,:)
109	real,SAVE,ALLOCATABLE :: dQ(:,:,:,:)
110
111
112	c champs de tansport en moyenne zonale
113	integer ntr,itr
114	parameter (ntr=5)
115
116	cym character*10 znom(ntr,nQ)
117	cym character*20 znoml(ntr,nQ)
118	cym character*10 zunites(ntr,nQ)
119	character*10,save :: znom(ntr,nQ)
120	character*20,save :: znoml(ntr,nQ)
121	character*10,save :: zunites(ntr,nQ)
122
123	INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5
124
125	character*3 ctrs(ntr)
126	data ctrs/' ','TOT','MMC','TRS','STN'/
127
128	real,SAVE,ALLOCATABLE :: zvQ(:,:,:,:),zvQtmp(:,:)
129	real,SAVE,ALLOCATABLE :: zavQ(:,:,:),psiQ(:,:,:)
130	real,SAVE,ALLOCATABLE :: zmasse(:,:),zamasse(:)
131
132	real,SAVE,ALLOCATABLE :: zv(:,:),psi(:,:)
133
134	integer i,j,l,iQ
135
136
137	c Initialisation du fichier contenant les moyennes zonales.
138	c ---------------------------------------------------------
139
140	character*10 infile
141
142	integer fileid
143	integer thoriid, zvertiid
144	save fileid
145
146	INTEGER,SAVE,ALLOCATABLE :: ndex3d(:)
147
148	C Variables locales
149	C
150	integer tau0
151	real zjulian
152	character*3 str
153	character*10 ctrac
154	integer ii,jj
155	integer zan, dayref
156	C
157	real,SAVE,ALLOCATABLE :: rlong(:),rlatg(:)
158	integer :: jjb,jje,jjn,ijb,ije
159	type(Request) :: Req
160
161	! definition du domaine d'ecriture pour le rebuild
162
163	INTEGER,DIMENSION(1) :: ddid
164	INTEGER,DIMENSION(1) :: dsg
165	INTEGER,DIMENSION(1) :: dsl
166	INTEGER,DIMENSION(1) :: dpf
167	INTEGER,DIMENSION(1) :: dpl
168	INTEGER,DIMENSION(1) :: dhs
169	INTEGER,DIMENSION(1) :: dhe
170
171	INTEGER :: bilan_dyn_domain_id
172
173
174	c=====================================================================
175	c Initialisation
176	c=====================================================================
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(jjb_v:jje_v))
218	ALLOCATE(rlatg(jjb_v:jje_v))
219
220	!$OMP END MASTER
221	!$OMP BARRIER
222	icum=0
223	c initialisation des fichiers
224	first=.false.
225	c ncum est la frequence de stokage en pas de temps
226	ncum=dt_cum/dt_app
227	if (abs(ncumdt_app-dt_cum).gt.1.e-5dt_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	stop
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	c Initialisation du fichier contenant les moyennes zonales.
254	c ---------------------------------------------------------
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(jjb:jje), jjn, rlatg(jjb:jje),
287	. 1, 1, 1, jjn,
288	. tau0, zjulian, dt_cum, thoriid, fileid,
289	. bilan_dyn_domain_id)
290
291	C
292	C Appel a histvert pour la grille verticale
293	C
294	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
295	. llm, presnivs, zvertiid)
296	C
297	C 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	c Declarations des champs avec dimension verticale
313	c 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	c Declarations pour les fonctions de courant
324	c 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	c Declarations pour les champs de transport d'air
333	c 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	c Declarations pour les fonctions de courant
341	c 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	c Declaration des champs 1D de transport en latitude
348	c 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	c print*,'8HISTDEF'
359	CALL histend(fileid)
360
361	!$OMP END MASTER
362	endif
363
364
365	c=====================================================================
366	c Calcul des champs dynamiques
367	c ----------------------------
368
369	jjb=jj_begin
370	jje=jj_end
371
372	c é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	c$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	c 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	c=====================================================================
405	c Cumul
406	c=====================================================================
407	c
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	c 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,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
446	flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
447	. +flux_u(:,jjb:jje,:)
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,:)=flux_v_cum(:,jjb:jje,:)
456	. +flux_v(:,jjb:jje,:)
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,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
467	. +Q(:,jjb:jje,:,iQ)*masse(:,jjb:jje,:)
468	ENDDO
469	!$OMP END DO NOWAIT
470	enddo
471
472	c=====================================================================
473	c FLUX ET TENDANCES
474	c=====================================================================
475
476	c Flux longitudinal
477	c -----------------
478	do iQ=1,nQ
479	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
480	do l=1,llm
481	do j=jjb,jje
482	do i=1,iim
483	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
484	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
485	enddo
486	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
487	enddo
488	enddo
489	!$OMP END DO NOWAIT
490	enddo
491
492	c flux méridien
493	c -------------
494	do iQ=1,nQ
495	call Register_Hallo_u(Q(1,1,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	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
512	enddo
513	enddo
514	enddo
515	!$OMP END DO NOWAIT
516	enddo
517
518
519	c tendances
520	c ---------
521
522	c 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	c 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	jjb=jj_begin
543	jje=jj_end
544
545	! do iQ=1,nQ
546	! do l=1,llm-1
547	! do j=jjb,jje
548	! do i=1,iip1
549	! ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
550	! dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
551	! dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
552	! enddo
553	! enddo
554	! enddo
555	! enddo
556
557	do iQ=1,nQ
558	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
559	do l=1,llm
560	IF (l<llm) THEN
561	do j=jjb,jje
562	do i=1,iip1
563	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
564	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
565	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
566	enddo
567	enddo
568	ENDIF
569	IF (l>2) THEN
570	do j=jjb,jje
571	do i=1,iip1
572	ww=-0.5w(i,j,l)(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
573	dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
574	enddo
575	enddo
576	ENDIF
577	enddo
578	!$OMP ENDDO NOWAIT
579	enddo
580	IF (prt_level > 5)
581	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
582	c=====================================================================
583	c PAS DE TEMPS D'ECRITURE
584	c=====================================================================
585	if (icum.eq.ncum) then
586	c=====================================================================
587
588	IF (prt_level > 5)
589	. WRITE(lunout,*)'Pas d ecriture'
590
591	jjb=jj_begin
592	jje=jj_end
593
594	c Normalisation
595	do iQ=1,nQ
596	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
597	do l=1,llm
598	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
599	. /masse_cum(:,jjb:jje,l)
600	enddo
601	!$OMP ENDDO NOWAIT
602	enddo
603
604	zz=1./REAL(ncum)
605
606	!$OMP MASTER
607	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
608	!$OMP END MASTER
609
610	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
611	DO l=1,llm
612	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
613	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
614	flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
615	dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
616	ENDDO
617	!$OMP ENDDO NOWAIT
618
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
628	jjb=jj_begin
629	jje=jj_end
630
631
632	c A retravailler eventuellement
633	c division de dQ par la masse pour revenir aux bonnes grandeurs
634	do iQ=1,nQ
635	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
636	DO l=1,llm
637	dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
638	ENDDO
639	!$OMP ENDDO NOWAIT
640	enddo
641
642	c=====================================================================
643	c Transport méridien
644	c=====================================================================
645
646	c cumul zonal des masses des mailles
647	c ----------------------------------
648	jjb=jj_begin
649	jje=jj_end
650	if (pole_sud) jje=jj_end-1
651
652	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
653	DO l=1,llm
654	zv(jjb:jje,l)=0.
655	zmasse(jjb:jje,l)=0.
656	ENDDO
657	!$OMP ENDDO NOWAIT
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,1,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
686	c print*,'3OK'
687	c --------------------------------------------------------------
688	c calcul de la moyenne zonale du transport :
689	c ------------------------------------------
690	c
691	c --
692	c TOT : la circulation totale [ vq ]
693	c
694	c - -
695	c MMC : mean meridional circulation [ v ] [ q ]
696	c
697	c ---- -- - -
698	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
699	c
700	c - * - * - - - -
701	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
702	c
703	c - -
704	c on utilise aussi l'intermediaire TMP : [ v q ]
705	c
706	c la variable zfactv transforme un transport meridien cumule
707	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
708	c
709	c --------------------------------------------------------------
710
711
712	c ----------------------------------------
713	c Transport dans le plan latitude-altitude
714	c ----------------------------------------
715
716	jjb=jj_begin
717	jje=jj_end
718	if (pole_sud) jje=jj_end-1
719
720	zvQ=0.
721	psiQ=0.
722	do iQ=1,nQ
723	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
724	do l=1,llm
725	zvQtmp(:,l)=0.
726	do j=jjb,jje
727	c print*,'j,l,iQ=',j,l,iQ
728	c Calcul des moyennes zonales du transort total et de zvQtmp
729	do i=1,iim
730	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
731	s +flux_vQ_cum(i,j,l,iQ)
732	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
733	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
734	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
735	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
736	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
737	enddo
738	c print*,'aOK'
739	c Decomposition
740	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
741	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
742	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
743	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
744	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
745	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
746	enddo
747	enddo
748	!$OMP ENDDO NOWAIT
749	c fonction de courant meridienne pour la quantite Q
750	!$OMP BARRIER
751	!$OMP MASTER
752	do l=llm,1,-1
753	do j=jjb,jje
754	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
755	enddo
756	enddo
757	!$OMP END MASTER
758	!$OMP BARRIER
759	enddo
760
761	c fonction de courant pour la circulation meridienne moyenne
762	!$OMP BARRIER
763	!$OMP MASTER
764	psi(jjb:jje,:)=0.
765	do l=llm,1,-1
766	do j=jjb,jje
767	psi(j,l)=psi(j,l+1)+zv(j,l)
768	zv(j,l)=zv(j,l)*zfactv(j,l)
769	enddo
770	enddo
771	!$OMP END MASTER
772	!$OMP BARRIER
773
774	c print*,'4OK'
775	c sorties proprement dites
776	!$OMP MASTER
777	if (i_sortie.eq.1) then
778	jjb=jj_begin
779	jje=jj_end
780	jjn=jj_nb
781	if (pole_sud) jje=jj_end-1
782	if (pole_sud) jjn=jj_nb-1
783	do iQ=1,nQ
784	do itr=1,ntr
785	call histwrite(fileid,znom(itr,iQ),itau,
786	s zvQ(jjb:jje,:,itr,iQ)
787	s ,jjn*llm,ndex3d)
788	enddo
789	call histwrite(fileid,'psi'//nom(iQ),
790	s itau,psiQ(jjb:jje,1:llm,iQ)
791	s ,jjn*llm,ndex3d)
792	enddo
793
794	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
795	s ,jjn*llm,ndex3d)
796	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
797	s ,jjn*llm,ndex3d)
798	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
799	call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
800	s jjn*llm,ndex3d)
801
802	endif
803
804
805	c -----------------
806	c Moyenne verticale
807	c -----------------
808
809	zamasse(jjb:jje)=0.
810	do l=1,llm
811	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
812	enddo
813
814	zavQ(jjb:jje,:,:)=0.
815	do iQ=1,nQ
816	do itr=2,ntr
817	do l=1,llm
818	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
819	s +zvQ(jjb:jje,l,itr,iQ)
820	s *zmasse(jjb:jje,l)
821	enddo
822	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
823	call histwrite(fileid,'a'//znom(itr,iQ),itau,
824	s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
825	enddo
826	enddo
827	!$OMP END MASTER
828	c on doit pouvoir tracer systematiquement la fonction de courant.
829
830	c=====================================================================
831	c/////////////////////////////////////////////////////////////////////
832	icum=0 !///////////////////////////////////////
833	endif ! icum.eq.ncum !///////////////////////////////////////
834	c/////////////////////////////////////////////////////////////////////
835	c=====================================================================
836
837	return
838	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format