Context Navigation

bilan_dyn_p.F @ 5394

Last change on this file since 5394 was 1885, checked in by Ehouarn Millour, 11 years ago
Make dyn3dpar/bilan_dyn_p.F work in MPI and OpenMP. EM
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.1 KB

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

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

Download in other formats:

Original Format