Context Navigation

source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn/bilan_dyn_p.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago

File size: 30.6 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
13
14
15	USE parallel_lmdz
16	USE mod_hallo
17	use misc_mod
18	use write_field_p
19	IMPLICIT NONE
20
21	!-----------------------------------------------------------------------
22	! INCLUDE 'dimensions.h'
23	!
24	! dimensions.h contient les dimensions du modele
25	! ndm est tel que iim=2**ndm
26	!-----------------------------------------------------------------------
27
28	INTEGER iim,jjm,llm,ndm
29
30	PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
31
32	!-----------------------------------------------------------------------
33	!
34	! $Header$
35	!
36	!
37	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
38	! veillez n'utiliser que des ! pour les commentaires
39	! et bien positionner les & des lignes de continuation
40	! (les placer en colonne 6 et en colonne 73)
41	!
42	!
43	!-----------------------------------------------------------------------
44	! INCLUDE 'paramet.h'
45
46	INTEGER iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
47	INTEGER kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
48	INTEGER ijmllm,mvar
49	INTEGER jcfil,jcfllm
50
51	PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3 &
52	& ,jjp1=jjm+1-1/jjm)
53	PARAMETER( llmp1 = llm+1, llmp2 = llm+2, llmm1 = llm-1 )
54	PARAMETER( kftd = iim/2 -ndm )
55	PARAMETER( ip1jm = iip1jjm, ip1jmp1= iip1jjp1 )
56	PARAMETER( ip1jmi1= ip1jm - iip1 )
57	PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
58	PARAMETER( mvar= ip1jmp1( 2llm+1) + ijmllm )
59	PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )
60
61	!-----------------------------------------------------------------------
62	!
63	! $Id: comconst.h 1437 2010-09-30 08:29:10Z emillour $
64	!
65	!-----------------------------------------------------------------------
66	! INCLUDE comconst.h
67
68	COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl, &
69	& iflag_top_bound,mode_top_bound
70	COMMON/comconstr/dtvr,daysec, &
71	& pi,dtphys,dtdiss,rad,r,kappa,cotot,unsim,g,omeg &
72	& ,dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta &
73	& ,dissip_pupstart ,tau_top_bound, &
74	& daylen,molmass, ihf
75	COMMON/cpdetvenus/cpp,nu_venus,t0_venus
76
77	INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
78	REAL dtvr ! dynamical time step (in s)
79	REAL daysec !length (in s) of a standard day
80	REAL pi ! something like 3.14159....
81	REAL dtphys ! (s) time step for the physics
82	REAL dtdiss ! (s) time step for the dissipation
83	REAL rad ! (m) radius of the planet
84	REAL r ! Reduced Gas constant r=R/mu
85	! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
86	REAL cpp ! Cp
87	REAL kappa ! kappa=R/Cp
88	REAL cotot
89	REAL unsim ! = 1./iim
90	REAL g ! (m/s2) gravity
91	REAL omeg ! (rad/s) rotation rate of the planet
92	! Dissipation factors, for Earth model:
93	REAL dissip_factz,dissip_zref !dissip_deltaz
94	! Dissipation factors, for other planets:
95	REAL dissip_fac_mid,dissip_fac_up,dissip_deltaz,dissip_hdelta
96	REAL dissip_pupstart
97	INTEGER iflag_top_bound,mode_top_bound
98	REAL tau_top_bound
99	REAL daylen ! length of solar day, in 'standard' day length
100	REAL molmass ! (g/mol) molar mass of the atmosphere
101
102	REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
103	REAL ihf ! (W/m2) intrinsic heat flux for giant planets
104
105
106	!-----------------------------------------------------------------------
107	!
108	! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
109	!
110	!-----------------------------------------------------------------------
111	! INCLUDE 'comvert.h'
112
113	COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm), &
114	& pa,preff,nivsigs(llm),nivsig(llm+1), &
115	& aps(llm),bps(llm),scaleheight,pseudoalt(llm)
116
117	common/comverti/disvert_type, pressure_exner
118
119	real ap ! hybrid pressure contribution at interlayers
120	real bp ! hybrid sigma contribution at interlayer
121	real presnivs ! (reference) pressure at mid-layers
122	real dpres
123	real pa ! reference pressure (Pa) at which hybrid coordinates
124	! become purely pressure
125	real preff ! reference surface pressure (Pa)
126	real nivsigs
127	real nivsig
128	real aps ! hybrid pressure contribution at mid-layers
129	real bps ! hybrid sigma contribution at mid-layers
130	real scaleheight ! atmospheric (reference) scale height (km)
131	real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
132	! preff and scaleheight
133
134	integer disvert_type ! type of vertical discretization:
135	! 1: Earth (default for planet_type==earth),
136	! automatic generation
137	! 2: Planets (default for planet_type!=earth),
138	! using 'z2sig.def' (or 'esasig.def) file
139
140	logical pressure_exner
141	! compute pressure inside layers using Exner function, else use mean
142	! of pressure values at interfaces
143
144	!-----------------------------------------------------------------------
145	!
146	! $Header$
147	!
148	!CDK comgeom2
149	COMMON/comgeom/ &
150	& cu(iip1,jjp1),cv(iip1,jjm),unscu2(iip1,jjp1),unscv2(iip1,jjm) , &
151	& aire(iip1,jjp1),airesurg(iip1,jjp1),aireu(iip1,jjp1) , &
152	& airev(iip1,jjm),unsaire(iip1,jjp1),apoln,apols , &
153	& unsairez(iip1,jjm),airuscv2(iip1,jjm),airvscu2(iip1,jjm) , &
154	& aireij1(iip1,jjp1),aireij2(iip1,jjp1),aireij3(iip1,jjp1) , &
155	& aireij4(iip1,jjp1),alpha1(iip1,jjp1),alpha2(iip1,jjp1) , &
156	& alpha3(iip1,jjp1),alpha4(iip1,jjp1),alpha1p2(iip1,jjp1) , &
157	& alpha1p4(iip1,jjp1),alpha2p3(iip1,jjp1),alpha3p4(iip1,jjp1) , &
158	& fext(iip1,jjm),constang(iip1,jjp1), rlatu(jjp1),rlatv(jjm), &
159	& rlonu(iip1),rlonv(iip1),cuvsurcv(iip1,jjm),cvsurcuv(iip1,jjm) , &
160	& cvusurcu(iip1,jjp1),cusurcvu(iip1,jjp1) , &
161	& cuvscvgam1(iip1,jjm),cuvscvgam2(iip1,jjm),cvuscugam1(iip1,jjp1), &
162	& cvuscugam2(iip1,jjp1),cvscuvgam(iip1,jjm),cuscvugam(iip1,jjp1) , &
163	& unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2 , &
164	& unsair_gam1(iip1,jjp1),unsair_gam2(iip1,jjp1) , &
165	& unsairz_gam(iip1,jjm),aivscu2gam(iip1,jjm),aiuscv2gam(iip1,jjm) &
166	& , xprimu(iip1),xprimv(iip1)
167
168
169	REAL &
170	& cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,apoln,apols,unsaire &
171	& ,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4 , &
172	& alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 , &
173	& fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2 , &
174	& cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1 , &
175	& unsapolnga2,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2 , &
176	& unsairz_gam,aivscu2gam,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu , &
177	& cusurcvu,xprimu,xprimv
178	!
179	! $Id: temps.h 1577 2011-10-20 15:06:47Z fairhead $
180	!
181	! ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
182	! veillez n'utiliser que des ! pour les commentaires
183	! et bien positionner les & des lignes de continuation
184	! (les placer en colonne 6 et en colonne 73)
185	!
186	!
187	! jD_ref = jour julien de la date de reference (lancement de l'experience)
188	! hD_ref = "heure" julienne de la date de reference
189	!-----------------------------------------------------------------------
190	! INCLUDE 'temps.h'
191
192	COMMON/temps_r/dt,jD_ref,jH_ref,start_time,hour_ini
193	COMMON/temps_i/day_ini,day_end,annee_ref,day_ref, &
194	& itau_dyn,itau_phy,itaufin
195	COMMON/temps_c/calend
196
197
198	INTEGER itaufin ! total number of dynamical steps for the run
199	INTEGER itau_dyn, itau_phy
200	INTEGER day_ini ! initial day # of simulation sequence
201	INTEGER day_end ! final day # ; i.e. day # when this simulation ends
202	INTEGER annee_ref
203	INTEGER day_ref
204	REAL dt ! (dynamics) time step (changes if doing Matsuno or LF step)
205	REAL jD_ref, jH_ref, start_time
206	CHARACTER (len=10) :: calend
207
208	! Additionnal Mars stuff:
209	real hour_ini ! initial fraction of day of simulation sequence (0=<hour_ini<1)
210
211	!-----------------------------------------------------------------------
212	!
213	! $Header$
214	!
215	!
216	! gestion des impressions de sorties et de débogage
217	! lunout: unité du fichier dans lequel se font les sorties
218	! (par defaut 6, la sortie standard)
219	! prt_level: niveau d'impression souhaité (0 = minimum)
220	!
221	INTEGER lunout, prt_level
222	COMMON /comprint/ lunout, prt_level
223
224	c====================================================================
225	c
226	c Sous-programme consacre à des diagnostics dynamiques de base
227	c
228	c
229	c De facon generale, les moyennes des scalaires Q sont ponderees par
230	c la masse.
231	c
232	c Les flux de masse sont eux simplement moyennes.
233	c
234	c====================================================================
235
236	c Arguments :
237	c ===========
238
239	integer ntrac
240	real dt_app,dt_cum
241	real ps(iip1,jjp1)
242	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
243	real flux_u(iip1,jjp1,llm)
244	real flux_v(iip1,jjm,llm)
245	real teta(iip1,jjp1,llm)
246	real phi(iip1,jjp1,llm)
247	real ucov(iip1,jjp1,llm)
248	real vcov(iip1,jjm,llm)
249	real trac(iip1,jjp1,llm,ntrac)
250
251	c Local :
252	c =======
253
254	integer,save :: icum,ncum
255	!$OMP THREADPRIVATE(icum,ncum)
256	logical,SAVE :: first=.true.
257	!$OMP THREADPRIVATE(first)
258
259	real zz,zqy
260	real,save :: zfactv(jjm,llm)
261
262	integer,parameter :: nQ=7
263
264
265	cym character*6 nom(nQ)
266	cym character*6 unites(nQ)
267	character(len=6),save :: nom(nQ)
268	character(len=6),save :: unites(nQ)
269
270	character(len=10) file
271	integer ifile
272	parameter (ifile=4)
273
274	integer,PARAMETER :: itemp=1,igeop=2,iecin=3,iang=4,iu=5
275	INTEGER,PARAMETER :: iovap=6,iun=7
276	integer,PARAMETER :: i_sortie=1
277
278	real,SAVE :: time=0.
279	integer,SAVE :: itau=0.
280	!$OMP THREADPRIVATE(time,itau)
281
282	real ww
283
284	c variables dynamiques intermédiaires
285	REAL,save :: vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
286	REAL,save :: ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
287	REAL,save :: massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
288	REAL,save :: vorpot(iip1,jjm,llm)
289	REAL,save :: w(iip1,jjp1,llm),ecin(iip1,jjp1,llm)
290	REAL,save ::convm(iip1,jjp1,llm)
291	REAL,save :: bern(iip1,jjp1,llm)
292
293	c champ contenant les scalaires advectés.
294	real,save :: Q(iip1,jjp1,llm,nQ)
295
296	c champs cumulés
297	real,save :: ps_cum(iip1,jjp1)
298	real,save :: masse_cum(iip1,jjp1,llm)
299	real,save :: flux_u_cum(iip1,jjp1,llm)
300	real,save :: flux_v_cum(iip1,jjm,llm)
301	real,save :: Q_cum(iip1,jjp1,llm,nQ)
302	real,save :: flux_uQ_cum(iip1,jjp1,llm,nQ)
303	real,save :: flux_vQ_cum(iip1,jjm,llm,nQ)
304	real,save :: flux_wQ_cum(iip1,jjp1,llm,nQ)
305	real,save :: dQ(iip1,jjp1,llm,nQ)
306
307
308	c champs de tansport en moyenne zonale
309	integer ntr,itr
310	parameter (ntr=5)
311
312	cym character*10 znom(ntr,nQ)
313	cym character*20 znoml(ntr,nQ)
314	cym character*10 zunites(ntr,nQ)
315	character*10,save :: znom(ntr,nQ)
316	character*20,save :: znoml(ntr,nQ)
317	character*10,save :: zunites(ntr,nQ)
318
319	INTEGER,PARAMETER :: iave=1,itot=2,immc=3,itrs=4,istn=5
320
321	character*3 ctrs(ntr)
322	data ctrs/' ','TOT','MMC','TRS','STN'/
323
324	real,save :: zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
325	real,save :: zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
326	real,save :: zmasse(jjm,llm),zamasse(jjm)
327
328	real,save :: zv(jjm,llm),psi(jjm,llm+1)
329
330	integer i,j,l,iQ
331
332
333	c Initialisation du fichier contenant les moyennes zonales.
334	c ---------------------------------------------------------
335
336	character*10 infile
337
338	integer fileid
339	integer thoriid, zvertiid
340	save fileid
341
342	integer,save :: ndex3d(jjm*llm)
343
344	C Variables locales
345	C
346	integer tau0
347	real zjulian
348	character*3 str
349	character*10 ctrac
350	integer ii,jj
351	integer zan, dayref
352	C
353	real,save :: rlong(jjm),rlatg(jjm)
354	integer :: jjb,jje,jjn,ijb,ije
355	type(Request),SAVE :: Req
356	!$OMP THREADPRIVATE(Req)
357
358	! definition du domaine d'ecriture pour le rebuild
359
360	INTEGER,DIMENSION(1) :: ddid
361	INTEGER,DIMENSION(1) :: dsg
362	INTEGER,DIMENSION(1) :: dsl
363	INTEGER,DIMENSION(1) :: dpf
364	INTEGER,DIMENSION(1) :: dpl
365	INTEGER,DIMENSION(1) :: dhs
366	INTEGER,DIMENSION(1) :: dhe
367
368	INTEGER :: bilan_dyn_domain_id
369
370
371	c=====================================================================
372	c Initialisation
373	c=====================================================================
374	if (adjust) return
375
376	time=time+dt_app
377	itau=itau+1
378
379	if (first) then
380
381	ndex3d=0
382
383	icum=0
384	c initialisation des fichiers
385	first=.false.
386	c ncum est la frequence de stokage en pas de temps
387	ncum=dt_cum/dt_app
388	if (abs(ncumdt_app-dt_cum).gt.1.e-5dt_app) then
389	WRITE(lunout,*)
390	. 'Pb : le pas de cumule doit etre multiple du pas'
391	WRITE(lunout,*)'dt_app=',dt_app
392	WRITE(lunout,*)'dt_cum=',dt_cum
393	stop
394	else
395	write(lunout,*) "bilan_dyn_p: ncum=",ncum
396	endif
397
398	! if (i_sortie.eq.1) then
399	! file='dynzon'
400	! if (mpi_rank==0) then
401	! call inigrads(ifile,1
402	! s ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
403	! s ,llm,presnivs,1.
404	! s ,dt_cum,file,'dyn_zon ')
405	! endif
406	! endif
407
408	!$OMP MASTER
409	nom(itemp)='T'
410	nom(igeop)='gz'
411	nom(iecin)='K'
412	nom(iang)='ang'
413	nom(iu)='u'
414	nom(iovap)='ovap'
415	nom(iun)='un'
416
417	unites(itemp)='K'
418	unites(igeop)='m2/s2'
419	unites(iecin)='m2/s2'
420	unites(iang)='ang'
421	unites(iu)='m/s'
422	unites(iovap)='kg/kg'
423	unites(iun)='un'
424
425
426	c Initialisation du fichier contenant les moyennes zonales.
427	c ---------------------------------------------------------
428
429	infile='dynzon'
430
431	zan = annee_ref
432	dayref = day_ref
433	CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
434	tau0 = itau_dyn
435
436	rlong=0.
437	rlatg=rlatv*180./pi
438
439	jjb=jj_begin
440	jje=jj_end
441	jjn=jj_nb
442	IF (pole_sud) THEN
443	jjn=jj_nb-1
444	jje=jj_end-1
445	ENDIF
446
447	ddid=(/ 2 /)
448	dsg=(/ jjm /)
449	dsl=(/ jjn /)
450	dpf=(/ jjb /)
451	dpl=(/ jje /)
452	dhs=(/ 0 /)
453	dhe=(/ 0 /)
454
455	call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
456	. 'box',bilan_dyn_domain_id)
457
458	call histbeg(trim(infile),
459	. 1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
460	. 1, 1, 1, jjn,
461	. tau0, zjulian, dt_cum, thoriid, fileid,
462	. bilan_dyn_domain_id)
463
464	C
465	C Appel a histvert pour la grille verticale
466	C
467	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
468	. llm, presnivs, zvertiid)
469	C
470	C Appels a histdef pour la definition des variables a sauvegarder
471	do iQ=1,nQ
472	do itr=1,ntr
473	if(itr.eq.1) then
474	znom(itr,iQ)=nom(iQ)
475	znoml(itr,iQ)=nom(iQ)
476	zunites(itr,iQ)=unites(iQ)
477	else
478	znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
479	znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
480	zunites(itr,iQ)='m/s * '//unites(iQ)
481	endif
482	enddo
483	enddo
484
485	c Declarations des champs avec dimension verticale
486	c print*,'1HISTDEF'
487	do iQ=1,nQ
488	do itr=1,ntr
489	IF (prt_level > 5)
490	. WRITE(lunout,*)'var ',itr,iQ
491	. ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
492	call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
493	. zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
494	. 32,'ave(X)',dt_cum,dt_cum)
495	enddo
496	c Declarations pour les fonctions de courant
497	c print*,'2HISTDEF'
498	call histdef(fileid,'psi'//nom(iQ)
499	. ,'stream fn. '//znoml(itot,iQ),
500	. zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
501	. 32,'ave(X)',dt_cum,dt_cum)
502	enddo
503
504
505	c Declarations pour les champs de transport d'air
506	c print*,'3HISTDEF'
507	call histdef(fileid, 'masse', 'masse',
508	. 'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
509	. 32, 'ave(X)', dt_cum, dt_cum)
510	call histdef(fileid, 'v', 'v',
511	. 'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
512	. 32, 'ave(X)', dt_cum, dt_cum)
513	c Declarations pour les fonctions de courant
514	c print*,'4HISTDEF'
515	call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
516	. 1,jjn,thoriid,llm,1,llm,zvertiid,
517	. 32,'ave(X)',dt_cum,dt_cum)
518
519
520	c Declaration des champs 1D de transport en latitude
521	c print*,'5HISTDEF'
522	do iQ=1,nQ
523	do itr=2,ntr
524	call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
525	. zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
526	. 32,'ave(X)',dt_cum,dt_cum)
527	enddo
528	enddo
529
530
531	c print*,'8HISTDEF'
532	CALL histend(fileid)
533
534	!$OMP END MASTER
535	!$OMP BARRIER
536	endif
537
538
539	c=====================================================================
540	c Calcul des champs dynamiques
541	c ----------------------------
542
543	jjb=jj_begin
544	jje=jj_end
545
546	c énergie cinétique
547	! ucont(:,jjb:jje,:)=0
548
549	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Req)
550	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Req)
551	call SendRequest(Req)
552	c$OMP BARRIER
553	call WaitRequest(Req)
554	c$OMP BARRIER
555
556	CALL covcont_p(llm,ucov,vcov,ucont,vcont)
557	CALL enercin_p(vcov,ucov,vcont,ucont,ecin)
558
559	c moment cinétique
560	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
561	do l=1,llm
562	ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
563	unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
564	enddo
565	!$OMP END DO
566
567	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
568	DO l=1,llm
569	Q(:,jjb:jje,l,itemp)=teta(:,jjb:jje,l)*pk(:,jjb:jje,l)/cpp
570	Q(:,jjb:jje,l,igeop)=phi(:,jjb:jje,l)
571	Q(:,jjb:jje,l,iecin)=ecin(:,jjb:jje,l)
572	Q(:,jjb:jje,l,iang)=ang(:,jjb:jje,l)
573	Q(:,jjb:jje,l,iu)=unat(:,jjb:jje,l)
574	Q(:,jjb:jje,l,iovap)=trac(:,jjb:jje,l,1)
575	Q(:,jjb:jje,l,iun)=1.
576	ENDDO
577	!$OMP END DO NOWAIT
578
579	c=====================================================================
580	c Cumul
581	c=====================================================================
582	c
583	if(icum.EQ.0) then
584	jjb=jj_begin
585	jje=jj_end
586
587	!$OMP MASTER
588	ps_cum(:,jjb:jje)=0.
589	!$OMP END MASTER
590	!$OMP BARRIER
591
592	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
593	DO l=1,llm
594	masse_cum(:,jjb:jje,l)=0.
595	flux_u_cum(:,jjb:jje,l)=0.
596	Q_cum(:,jjb:jje,l,:)=0.
597	flux_uQ_cum(:,jjb:jje,l,:)=0.
598	if (pole_sud) jje=jj_end-1
599	flux_v_cum(:,jjb:jje,l)=0.
600	flux_vQ_cum(:,jjb:jje,l,:)=0.
601	ENDDO
602	!$OMP END DO NOWAIT
603	endif
604
605	IF (prt_level > 5)
606	. WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
607	icum=icum+1
608
609	c accumulation des flux de masse horizontaux
610	jjb=jj_begin
611	jje=jj_end
612
613	!$OMP MASTER
614	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
615	!$OMP END MASTER
616	!$OMP BARRIER
617
618	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
619	DO l=1,llm
620	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)+masse(:,jjb:jje,l)
621	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)
622	. +flux_u(:,jjb:jje,l)
623	ENDDO
624	!$OMP END DO NOWAIT
625
626	if (pole_sud) jje=jj_end-1
627
628	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
629	DO l=1,llm
630	flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)
631	. +flux_v(:,jjb:jje,l)
632	ENDDO
633	!$OMP END DO NOWAIT
634
635	jjb=jj_begin
636	jje=jj_end
637
638	do iQ=1,nQ
639	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
640	DO l=1,llm
641	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
642	. +Q(:,jjb:jje,l,iQ)*masse(:,jjb:jje,l)
643	ENDDO
644	!$OMP END DO NOWAIT
645	enddo
646
647	c=====================================================================
648	c FLUX ET TENDANCES
649	c=====================================================================
650
651	c Flux longitudinal
652	c -----------------
653	do iQ=1,nQ
654	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
655	do l=1,llm
656	do j=jjb,jje
657	do i=1,iim
658	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
659	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
660	enddo
661	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
662	enddo
663	enddo
664	!$OMP END DO NOWAIT
665	enddo
666
667	c flux méridien
668	c -------------
669	do iQ=1,nQ
670	call Register_Hallo(Q(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
671	enddo
672	call SendRequest(Req)
673	!$OMP BARRIER
674	call WaitRequest(Req)
675	!$OMP BARRIER
676
677	jjb=jj_begin
678	jje=jj_end
679	if (pole_sud) jje=jj_end-1
680
681	do iQ=1,nQ
682	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
683	do l=1,llm
684	do j=jjb,jje
685	do i=1,iip1
686	flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
687	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
688	enddo
689	enddo
690	enddo
691	!$OMP END DO NOWAIT
692	enddo
693
694
695	c tendances
696	c ---------
697
698	c convergence horizontale
699	call Register_Hallo(flux_uQ_cum,ip1jmp1,llm,2,2,2,2,Req)
700	call Register_Hallo(flux_vQ_cum,ip1jm,llm,2,2,2,2,Req)
701	call SendRequest(Req)
702	!$OMP BARRIER
703	call WaitRequest(Req)
704	c$OMP BARRIER
705
706	call convflu_p(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
707
708	c calcul de la vitesse verticale
709	call Register_Hallo(flux_u_cum,ip1jmp1,llm,2,2,2,2,Req)
710	call Register_Hallo(flux_v_cum,ip1jm,llm,2,2,2,2,Req)
711	call SendRequest(Req)
712	!$OMP BARRIER
713	call WaitRequest(Req)
714	c$OMP BARRIER
715
716	call convmas_p(flux_u_cum,flux_v_cum,convm)
717	CALL vitvert_p(convm,w)
718	!$OMP BARRIER
719
720	jjb=jj_begin
721	jje=jj_end
722
723	do iQ=1,nQ
724	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
725	do l=1,llm
726	IF (l<llm) THEN
727	do j=jjb,jje
728	do i=1,iip1
729	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
730	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
731	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
732	enddo
733	enddo
734	ENDIF
735	IF (l>2) THEN
736	do j=jjb,jje
737	do i=1,iip1
738	ww=-0.5w(i,j,l)(Q(i,j,l-1,iQ)+Q(i,j,l,iQ))
739	dQ(i,j,l,iQ)=dQ(i,j,l,iQ)+ww
740	enddo
741	enddo
742	ENDIF
743	enddo
744	!$OMP ENDDO NOWAIT
745	enddo
746	IF (prt_level > 5)
747	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
748	c=====================================================================
749	c PAS DE TEMPS D'ECRITURE
750	c=====================================================================
751	if (icum.eq.ncum) then
752	c=====================================================================
753
754	IF (prt_level > 5)
755	. WRITE(lunout,*)'Pas d ecriture'
756
757	jjb=jj_begin
758	jje=jj_end
759
760	c Normalisation
761	do iQ=1,nQ
762	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
763	do l=1,llm
764	Q_cum(:,jjb:jje,l,iQ)=Q_cum(:,jjb:jje,l,iQ)
765	. /masse_cum(:,jjb:jje,l)
766	enddo
767	!$OMP ENDDO NOWAIT
768	enddo
769
770	zz=1./REAL(ncum)
771
772	!$OMP MASTER
773	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
774	!$OMP END MASTER
775
776	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
777	DO l=1,llm
778	masse_cum(:,jjb:jje,l)=masse_cum(:,jjb:jje,l)*zz
779	flux_u_cum(:,jjb:jje,l)=flux_u_cum(:,jjb:jje,l)*zz
780	flux_uQ_cum(:,jjb:jje,l,:)=flux_uQ_cum(:,jjb:jje,l,:)*zz
781	dQ(:,jjb:jje,l,:)=dQ(:,jjb:jje,l,:)*zz
782	ENDDO
783	!$OMP ENDDO NOWAIT
784
785
786	IF (pole_sud) jje=jj_end-1
787	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
788	DO l=1,llm
789	flux_v_cum(:,jjb:jje,l)=flux_v_cum(:,jjb:jje,l)*zz
790	flux_vQ_cum(:,jjb:jje,l,:)=flux_vQ_cum(:,jjb:jje,l,:)*zz
791	ENDDO
792	!$OMP ENDDO
793
794	jjb=jj_begin
795	jje=jj_end
796
797
798	c A retravailler eventuellement
799	c division de dQ par la masse pour revenir aux bonnes grandeurs
800	do iQ=1,nQ
801	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
802	DO l=1,llm
803	dQ(:,jjb:jje,l,iQ)=dQ(:,jjb:jje,l,iQ)/masse_cum(:,jjb:jje,l)
804	ENDDO
805	!$OMP ENDDO NOWAIT
806	enddo
807
808	c=====================================================================
809	c Transport méridien
810	c=====================================================================
811
812	c cumul zonal des masses des mailles
813	c ----------------------------------
814	jjb=jj_begin
815	jje=jj_end
816	if (pole_sud) jje=jj_end-1
817
818	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
819	DO l=1,llm
820	zv(jjb:jje,l)=0.
821	zmasse(jjb:jje,l)=0.
822	ENDDO
823	!$OMP ENDDO NOWAIT
824
825	call Register_Hallo(masse_cum,ip1jmp1,llm,1,1,1,1,Req)
826	do iQ=1,nQ
827	call Register_Hallo(Q_cum(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
828	enddo
829
830	call SendRequest(Req)
831	!$OMP BARRIER
832	call WaitRequest(Req)
833	c$OMP BARRIER
834
835	call massbar_p(masse_cum,massebx,masseby)
836
837	jjb=jj_begin
838	jje=jj_end
839	if (pole_sud) jje=jj_end-1
840
841	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
842	do l=1,llm
843	do j=jjb,jje
844	do i=1,iim
845	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
846	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
847	enddo
848	zfactv(j,l)=cv(1,j)/zmasse(j,l)
849	enddo
850	enddo
851	!$OMP ENDDO
852
853	c print*,'3OK'
854	c --------------------------------------------------------------
855	c calcul de la moyenne zonale du transport :
856	c ------------------------------------------
857	c
858	c --
859	c TOT : la circulation totale [ vq ]
860	c
861	c - -
862	c MMC : mean meridional circulation [ v ] [ q ]
863	c
864	c ---- -- - -
865	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
866	c
867	c - * - * - - - -
868	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
869	c
870	c - -
871	c on utilise aussi l'intermediaire TMP : [ v q ]
872	c
873	c la variable zfactv transforme un transport meridien cumule
874	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
875	c
876	c --------------------------------------------------------------
877
878
879	c ----------------------------------------
880	c Transport dans le plan latitude-altitude
881	c ----------------------------------------
882
883	jjb=jj_begin
884	jje=jj_end
885	if (pole_sud) jje=jj_end-1
886
887	zvQ=0.
888	psiQ=0.
889	do iQ=1,nQ
890	!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
891	do l=1,llm
892	zvQtmp(:,l)=0.
893	do j=jjb,jje
894	c print*,'j,l,iQ=',j,l,iQ
895	c Calcul des moyennes zonales du transort total et de zvQtmp
896	do i=1,iim
897	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
898	s +flux_vQ_cum(i,j,l,iQ)
899	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
900	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
901	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
902	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
903	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
904	enddo
905	c print*,'aOK'
906	c Decomposition
907	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
908	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
909	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
910	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
911	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
912	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
913	enddo
914	enddo
915	!$OMP ENDDO NOWAIT
916	c fonction de courant meridienne pour la quantite Q
917	!$OMP BARRIER
918	!$OMP MASTER
919	do l=llm,1,-1
920	do j=jjb,jje
921	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
922	enddo
923	enddo
924	!$OMP END MASTER
925	!$OMP BARRIER
926	enddo ! of do iQ=1,nQ
927
928	c fonction de courant pour la circulation meridienne moyenne
929	!$OMP BARRIER
930	!$OMP MASTER
931	psi(jjb:jje,:)=0.
932	do l=llm,1,-1
933	do j=jjb,jje
934	psi(j,l)=psi(j,l+1)+zv(j,l)
935	zv(j,l)=zv(j,l)*zfactv(j,l)
936	enddo
937	enddo
938	!$OMP END MASTER
939	!$OMP BARRIER
940
941	c print*,'4OK'
942	c sorties proprement dites
943	!$OMP MASTER
944	if (i_sortie.eq.1) then
945	jjb=jj_begin
946	jje=jj_end
947	jjn=jj_nb
948	if (pole_sud) jje=jj_end-1
949	if (pole_sud) jjn=jj_nb-1
950
951	do iQ=1,nQ
952	do itr=1,ntr
953	call histwrite(fileid,znom(itr,iQ),itau,
954	s zvQ(jjb:jje,:,itr,iQ)
955	s ,jjn*llm,ndex3d)
956	enddo
957	call histwrite(fileid,'psi'//nom(iQ),
958	s itau,psiQ(jjb:jje,1:llm,iQ)
959	s ,jjn*llm,ndex3d)
960	enddo
961	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
962	s ,jjn*llm,ndex3d)
963	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
964	s ,jjn*llm,ndex3d)
965	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
966	call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
967	s jjn*llm,ndex3d)
968
969	endif
970
971
972	c -----------------
973	c Moyenne verticale
974	c -----------------
975
976	zamasse(jjb:jje)=0.
977	do l=1,llm
978	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
979	enddo
980
981	zavQ(jjb:jje,:,:)=0.
982	do iQ=1,nQ
983	do itr=2,ntr
984	do l=1,llm
985	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
986	s +zvQ(jjb:jje,l,itr,iQ)
987	s *zmasse(jjb:jje,l)
988	enddo
989	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
990	call histwrite(fileid,'a'//znom(itr,iQ),itau,
991	s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
992	enddo
993	enddo
994	!$OMP END MASTER
995	!$OMP BARRIER
996	c on doit pouvoir tracer systematiquement la fonction de courant.
997
998	c=====================================================================
999	c/////////////////////////////////////////////////////////////////////
1000	icum=0 !///////////////////////////////////////
1001	endif ! icum.eq.ncum !///////////////////////////////////////
1002	c/////////////////////////////////////////////////////////////////////
1003	c=====================================================================
1004	return
1005	end
1006

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

Download in other formats: