Context Navigation

bilan_dyn_loc.F @ 1823

Last change on this file since 1823 was 1823, checked in by Ehouarn Millour, 11 years ago
Remplacement de parallel.F90 (en conflit avec orchidée) par parallel_lmdz.F90. UG ......................................... Renaming parallel.F90 (conflicting with orchidée) into parallel_lmdz.F90. UG
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_lmdz
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