source: trunk/libf/dyn3d/top_bound.F @ 128

      SUBROUTINE top_bound( vcov,ucov,teta,phi,masse, du,dv,dh )
      IMPLICIT NONE
c
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comvert.h"
#include "comgeom2.h"


c ..  DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
C     F. LOTT DEC. 2006
c                                 (  10/12/06  )

c=======================================================================
c
c   Auteur:  F. LOTT  
c   -------
c
c   Objet:
c   ------
c
c   Dissipation linéaire (ex top_bound de la physique)
c
c=======================================================================
c-----------------------------------------------------------------------
c   Declarations:
c   -------------

! #include "comgeom.h"
#include "comdissipn.h"

c   Arguments:
c   ----------

      REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
      REAL phi(iip1,jjp1,llm)                  ! geopotentiel
      REAL masse(iip1,jjp1,llm)
      REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)

c   Local:
c   ------

      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
      
      INTEGER NDAMP
      PARAMETER (NDAMP=4)
      integer i
      REAL,SAVE :: rdamp(llm)
!     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/) 

      LOGICAL,SAVE :: first=.true.

      REAL zkm
      INTEGER j,l


C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
      
      if (iflag_top_bound.eq.0) return

      if (first) then
         if (iflag_top_bound.eq.1) then
! couche eponge dans les 4 dernieres couches du modele
             rdamp(:)=0.
             rdamp(llm)=tau_top_bound
             rdamp(llm-1)=tau_top_bound/2.
             rdamp(llm-2)=tau_top_bound/4.
             rdamp(llm-3)=tau_top_bound/8.
         else if (iflag_top_bound.eq.2) then
! couche eponge dans toutes les couches de pression plus faible que
! 100 fois la pression de la derniere couche
             rdamp(:)=tau_top_bound
     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
         endif
         first=.false.
         print*,'TOP_BOUND mode',mode_top_bound
         print*,'Coeffs pour la couche eponge a l equateur'
         print*,'p (Pa)  z(km)  tau (s)'
         do l=1,llm
           if (rdamp(l).ne.0.) then
            zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
          print*,presnivs(l),zkm,1./rdamp(l)
           endif
         enddo
      endif

      CALL massbar(masse,massebx,masseby)

c   mode = 0 : pas de sponge
c   mode = 1 : u et v -> 0
c   mode = 2 : u et v -> moyenne zonale
c   mode = 3 : u, v et h -> moyenne zonale

      if (mode_top_bound.ge.2) then
       do l=1,llm
        do j=1,jjm
          vzon(j,l)=0.
          zm=0.
          do i=1,iim
! Rm: on peut travailler directement avec la moyenne zonale de vcov
! plutot qu'avec celle de v car le coefficient cv qui relie les deux
! ne varie qu'en latitude
            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
            zm=zm+masseby(i,j,l)
          enddo
          vzon(j,l)=vzon(j,l)/zm
        enddo
       enddo

       do l=1,llm
        do j=2,jjm
          uzon(j,l)=0.
          zm=0.
          do i=1,iim
            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
            zm=zm+massebx(i,j,l)
          enddo
          uzon(j,l)=uzon(j,l)/zm
        enddo
       enddo
      else
          vzon(:,:)=0.
          uzon(:,:)=0.
      endif

      if (mode_top_bound.ge.3) then
       do l=1,llm
        do j=2,jjm
          zm=0.
          tzon(j,l)=0.
          do i=1,iim
            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
            zm=zm+masse(i,j,l)
          enddo
          tzon(j,l)=tzon(j,l)/zm
        enddo
       enddo
      endif

C   AMORTISSEMENTS LINEAIRES:

      if (mode_top_bound.ge.1) then
       do l=1,llm
        do i=1,iip1
          do j=1,jjm
            dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
          enddo
        enddo
       enddo

       do l=1,llm
        do i=1,iip1
          do j=2,jjm
            du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
          enddo
        enddo
       enddo
      endif

      if (mode_top_bound.ge.3) then
       do l=1,llm
        do i=1,iip1
          do j=2,jjm
            dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
          enddo
        enddo
       enddo
      endif
      
      RETURN
      END