/[lmdze]/trunk/dyn3d/bilan_dyn.f
ViewVC logotype

Annotation of /trunk/dyn3d/bilan_dyn.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 30 - (hide annotations)
Thu Apr 1 09:07:28 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/bilan_dyn.f
File size: 16378 byte(s)
Imported Source files of the external library "IOIPSL_Lionel" into
"libf/IOIPSL".

Split "cray.f90" into "scopy.f90" and "ssum.f90".

Rewrote "leapfrog" in order to have a clearer algorithmic structure.

1 guez 3 !
2     ! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/bilan_dyn.F,v 1.5 2005/03/16 10:12:17 fairhead Exp $
3     !
4     SUBROUTINE bilan_dyn (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 guez 30 USE histcom
13     use calendar
14     use histwrite_m
15 guez 3 use dimens_m
16     use paramet_m
17     use comconst
18     use comvert
19     use comgeom, only: constang_2d, cu_2d, cv_2d, rlatv
20     use temps
21     use iniprint
22     use inigrads_m, only: inigrads
23    
24     IMPLICIT NONE
25    
26    
27     c====================================================================
28     c
29     c Sous-programme consacre à des diagnostics dynamiques de base
30     c
31     c
32     c De facon generale, les moyennes des scalaires Q sont ponderees par
33     c la masse.
34     c
35     c Les flux de masse sont eux simplement moyennes.
36     c
37     c====================================================================
38    
39     c Arguments :
40     c ===========
41    
42     integer ntrac
43     real dt_app,dt_cum
44     real ps(iip1,jjp1)
45     real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
46     real flux_u(iip1,jjp1,llm)
47     real flux_v(iip1,jjm,llm)
48     real teta(iip1,jjp1,llm)
49     real phi(iip1,jjp1,llm)
50     real ucov(iip1,jjp1,llm)
51     real vcov(iip1,jjm,llm)
52     real trac(iip1,jjp1,llm,ntrac)
53    
54     c Local :
55     c =======
56    
57     integer icum,ncum
58     logical first
59     real zz,zqy,zfactv(jjm,llm)
60    
61     integer nQ
62     parameter (nQ=7)
63    
64    
65     cym character*6 nom(nQ)
66     cym character*6 unites(nQ)
67     character*6,save :: nom(nQ)
68     character*6,save :: unites(nQ)
69    
70     character*10 file
71     integer ifile
72     parameter (ifile=4)
73    
74     integer itemp,igeop,iecin,iang,iu,iovap,iun
75     integer i_sortie
76    
77     save first,icum,ncum
78     save itemp,igeop,iecin,iang,iu,iovap,iun
79     save i_sortie
80    
81     real time
82     integer itau
83     save time,itau
84     data time,itau/0.,0/
85    
86     data first/.true./
87     data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
88     data i_sortie/1/
89    
90     real ww
91    
92     c variables dynamiques intermédiaires
93     REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
94     REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
95     REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
96     REAL vorpot(iip1,jjm,llm)
97     REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
98     REAL bern(iip1,jjp1,llm)
99    
100     c champ contenant les scalaires advectés.
101     real Q(iip1,jjp1,llm,nQ)
102    
103     c champs cumulés
104     real ps_cum(iip1,jjp1)
105     real masse_cum(iip1,jjp1,llm)
106     real flux_u_cum(iip1,jjp1,llm)
107     real flux_v_cum(iip1,jjm,llm)
108     real Q_cum(iip1,jjp1,llm,nQ)
109     real flux_uQ_cum(iip1,jjp1,llm,nQ)
110     real flux_vQ_cum(iip1,jjm,llm,nQ)
111     real flux_wQ_cum(iip1,jjp1,llm,nQ)
112     real dQ(iip1,jjp1,llm,nQ)
113    
114     save ps_cum,masse_cum,flux_u_cum,flux_v_cum
115     save Q_cum,flux_uQ_cum,flux_vQ_cum
116    
117     c champs de tansport en moyenne zonale
118     integer ntr,itr
119     parameter (ntr=5)
120    
121     cym character*10 znom(ntr,nQ)
122     cym character*20 znoml(ntr,nQ)
123     cym character*10 zunites(ntr,nQ)
124     character*10,save :: znom(ntr,nQ)
125     character*20,save :: znoml(ntr,nQ)
126     character*10,save :: zunites(ntr,nQ)
127    
128     integer iave,itot,immc,itrs,istn
129     data iave,itot,immc,itrs,istn/1,2,3,4,5/
130     character*3 ctrs(ntr)
131     data ctrs/' ','TOT','MMC','TRS','STN'/
132    
133     real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
134     real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
135     real zmasse(jjm,llm),zamasse(jjm)
136    
137     real zv(jjm,llm),psi(jjm,llm+1)
138    
139     integer i,j,l,iQ
140    
141    
142     c Initialisation du fichier contenant les moyennes zonales.
143     c ---------------------------------------------------------
144    
145     integer fileid
146     integer thoriid, zvertiid
147     save fileid
148    
149     integer ndex3d(jjm*llm)
150    
151     C Variables locales
152     C
153     real zjulian
154     character*3 str
155     character*10 ctrac
156     integer ii,jj
157     integer zan, dayref
158     C
159     real rlong(jjm),rlatg(jjm)
160    
161     !!print *, "Call sequence information: bilan_dyn"
162    
163     c=====================================================================
164     c Initialisation
165     c=====================================================================
166    
167     time=time+dt_app
168     itau=itau+1
169    
170     if (first) then
171    
172    
173     icum=0
174     c initialisation des fichiers
175     first=.false.
176     c ncum est la frequence de stokage en pas de temps
177     ncum=dt_cum/dt_app
178     if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
179 guez 12 print *,
180 guez 3 . 'Pb : le pas de cumule doit etre multiple du pas'
181 guez 12 print *,'dt_app=',dt_app
182     print *,'dt_cum=',dt_cum
183 guez 3 stop
184     endif
185    
186     if (i_sortie.eq.1) then
187     file='dynzon'
188 guez 20 call inigrads(ifile ,(/0./),180./pi,0.,0.,rlatv,-90.,90.,
189     $ 180./pi ,presnivs,1. ,dt_cum,file,'dyn_zon ')
190 guez 3 endif
191    
192     nom(itemp)='T'
193     nom(igeop)='gz'
194     nom(iecin)='K'
195     nom(iang)='ang'
196     nom(iu)='u'
197     nom(iovap)='ovap'
198     nom(iun)='un'
199    
200     unites(itemp)='K'
201     unites(igeop)='m2/s2'
202     unites(iecin)='m2/s2'
203     unites(iang)='ang'
204     unites(iu)='m/s'
205     unites(iovap)='kg/kg'
206     unites(iun)='un'
207    
208    
209     c Initialisation du fichier contenant les moyennes zonales.
210     c ---------------------------------------------------------
211    
212     zan = annee_ref
213     dayref = day_ref
214     CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
215    
216     rlong=0.
217     rlatg=rlatv*180./pi
218    
219 guez 27 call histbeg_totreg('dynzon', rlong(:1), rlatg,
220 guez 3 . 1, 1, 1, jjm,
221 guez 27 . itau_dyn, zjulian, dt_cum, thoriid, fileid)
222 guez 3
223     C
224     C Appel a histvert pour la grille verticale
225     C
226     call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
227     . llm, presnivs, zvertiid)
228     C
229     C Appels a histdef pour la definition des variables a sauvegarder
230     do iQ=1,nQ
231     do itr=1,ntr
232     if(itr.eq.1) then
233     znom(itr,iQ)=nom(iQ)
234     znoml(itr,iQ)=nom(iQ)
235     zunites(itr,iQ)=unites(iQ)
236     else
237     znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
238     znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
239     zunites(itr,iQ)='m/s * '//unites(iQ)
240     endif
241     enddo
242     enddo
243    
244     c Declarations des champs avec dimension verticale
245     c print*,'1HISTDEF'
246     do iQ=1,nQ
247     do itr=1,ntr
248     IF (prt_level > 5)
249 guez 12 . print *,'var ',itr,iQ
250 guez 3 . ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
251     call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
252     . zunites(itr,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
253 guez 15 . 'ave(X)',dt_cum,dt_cum)
254 guez 3 enddo
255     c Declarations pour les fonctions de courant
256     c print*,'2HISTDEF'
257     call histdef(fileid,'psi'//nom(iQ)
258     . ,'stream fn. '//znoml(itot,iQ),
259     . zunites(itot,iQ),1,jjm,thoriid,llm,1,llm,zvertiid,
260 guez 15 . 'ave(X)',dt_cum,dt_cum)
261 guez 3 enddo
262    
263    
264     c Declarations pour les champs de transport d'air
265     c print*,'3HISTDEF'
266     call histdef(fileid, 'masse', 'masse',
267     . 'kg', 1, jjm, thoriid, llm, 1, llm, zvertiid,
268 guez 15 . 'ave(X)', dt_cum, dt_cum)
269 guez 3 call histdef(fileid, 'v', 'v',
270     . 'm/s', 1, jjm, thoriid, llm, 1, llm, zvertiid,
271 guez 15 . 'ave(X)', dt_cum, dt_cum)
272 guez 3 c Declarations pour les fonctions de courant
273     c print*,'4HISTDEF'
274     call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
275     . 1,jjm,thoriid,llm,1,llm,zvertiid,
276 guez 15 . 'ave(X)',dt_cum,dt_cum)
277 guez 3
278    
279     c Declaration des champs 1D de transport en latitude
280     c print*,'5HISTDEF'
281     do iQ=1,nQ
282     do itr=2,ntr
283     call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
284     . zunites(itr,iQ),1,jjm,thoriid,1,1,1,-99,
285 guez 15 . 'ave(X)',dt_cum,dt_cum)
286 guez 3 enddo
287     enddo
288    
289    
290     c print*,'8HISTDEF'
291     CALL histend(fileid)
292    
293    
294     endif
295    
296    
297     c=====================================================================
298     c Calcul des champs dynamiques
299     c ----------------------------
300    
301     c énergie cinétique
302     ucont(:,:,:)=0
303     CALL covcont(llm,ucov,vcov,ucont,vcont)
304     CALL enercin(vcov,ucov,vcont,ucont,ecin)
305    
306     c moment cinétique
307     do l=1,llm
308     ang(:,:,l)=ucov(:,:,l)+constang_2d(:,:)
309     unat(:,:,l)=ucont(:,:,l)*cu_2d(:,:)
310     enddo
311    
312     Q(:,:,:,itemp)=teta(:,:,:)*pk(:,:,:)/cpp
313     Q(:,:,:,igeop)=phi(:,:,:)
314     Q(:,:,:,iecin)=ecin(:,:,:)
315     Q(:,:,:,iang)=ang(:,:,:)
316     Q(:,:,:,iu)=unat(:,:,:)
317     Q(:,:,:,iovap)=trac(:,:,:,1)
318     Q(:,:,:,iun)=1.
319    
320    
321     c=====================================================================
322     c Cumul
323     c=====================================================================
324     c
325     if(icum.EQ.0) then
326     ps_cum=0.
327     masse_cum=0.
328     flux_u_cum=0.
329     flux_v_cum=0.
330     Q_cum=0.
331     flux_vQ_cum=0.
332     flux_uQ_cum=0.
333     endif
334    
335     IF (prt_level > 5)
336 guez 12 . print *,'dans bilan_dyn ',icum,'->',icum+1
337 guez 3 icum=icum+1
338    
339     c accumulation des flux de masse horizontaux
340     ps_cum=ps_cum+ps
341     masse_cum=masse_cum+masse
342     flux_u_cum=flux_u_cum+flux_u
343     flux_v_cum=flux_v_cum+flux_v
344     do iQ=1,nQ
345     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)+Q(:,:,:,iQ)*masse(:,:,:)
346     enddo
347    
348     c=====================================================================
349     c FLUX ET TENDANCES
350     c=====================================================================
351    
352     c Flux longitudinal
353     c -----------------
354     do iQ=1,nQ
355     do l=1,llm
356     do j=1,jjp1
357     do i=1,iim
358     flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
359     s +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
360     enddo
361     flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
362     enddo
363     enddo
364     enddo
365    
366     c flux méridien
367     c -------------
368     do iQ=1,nQ
369     do l=1,llm
370     do j=1,jjm
371     do i=1,iip1
372     flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
373     s +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
374     enddo
375     enddo
376     enddo
377     enddo
378    
379    
380     c tendances
381     c ---------
382    
383     c convergence horizontale
384     call convflu(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
385    
386     c calcul de la vitesse verticale
387     call convmas(flux_u_cum,flux_v_cum,convm)
388     CALL vitvert(convm,w)
389    
390     do iQ=1,nQ
391     do l=1,llm-1
392     do j=1,jjp1
393     do i=1,iip1
394     ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
395     dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
396     dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
397     enddo
398     enddo
399     enddo
400     enddo
401     IF (prt_level > 5)
402 guez 12 . print *,'Apres les calculs fait a chaque pas'
403 guez 3 c=====================================================================
404     c PAS DE TEMPS D'ECRITURE
405     c=====================================================================
406     if (icum.eq.ncum) then
407     c=====================================================================
408    
409     IF (prt_level > 5)
410 guez 12 . print *,'Pas d ecriture'
411 guez 3
412     c Normalisation
413     do iQ=1,nQ
414     Q_cum(:,:,:,iQ)=Q_cum(:,:,:,iQ)/masse_cum(:,:,:)
415     enddo
416     zz=1./float(ncum)
417     ps_cum=ps_cum*zz
418     masse_cum=masse_cum*zz
419     flux_u_cum=flux_u_cum*zz
420     flux_v_cum=flux_v_cum*zz
421     flux_uQ_cum=flux_uQ_cum*zz
422     flux_vQ_cum=flux_vQ_cum*zz
423     dQ=dQ*zz
424    
425    
426     c A retravailler eventuellement
427     c division de dQ par la masse pour revenir aux bonnes grandeurs
428     do iQ=1,nQ
429     dQ(:,:,:,iQ)=dQ(:,:,:,iQ)/masse_cum(:,:,:)
430     enddo
431    
432     c=====================================================================
433     c Transport méridien
434     c=====================================================================
435    
436     c cumul zonal des masses des mailles
437     c ----------------------------------
438     zv=0.
439     zmasse=0.
440     call massbar(masse_cum,massebx,masseby)
441     do l=1,llm
442     do j=1,jjm
443     do i=1,iim
444     zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
445     zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
446     enddo
447     zfactv(j,l)=cv_2d(1,j)/zmasse(j,l)
448     enddo
449     enddo
450    
451     c print*,'3OK'
452     c --------------------------------------------------------------
453     c calcul de la moyenne zonale du transport :
454     c ------------------------------------------
455     c
456     c --
457     c TOT : la circulation totale [ vq ]
458     c
459     c - -
460     c MMC : mean meridional circulation [ v ] [ q ]
461     c
462     c ---- -- - -
463     c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
464     c
465     c - * - * - - - -
466     c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
467     c
468     c - -
469     c on utilise aussi l'intermediaire TMP : [ v q ]
470     c
471     c la variable zfactv transforme un transport meridien cumule
472     c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
473     c
474     c --------------------------------------------------------------
475    
476    
477     c ----------------------------------------
478     c Transport dans le plan latitude-altitude
479     c ----------------------------------------
480    
481     zvQ=0.
482     psiQ=0.
483     do iQ=1,nQ
484     zvQtmp=0.
485     do l=1,llm
486     do j=1,jjm
487     c print*,'j,l,iQ=',j,l,iQ
488     c Calcul des moyennes zonales du transort total et de zvQtmp
489     do i=1,iim
490     zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
491     s +flux_vQ_cum(i,j,l,iQ)
492     zqy= 0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
493     s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
494     zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
495     s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
496     zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
497     enddo
498     c print*,'aOK'
499     c Decomposition
500     zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
501     zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
502     zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
503     zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
504     zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
505     zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
506     enddo
507     enddo
508     c fonction de courant meridienne pour la quantite Q
509     do l=llm,1,-1
510     do j=1,jjm
511     psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
512     enddo
513     enddo
514     enddo
515    
516     c fonction de courant pour la circulation meridienne moyenne
517     psi=0.
518     do l=llm,1,-1
519     do j=1,jjm
520     psi(j,l)=psi(j,l+1)+zv(j,l)
521     zv(j,l)=zv(j,l)*zfactv(j,l)
522     enddo
523     enddo
524    
525     c print*,'4OK'
526     c sorties proprement dites
527     if (i_sortie.eq.1) then
528     do iQ=1,nQ
529     do itr=1,ntr
530 guez 15 call histwrite(fileid,znom(itr,iQ),itau,zvQ(:,:,itr,iQ))
531 guez 3 enddo
532 guez 15 call histwrite(fileid,'psi'//nom(iQ),itau,psiQ(:,1:llm,iQ))
533 guez 3 enddo
534    
535 guez 15 call histwrite(fileid,'masse',itau,zmasse)
536     call histwrite(fileid,'v',itau,zv)
537 guez 3 psi=psi*1.e-9
538 guez 15 call histwrite(fileid,'psi',itau,psi(:,1:llm))
539 guez 3
540     endif
541    
542    
543     c -----------------
544     c Moyenne verticale
545     c -----------------
546    
547     zamasse=0.
548     do l=1,llm
549     zamasse(:)=zamasse(:)+zmasse(:,l)
550     enddo
551     zavQ=0.
552     do iQ=1,nQ
553     do itr=2,ntr
554     do l=1,llm
555     zavQ(:,itr,iQ)=zavQ(:,itr,iQ)+zvQ(:,l,itr,iQ)*zmasse(:,l)
556     enddo
557     zavQ(:,itr,iQ)=zavQ(:,itr,iQ)/zamasse(:)
558 guez 15 call histwrite(fileid,'a'//znom(itr,iQ),itau,zavQ(:,itr,iQ))
559 guez 3 enddo
560     enddo
561    
562     c on doit pouvoir tracer systematiquement la fonction de courant.
563    
564     c=====================================================================
565     c/////////////////////////////////////////////////////////////////////
566     icum=0 !///////////////////////////////////////
567     endif ! icum.eq.ncum !///////////////////////////////////////
568     c/////////////////////////////////////////////////////////////////////
569     c=====================================================================
570    
571     return
572     end

  ViewVC Help
Powered by ViewVC 1.1.21