Context Navigation

bilan_dyn_p.F

Last change on this file was 3990, checked in by millour, 8 years ago
An updated LMDZ5 (based on LMDZ rev 2786) to play with. EM
File size: 22.3 KB

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

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

Download in other formats:

Original Format