source: LMDZ4/branches/LMDZ4-dev-20091210/libf/dyn3dpar/bilan_dyn_p.F @ 5444

Last change on this file since 5444 was 1186, checked in by Ehouarn Millour, 16 years ago

Cleanup around IOIPSL, so that LMDZ dynamics may be used without IOIPSL.

  • moved ersatz IOIPSL routines (ioipsl_* , taken from IOIPSLv2_1_8, so that 'getin' function may be used even if not using the IOIPSL library) from dyn3d/dyn3dpar to bibio.
  • enclosed 'use ioipsl' instruction with #ifdef CPP_IOIPSL cpp keys.

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.9 KB
Line 
1!
2! $Id: bilan_dyn_p.F 1186 2009-06-18 09:20:44Z fhourdin $
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
7c   AFAIRE
8c   Prevoir en champ nq+1 le diagnostique de l'energie
9c   en faisant Qzon=Cv T + L * ...
10c             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
29c====================================================================
30c
31c   Sous-programme consacre à des diagnostics dynamiques de base
32c
33c
34c   De facon generale, les moyennes des scalaires Q sont ponderees par
35c   la masse.
36c
37c   Les flux de masse sont eux simplement moyennes.
38c
39c====================================================================
40
41c   Arguments :
42c   ===========
43
44      integer ntrac
45      real dt_app,dt_cum
46      real ps(iip1,jjp1)
47      real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
48      real flux_u(iip1,jjp1,llm)
49      real flux_v(iip1,jjm,llm)
50      real teta(iip1,jjp1,llm)
51      real phi(iip1,jjp1,llm)
52      real ucov(iip1,jjp1,llm)
53      real vcov(iip1,jjm,llm)
54      real trac(iip1,jjp1,llm,ntrac)
55
56c   Local :
57c   =======
58
59      integer icum,ncum
60      logical first
61      real zz,zqy,zfactv(jjm,llm)
62
63      integer nQ
64      parameter (nQ=7)
65
66
67cym      character*6 nom(nQ)
68cym      character*6 unites(nQ)
69      character*6,save :: nom(nQ)
70      character*6,save :: unites(nQ)
71
72      character*10 file
73      integer ifile
74      parameter (ifile=4)
75
76      integer itemp,igeop,iecin,iang,iu,iovap,iun
77      integer i_sortie
78
79      save first,icum,ncum
80      save itemp,igeop,iecin,iang,iu,iovap,iun
81      save i_sortie
82
83      real time
84      integer itau
85      save time,itau
86      data time,itau/0.,0/
87
88      data first/.true./
89      data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
90      data i_sortie/1/
91
92      real ww
93
94c   variables dynamiques intermédiaires
95      REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
96      REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
97      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
98      REAL vorpot(iip1,jjm,llm)
99      REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
100      REAL bern(iip1,jjp1,llm)
101
102c   champ contenant les scalaires advectés.
103      real Q(iip1,jjp1,llm,nQ)
104   
105c   champs cumulés
106      real ps_cum(iip1,jjp1)
107      real masse_cum(iip1,jjp1,llm)
108      real flux_u_cum(iip1,jjp1,llm)
109      real flux_v_cum(iip1,jjm,llm)
110      real Q_cum(iip1,jjp1,llm,nQ)
111      real flux_uQ_cum(iip1,jjp1,llm,nQ)
112      real flux_vQ_cum(iip1,jjm,llm,nQ)
113      real flux_wQ_cum(iip1,jjp1,llm,nQ)
114      real dQ(iip1,jjp1,llm,nQ)
115
116      save ps_cum,masse_cum,flux_u_cum,flux_v_cum
117      save Q_cum,flux_uQ_cum,flux_vQ_cum
118
119c   champs de tansport en moyenne zonale
120      integer ntr,itr
121      parameter (ntr=5)
122
123cym      character*10 znom(ntr,nQ)
124cym      character*20 znoml(ntr,nQ)
125cym      character*10 zunites(ntr,nQ)
126      character*10,save :: znom(ntr,nQ)
127      character*20,save :: znoml(ntr,nQ)
128      character*10,save :: zunites(ntr,nQ)
129
130      integer iave,itot,immc,itrs,istn
131      data iave,itot,immc,itrs,istn/1,2,3,4,5/
132      character*3 ctrs(ntr)
133      data ctrs/'  ','TOT','MMC','TRS','STN'/
134
135      real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
136      real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
137      real zmasse(jjm,llm),zamasse(jjm)
138
139      real zv(jjm,llm),psi(jjm,llm+1)
140
141      integer i,j,l,iQ
142
143
144c   Initialisation du fichier contenant les moyennes zonales.
145c   ---------------------------------------------------------
146
147      character*10 infile
148
149      integer fileid
150      integer thoriid, zvertiid
151      save fileid
152
153      integer ndex3d(jjm*llm)
154
155C   Variables locales
156C
157      integer tau0
158      real zjulian
159      character*3 str
160      character*10 ctrac
161      integer ii,jj
162      integer zan, dayref
163C
164      real rlong(jjm),rlatg(jjm)
165      integer :: jjb,jje,jjn,ijb,ije
166      type(Request) :: Req
167
168! definition du domaine d'ecriture pour le rebuild
169
170      INTEGER,DIMENSION(1) :: ddid
171      INTEGER,DIMENSION(1) :: dsg
172      INTEGER,DIMENSION(1) :: dsl
173      INTEGER,DIMENSION(1) :: dpf
174      INTEGER,DIMENSION(1) :: dpl
175      INTEGER,DIMENSION(1) :: dhs
176      INTEGER,DIMENSION(1) :: dhe
177     
178      INTEGER :: bilan_dyn_domain_id
179
180
181c=====================================================================
182c   Initialisation
183c=====================================================================
184      ndex3d=0
185      if (adjust) return
186     
187      time=time+dt_app
188      itau=itau+1
189
190      if (first) then
191
192
193        icum=0
194c       initialisation des fichiers
195        first=.false.
196c   ncum est la frequence de stokage en pas de temps
197        ncum=dt_cum/dt_app
198        if (abs(ncum*dt_app-dt_cum).gt.1.e-5*dt_app) then
199           WRITE(lunout,*)
200     .            'Pb : le pas de cumule doit etre multiple du pas'
201           WRITE(lunout,*)'dt_app=',dt_app
202           WRITE(lunout,*)'dt_cum=',dt_cum
203           stop
204        endif
205
206        if (i_sortie.eq.1) then
207         file='dynzon'
208         if (mpi_rank==0) then
209         call inigrads(ifile,1
210     s  ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
211     s  ,llm,presnivs,1.
212     s  ,dt_cum,file,'dyn_zon ')
213         endif
214        endif
215
216        nom(itemp)='T'
217        nom(igeop)='gz'
218        nom(iecin)='K'
219        nom(iang)='ang'
220        nom(iu)='u'
221        nom(iovap)='ovap'
222        nom(iun)='un'
223
224        unites(itemp)='K'
225        unites(igeop)='m2/s2'
226        unites(iecin)='m2/s2'
227        unites(iang)='ang'
228        unites(iu)='m/s'
229        unites(iovap)='kg/kg'
230        unites(iun)='un'
231
232
233c   Initialisation du fichier contenant les moyennes zonales.
234c   ---------------------------------------------------------
235
236      infile='dynzon'
237
238      zan = annee_ref
239      dayref = day_ref
240      CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
241      tau0 = itau_dyn
242     
243      rlong=0.
244      rlatg=rlatv*180./pi
245
246      jjb=jj_begin
247      jje=jj_end
248      jjn=jj_nb
249      IF (pole_sud) THEN
250        jjn=jj_nb-1
251        jje=jj_end-1
252      ENDIF
253
254      ddid=(/ 2 /)
255      dsg=(/ jjm /)
256      dsl=(/ jjn /)
257      dpf=(/ jjb /)
258      dpl=(/ jje /)
259      dhs=(/ 0 /)
260      dhe=(/ 0 /)
261
262      call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
263     .                 'box',bilan_dyn_domain_id)
264       
265      call histbeg(trim(infile),
266     .             1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
267     .             1, 1, 1, jjn,
268     .             tau0, zjulian, dt_cum, thoriid, fileid,
269     .             bilan_dyn_domain_id)
270
271C
272C  Appel a histvert pour la grille verticale
273C
274      call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
275     .              llm, presnivs, zvertiid)
276C
277C  Appels a histdef pour la definition des variables a sauvegarder
278      do iQ=1,nQ
279         do itr=1,ntr
280            if(itr.eq.1) then
281               znom(itr,iQ)=nom(iQ)
282               znoml(itr,iQ)=nom(iQ)
283               zunites(itr,iQ)=unites(iQ)
284            else
285               znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
286               znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
287               zunites(itr,iQ)='m/s * '//unites(iQ)
288            endif
289         enddo
290      enddo
291
292c   Declarations des champs avec dimension verticale
293c      print*,'1HISTDEF'
294      do iQ=1,nQ
295         do itr=1,ntr
296      IF (prt_level > 5)
297     . WRITE(lunout,*)'var ',itr,iQ
298     .      ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
299            call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
300     .        zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
301     .        32,'ave(X)',dt_cum,dt_cum)
302         enddo
303c   Declarations pour les fonctions de courant
304c      print*,'2HISTDEF'
305          call histdef(fileid,'psi'//nom(iQ)
306     .      ,'stream fn. '//znoml(itot,iQ),
307     .      zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
308     .      32,'ave(X)',dt_cum,dt_cum)
309      enddo
310
311
312c   Declarations pour les champs de transport d'air
313c      print*,'3HISTDEF'
314      call histdef(fileid, 'masse', 'masse',
315     .             'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
316     .             32, 'ave(X)', dt_cum, dt_cum)
317      call histdef(fileid, 'v', 'v',
318     .             'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
319     .             32, 'ave(X)', dt_cum, dt_cum)
320c   Declarations pour les fonctions de courant
321c      print*,'4HISTDEF'
322          call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
323     .      1,jjn,thoriid,llm,1,llm,zvertiid,
324     .      32,'ave(X)',dt_cum,dt_cum)
325
326
327c   Declaration des champs 1D de transport en latitude
328c      print*,'5HISTDEF'
329      do iQ=1,nQ
330         do itr=2,ntr
331            call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
332     .        zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
333     .        32,'ave(X)',dt_cum,dt_cum)
334         enddo
335      enddo
336
337
338c      print*,'8HISTDEF'
339               CALL histend(fileid)
340
341
342      endif
343
344
345c=====================================================================
346c   Calcul des champs dynamiques
347c   ----------------------------
348
349      jjb=jj_begin
350      jje=jj_end
351   
352c   énergie cinétique
353      ucont(:,jjb:jje,:)=0
354
355      call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Req)
356      call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Req)
357      call SendRequest(Req)
358      call WaitRequest(Req)
359
360      CALL covcont_p(llm,ucov,vcov,ucont,vcont)
361      CALL enercin_p(vcov,ucov,vcont,ucont,ecin)
362
363c   moment cinétique
364      do l=1,llm
365         ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
366         unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
367      enddo
368
369      Q(:,jjb:jje,:,itemp)=teta(:,jjb:jje,:)*pk(:,jjb:jje,:)/cpp
370      Q(:,jjb:jje,:,igeop)=phi(:,jjb:jje,:)
371      Q(:,jjb:jje,:,iecin)=ecin(:,jjb:jje,:)
372      Q(:,jjb:jje,:,iang)=ang(:,jjb:jje,:)
373      Q(:,jjb:jje,:,iu)=unat(:,jjb:jje,:)
374      Q(:,jjb:jje,:,iovap)=trac(:,jjb:jje,:,1)
375      Q(:,jjb:jje,:,iun)=1.
376
377
378c=====================================================================
379c   Cumul
380c=====================================================================
381c
382      if(icum.EQ.0) then
383         jjb=jj_begin
384         jje=jj_end
385
386         ps_cum(:,jjb:jje)=0.
387         masse_cum(:,jjb:jje,:)=0.
388         flux_u_cum(:,jjb:jje,:)=0.
389         Q_cum(:,jjb:jje,:,:)=0.
390         flux_uQ_cum(:,jjb:jje,:,:)=0.
391         flux_v_cum(:,jjb:jje,:)=0.
392         if (pole_sud) jje=jj_end-1
393         flux_v_cum(:,jjb:jje,:)=0.
394         flux_vQ_cum(:,jjb:jje,:,:)=0.
395      endif
396
397      IF (prt_level > 5)
398     . WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
399      icum=icum+1
400
401c   accumulation des flux de masse horizontaux
402      jjb=jj_begin
403      jje=jj_end
404
405      ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
406      masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
407      flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
408     .                       +flux_u(:,jjb:jje,:)
409      if (pole_sud) jje=jj_end-1
410      flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)
411     .                         +flux_v(:,jjb:jje,:)
412
413      jjb=jj_begin
414      jje=jj_end
415
416      do iQ=1,nQ
417        Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
418     .                       +Q(:,jjb:jje,:,iQ)*masse(:,jjb:jje,:)
419      enddo
420
421c=====================================================================
422c  FLUX ET TENDANCES
423c=====================================================================
424
425c   Flux longitudinal
426c   -----------------
427      do iQ=1,nQ
428         do l=1,llm
429            do j=jjb,jje
430               do i=1,iim
431                  flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
432     s            +flux_u(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
433               enddo
434               flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
435            enddo
436         enddo
437      enddo
438
439c    flux méridien
440c    -------------
441      do iQ=1,nQ
442        call Register_Hallo(Q(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
443      enddo
444      call SendRequest(Req)
445      call WaitRequest(Req)
446     
447      jjb=jj_begin
448      jje=jj_end
449      if (pole_sud) jje=jj_end-1
450     
451      do iQ=1,nQ
452         do l=1,llm
453            do j=jjb,jje
454               do i=1,iip1
455                  flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
456     s            +flux_v(i,j,l)*0.5*(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
457               enddo
458            enddo
459         enddo
460      enddo
461
462
463c    tendances
464c    ---------
465
466c   convergence horizontale
467      call Register_Hallo(flux_uQ_cum,ip1jmp1,llm,2,2,2,2,Req)
468      call Register_Hallo(flux_vQ_cum,ip1jm,llm,2,2,2,2,Req)
469      call SendRequest(Req)
470      call WaitRequest(Req)
471
472      call  convflu_p(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
473
474c   calcul de la vitesse verticale
475      call Register_Hallo(flux_u_cum,ip1jmp1,llm,2,2,2,2,Req)
476      call Register_Hallo(flux_v_cum,ip1jm,llm,2,2,2,2,Req)
477      call SendRequest(Req)
478      call WaitRequest(Req)
479
480      call convmas_p(flux_u_cum,flux_v_cum,convm)
481      CALL vitvert_p(convm,w)
482
483      jjb=jj_begin
484      jje=jj_end
485
486      do iQ=1,nQ
487         do l=1,llm-1
488            do j=jjb,jje
489               do i=1,iip1
490                  ww=-0.5*w(i,j,l+1)*(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
491                  dQ(i,j,l  ,iQ)=dQ(i,j,l  ,iQ)-ww
492                  dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
493               enddo
494            enddo
495         enddo
496      enddo
497      IF (prt_level > 5)
498     . WRITE(lunout,*)'Apres les calculs fait a chaque pas'
499c=====================================================================
500c   PAS DE TEMPS D'ECRITURE
501c=====================================================================
502      if (icum.eq.ncum) then
503c=====================================================================
504
505      IF (prt_level > 5)
506     . WRITE(lunout,*)'Pas d ecriture'
507
508c   Normalisation
509      do iQ=1,nQ
510         Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
511     .                        /masse_cum(:,jjb:jje,:)
512      enddo
513      zz=1./float(ncum)
514
515      jjb=jj_begin
516      jje=jj_end
517
518      ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
519      masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)*zz
520      flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)*zz
521      flux_uQ_cum(:,jjb:jje,:,:)=flux_uQ_cum(:,jjb:jje,:,:)*zz
522      dQ(:,jjb:jje,:,:)=dQ(:,jjb:jje,:,:)*zz
523     
524      IF (pole_sud) jje=jj_end-1
525      flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)*zz
526      flux_vQ_cum(:,jjb:jje,:,:)=flux_vQ_cum(:,jjb:jje,:,:)*zz
527
528      jjb=jj_begin
529      jje=jj_end
530
531
532c   A retravailler eventuellement
533c   division de dQ par la masse pour revenir aux bonnes grandeurs
534      do iQ=1,nQ
535         dQ(:,jjb:jje,:,iQ)=dQ(:,jjb:jje,:,iQ)/masse_cum(:,jjb:jje,:)
536      enddo
537 
538c=====================================================================
539c   Transport méridien
540c=====================================================================
541
542c   cumul zonal des masses des mailles
543c   ----------------------------------
544      jjb=jj_begin
545      jje=jj_end
546      if (pole_sud) jje=jj_end-1
547
548      zv(jjb:jje,:)=0.
549      zmasse(jjb:jje,:)=0.
550
551      call Register_Hallo(masse_cum,ip1jmp1,llm,1,1,1,1,Req)
552      do iQ=1,nQ
553        call Register_Hallo(Q_cum(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
554      enddo
555
556      call SendRequest(Req)
557      call WaitRequest(Req)
558
559      call massbar_p(masse_cum,massebx,masseby)
560     
561      jjb=jj_begin
562      jje=jj_end
563      if (pole_sud) jje=jj_end-1
564     
565      do l=1,llm
566         do j=jjb,jje
567            do i=1,iim
568               zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
569               zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
570            enddo
571            zfactv(j,l)=cv(1,j)/zmasse(j,l)
572         enddo
573      enddo
574
575c     print*,'3OK'
576c   --------------------------------------------------------------
577c   calcul de la moyenne zonale du transport :
578c   ------------------------------------------
579c
580c                                     --
581c TOT : la circulation totale       [ vq ]
582c
583c                                      -     -
584c MMC : mean meridional circulation [ v ] [ q ]
585c
586c                                     ----      --       - -
587c TRS : transitoires                [ v'q'] = [ vq ] - [ v q ]
588c
589c                                     - * - *       - -       -     -
590c STT : stationaires                [ v   q   ] = [ v q ] - [ v ] [ q ]
591c
592c                                              - -
593c    on utilise aussi l'intermediaire TMP :  [ v q ]
594c
595c    la variable zfactv transforme un transport meridien cumule
596c    en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
597c
598c   --------------------------------------------------------------
599
600
601c   ----------------------------------------
602c   Transport dans le plan latitude-altitude
603c   ----------------------------------------
604
605      jjb=jj_begin
606      jje=jj_end
607      if (pole_sud) jje=jj_end-1
608     
609      zvQ=0.
610      psiQ=0.
611      do iQ=1,nQ
612         zvQtmp=0.
613         do l=1,llm
614            do j=jjb,jje
615c              print*,'j,l,iQ=',j,l,iQ
616c   Calcul des moyennes zonales du transort total et de zvQtmp
617               do i=1,iim
618                  zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
619     s                            +flux_vQ_cum(i,j,l,iQ)
620                  zqy=      0.5*(Q_cum(i,j,l,iQ)*masse_cum(i,j,l)+
621     s                           Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
622                  zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
623     s             /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
624                  zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
625               enddo
626c              print*,'aOK'
627c   Decomposition
628               zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
629               zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
630               zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
631               zvQ(j,l,immc,iQ)=zv(j,l)*zvQ(j,l,iave,iQ)*zfactv(j,l)
632               zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
633               zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
634            enddo
635         enddo
636c   fonction de courant meridienne pour la quantite Q
637         do l=llm,1,-1
638            do j=jjb,jje
639               psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
640            enddo
641         enddo
642      enddo
643
644c   fonction de courant pour la circulation meridienne moyenne
645      psi(jjb:jje,:)=0.
646      do l=llm,1,-1
647         do j=jjb,jje
648            psi(j,l)=psi(j,l+1)+zv(j,l)
649            zv(j,l)=zv(j,l)*zfactv(j,l)
650         enddo
651      enddo
652
653c     print*,'4OK'
654c   sorties proprement dites
655      if (i_sortie.eq.1) then
656      jjb=jj_begin
657      jje=jj_end
658      jjn=jj_nb
659      if (pole_sud) jje=jj_end-1
660      if (pole_sud) jjn=jj_nb-1
661     
662      do iQ=1,nQ
663         do itr=1,ntr
664            call histwrite(fileid,znom(itr,iQ),itau,
665     s                     zvQ(jjb:jje,:,itr,iQ)
666     s                     ,jjn*llm,ndex3d)
667         enddo
668         call histwrite(fileid,'psi'//nom(iQ),
669     s                  itau,psiQ(jjb:jje,1:llm,iQ)
670     s                  ,jjn*llm,ndex3d)
671      enddo
672
673      call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
674     s   ,jjn*llm,ndex3d)
675      call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
676     s   ,jjn*llm,ndex3d)
677      psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
678      call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
679     s               jjn*llm,ndex3d)
680
681      endif
682
683
684c   -----------------
685c   Moyenne verticale
686c   -----------------
687
688      zamasse(jjb:jje)=0.
689      do l=1,llm
690         zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
691      enddo
692     
693      zavQ(jjb:jje,:,:)=0.
694      do iQ=1,nQ
695         do itr=2,ntr
696            do l=1,llm
697               zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
698     s                             +zvQ(jjb:jje,l,itr,iQ)
699     s                             *zmasse(jjb:jje,l)
700            enddo
701            zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
702            call histwrite(fileid,'a'//znom(itr,iQ),itau,
703     s                     zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
704         enddo
705      enddo
706
707c     on doit pouvoir tracer systematiquement la fonction de courant.
708
709c=====================================================================
710c/////////////////////////////////////////////////////////////////////
711      icum=0                  !///////////////////////////////////////
712      endif ! icum.eq.ncum    !///////////////////////////////////////
713c/////////////////////////////////////////////////////////////////////
714c=====================================================================
715
716      return
717      end
Note: See TracBrowser for help on using the repository browser.