Context Navigation

bilan_dyn_loc.F @ 1632

Last change on this file since 1632 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

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
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	flux_v_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,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
447	flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
448	. +flux_u(:,jjb:jje,:)
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,:)=flux_v_cum(:,jjb:jje,:)
457	. +flux_v(:,jjb:jje,:)
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,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
468	. +Q(:,jjb:jje,:,iQ)*masse(:,jjb:jje,:)
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,1,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 END DO NOWAIT
517	enddo
518
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	jjb=jj_begin
544	jje=jj_end
545
546	! do iQ=1,nQ
547	! do l=1,llm-1
548	! do j=jjb,jje
549	! do i=1,iip1
550	! ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
551	! dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
552	! dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
553	! enddo
554	! enddo
555	! enddo
556	! enddo
557
558	do iQ=1,nQ
559	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
560	do l=1,llm
561	IF (l<llm) THEN
562	do j=jjb,jje
563	do i=1,iip1
564	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
565	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
566	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
567	enddo
568	enddo
569	ENDIF
570	IF (l>2) THEN
571	do j=jjb,jje
572	do i=1,iip1
573	ww=-0.5w(i,j,l)(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
574	dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
575	enddo
576	enddo
577	ENDIF
578	enddo
579	!$OMP ENDDO NOWAIT
580	enddo
581	IF (prt_level > 5)
582	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
583	c=====================================================================
584	c PAS DE TEMPS D'ECRITURE
585	c=====================================================================
586	if (icum.eq.ncum) then
587	c=====================================================================
588
589	IF (prt_level > 5)
590	. WRITE(lunout,*)'Pas d ecriture'
591
592	jjb=jj_begin
593	jje=jj_end
594
595	c Normalisation
596	do iQ=1,nQ
597	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
598	do l=1,llm
599	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
600	. /masse_cum(:,jjb:jje,l)
601	enddo
602	!$OMP ENDDO NOWAIT
603	enddo
604
605	zz=1./REAL(ncum)
606
607	!$OMP MASTER
608	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
609	!$OMP END MASTER
610
611	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
612	DO l=1,llm
613	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
614	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
615	flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
616	dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
617	ENDDO
618	!$OMP ENDDO NOWAIT
619
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
629	jjb=jj_begin
630	jje=jj_end
631
632
633	c A retravailler eventuellement
634	c division de dQ par la masse pour revenir aux bonnes grandeurs
635	do iQ=1,nQ
636	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
637	DO l=1,llm
638	dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
639	ENDDO
640	!$OMP ENDDO NOWAIT
641	enddo
642
643	c=====================================================================
644	c Transport méridien
645	c=====================================================================
646
647	c cumul zonal des masses des mailles
648	c ----------------------------------
649	jjb=jj_begin
650	jje=jj_end
651	if (pole_sud) jje=jj_end-1
652
653	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
654	DO l=1,llm
655	zv(jjb:jje,l)=0.
656	zmasse(jjb:jje,l)=0.
657	ENDDO
658	!$OMP ENDDO NOWAIT
659
660	call Register_Hallo_u(masse_cum,llm,1,1,1,1,Req)
661	do iQ=1,nQ
662	call Register_Hallo_u(Q_cum(1,1,1,iQ),llm,0,1,1,0,Req)
663	enddo
664
665	call SendRequest(Req)
666	!$OMP BARRIER
667	call WaitRequest(Req)
668
669	call massbar_loc(masse_cum,massebx,masseby)
670
671	jjb=jj_begin
672	jje=jj_end
673	if (pole_sud) jje=jj_end-1
674
675	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
676	do l=1,llm
677	do j=jjb,jje
678	do i=1,iim
679	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
680	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
681	enddo
682	zfactv(j,l)=cv(1,j)/zmasse(j,l)
683	enddo
684	enddo
685	!$OMP ENDDO NOWAIT
686
687	c print*,'3OK'
688	c --------------------------------------------------------------
689	c calcul de la moyenne zonale du transport :
690	c ------------------------------------------
691	c
692	c --
693	c TOT : la circulation totale [ vq ]
694	c
695	c - -
696	c MMC : mean meridional circulation [ v ] [ q ]
697	c
698	c ---- -- - -
699	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
700	c
701	c - * - * - - - -
702	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
703	c
704	c - -
705	c on utilise aussi l'intermediaire TMP : [ v q ]
706	c
707	c la variable zfactv transforme un transport meridien cumule
708	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
709	c
710	c --------------------------------------------------------------
711
712
713	c ----------------------------------------
714	c Transport dans le plan latitude-altitude
715	c ----------------------------------------
716
717	jjb=jj_begin
718	jje=jj_end
719	if (pole_sud) jje=jj_end-1
720
721	zvQ=0.
722	psiQ=0.
723	do iQ=1,nQ
724	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
725	do l=1,llm
726	zvQtmp(:,l)=0.
727	do j=jjb,jje
728	c print*,'j,l,iQ=',j,l,iQ
729	c Calcul des moyennes zonales du transort total et de zvQtmp
730	do i=1,iim
731	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
732	s +flux_vQ_cum(i,j,l,iQ)
733	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
734	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
735	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
736	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
737	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
738	enddo
739	c print*,'aOK'
740	c Decomposition
741	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
742	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
743	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
744	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
745	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
746	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
747	enddo
748	enddo
749	!$OMP ENDDO NOWAIT
750	c fonction de courant meridienne pour la quantite Q
751	!$OMP BARRIER
752	!$OMP MASTER
753	do l=llm,1,-1
754	do j=jjb,jje
755	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
756	enddo
757	enddo
758	!$OMP END MASTER
759	!$OMP BARRIER
760	enddo
761
762	c fonction de courant pour la circulation meridienne moyenne
763	!$OMP BARRIER
764	!$OMP MASTER
765	psi(jjb:jje,:)=0.
766	do l=llm,1,-1
767	do j=jjb,jje
768	psi(j,l)=psi(j,l+1)+zv(j,l)
769	zv(j,l)=zv(j,l)*zfactv(j,l)
770	enddo
771	enddo
772	!$OMP END MASTER
773	!$OMP BARRIER
774
775	c print*,'4OK'
776	c sorties proprement dites
777	!$OMP MASTER
778	if (i_sortie.eq.1) then
779	jjb=jj_begin
780	jje=jj_end
781	jjn=jj_nb
782	if (pole_sud) jje=jj_end-1
783	if (pole_sud) jjn=jj_nb-1
784	do iQ=1,nQ
785	do itr=1,ntr
786	call histwrite(fileid,znom(itr,iQ),itau,
787	s zvQ(jjb:jje,:,itr,iQ)
788	s ,jjn*llm,ndex3d)
789	enddo
790	call histwrite(fileid,'psi'//nom(iQ),
791	s itau,psiQ(jjb:jje,1:llm,iQ)
792	s ,jjn*llm,ndex3d)
793	enddo
794
795	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
796	s ,jjn*llm,ndex3d)
797	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
798	s ,jjn*llm,ndex3d)
799	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
800	call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
801	s jjn*llm,ndex3d)
802
803	endif
804
805
806	c -----------------
807	c Moyenne verticale
808	c -----------------
809
810	zamasse(jjb:jje)=0.
811	do l=1,llm
812	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
813	enddo
814
815	zavQ(jjb:jje,:,:)=0.
816	do iQ=1,nQ
817	do itr=2,ntr
818	do l=1,llm
819	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
820	s +zvQ(jjb:jje,l,itr,iQ)
821	s *zmasse(jjb:jje,l)
822	enddo
823	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
824	call histwrite(fileid,'a'//znom(itr,iQ),itau,
825	s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
826	enddo
827	enddo
828	!$OMP END MASTER
829	c on doit pouvoir tracer systematiquement la fonction de courant.
830
831	c=====================================================================
832	c/////////////////////////////////////////////////////////////////////
833	icum=0 !///////////////////////////////////////
834	endif ! icum.eq.ncum !///////////////////////////////////////
835	c/////////////////////////////////////////////////////////////////////
836	c=====================================================================
837
838	return
839	end

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

Download in other formats:

Original Format