Context Navigation

bilan_dyn_loc.F @ 3502

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

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

Download in other formats:

Original Format