/[lmdze]/trunk/libf/dyn3d/bilan_dyn.f90
ViewVC logotype

Diff of /trunk/libf/dyn3d/bilan_dyn.f90

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.39  
changed lines
  Added in v.40

  ViewVC Help
Powered by ViewVC 1.1.21