GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: dyn3d/bilan_dyn.F Lines: 0 185 0.0 %
Date: 2023-06-30 12:51:15 Branches: 0 248 0.0 %

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