Context Navigation

bilan_dyn_p.F @ 3452

Last change on this file since 3452 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 22.2 KB

Line
1	!
2	! $Id: bilan_dyn_p.F 1907 2013-11-26 13:10:46Z lguez $
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 comvert_mod, ONLY: presnivs
20	USE comconst_mod, ONLY: cpp,pi
21	USE temps_mod, ONLY: annee_ref,day_ref,itau_dyn
22	IMPLICIT NONE
23
24	#include "dimensions.h"
25	#include "paramet.h"
26	#include "comgeom2.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
811

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

Download in other formats:

Original Format