source: LMDZ5/trunk/libf/dyn3dmem/bilan_dyn_p.F @ 1632

Last change on this file since 1632 was 1632, checked in by Laurent Fairhead, 13 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

File size: 20.0 KB
Line 
1!
2! $Id: bilan_dyn_p.F 1299 2010-01-20 14:27:21Z fairhead $
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./REAL(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.