GCC Code Coverage Report


Directory: ./
File: dyn/bilan_dyn.f
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 185 0.0%
Branches: 0 248 0.0%

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