source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/dyn/sponge_mod_p.f90 @ 224

module sponge_mod_p

implicit none


! sponge parameters (set/read via conf_gcm.F)
logical,save :: callsponge  ! do we use a sponge on upper layers
integer,save :: mode_sponge ! sponge mode
integer,save :: nsponge ! number of sponge layers 
real,save :: tetasponge  ! sponge time scale (s) at topmost layer


contains

      subroutine sponge_p(ucov,vcov,h,ps,dt,mode)

! Sponge routine: Quench ucov, vcov and potential temperature near the
!                 top of the model
! Depending on 'mode' relaxation of variables is towards:
! mode = 0 : h -> h_mean , ucov -> 0 , vcov -> 0
! mode = 1 : h -> h_mean , ucov -> ucov_mean , vcov -> 0
! mode >= 2 : h -> h_mean , ucov -> ucov_mean , vcov -> vcov_mean
! Number of layer over which sponge is applied is 'nsponge' (read from def file)
! Time scale for quenching at top level is given by 'tetasponge' (read from
! def file) and doubles as level indexes decrease.
! Quenching is modeled as: A(t)=Am+A0exp(-lambda*t)
! where Am is the zonal average of the field (or zero), and lambda the inverse
! of the characteristic quenching/relaxation time scale
! Thus, assuming Am to be time-independent, field at time t+dt is given by:
! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*dt))

      USE Write_Field_p
      use parallel_lmdz, only: pole_sud,pole_nord,jj_begin,jj_end
      implicit none
!-----------------------------------------------------------------------
!   INCLUDE 'dimensions.h'
!
!   dimensions.h contient les dimensions du modele
!   ndm est tel que iim=2**ndm
!-----------------------------------------------------------------------

      INTEGER iim,jjm,llm,ndm

      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)

!-----------------------------------------------------------------------
!
! $Header$
!
!
!  ATTENTION!!!!: ce fichier include est compatible format fixe/format libre
!                 veillez  n'utiliser que des ! pour les commentaires
!                 et  bien positionner les & des lignes de continuation
!                 (les placer en colonne 6 et en colonne 73)
!
!
!-----------------------------------------------------------------------
!   INCLUDE 'paramet.h'

      INTEGER  iip1,iip2,iip3,jjp1,llmp1,llmp2,llmm1
      INTEGER  kftd,ip1jm,ip1jmp1,ip1jmi1,ijp1llm
      INTEGER  ijmllm,mvar
      INTEGER jcfil,jcfllm

      PARAMETER( iip1= iim+1,iip2=iim+2,iip3=iim+3                       &
     &    ,jjp1=jjm+1-1/jjm)
      PARAMETER( llmp1 = llm+1,  llmp2 = llm+2, llmm1 = llm-1 )
      PARAMETER( kftd  = iim/2 -ndm )
      PARAMETER( ip1jm  = iip1*jjm,  ip1jmp1= iip1*jjp1 )
      PARAMETER( ip1jmi1= ip1jm - iip1 )
      PARAMETER( ijp1llm= ip1jmp1 * llm, ijmllm= ip1jm * llm )
      PARAMETER( mvar= ip1jmp1*( 2*llm+1) + ijmllm )
      PARAMETER( jcfil=jjm/2+5, jcfllm=jcfil*llm )

!-----------------------------------------------------------------------
!
! $Header$
!
!-----------------------------------------------------------------------
! INCLUDE comdissip.h

      COMMON/comdissip/                                                 &
     &    niterdis,coefdis,tetavel,tetatemp,gamdissip


      INTEGER niterdis

      REAL tetavel,tetatemp,coefdis,gamdissip

!-----------------------------------------------------------------------
!
! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
!
!-----------------------------------------------------------------------
!   INCLUDE 'comvert.h'

      COMMON/comvertr/ap(llm+1),bp(llm+1),presnivs(llm),dpres(llm),     &
     &               pa,preff,nivsigs(llm),nivsig(llm+1),               &
     &               aps(llm),bps(llm),scaleheight,pseudoalt(llm)

      common/comverti/disvert_type, pressure_exner

      real ap     ! hybrid pressure contribution at interlayers
      real bp     ! hybrid sigma contribution at interlayer
      real presnivs ! (reference) pressure at mid-layers
      real dpres
      real pa     ! reference pressure (Pa) at which hybrid coordinates
                  ! become purely pressure
      real preff  ! reference surface pressure (Pa)
      real nivsigs
      real nivsig
      real aps    ! hybrid pressure contribution at mid-layers
      real bps    ! hybrid sigma contribution at mid-layers
      real scaleheight ! atmospheric (reference) scale height (km)
      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
                     ! preff and scaleheight

      integer disvert_type ! type of vertical discretization:
                           ! 1: Earth (default for planet_type==earth),
                           !     automatic generation
                           ! 2: Planets (default for planet_type!=earth),
                           !     using 'z2sig.def' (or 'esasig.def) file

      logical pressure_exner
!     compute pressure inside layers using Exner function, else use mean
!     of pressure values at interfaces

 !-----------------------------------------------------------------------
!
! $Header$
!
!CDK comgeom2
      COMMON/comgeom/                                                   &
     & cu(iip1,jjp1),cv(iip1,jjm),unscu2(iip1,jjp1),unscv2(iip1,jjm)  , &
     & aire(iip1,jjp1),airesurg(iip1,jjp1),aireu(iip1,jjp1)           , &
     & airev(iip1,jjm),unsaire(iip1,jjp1),apoln,apols                 , &
     & unsairez(iip1,jjm),airuscv2(iip1,jjm),airvscu2(iip1,jjm)       , &
     & aireij1(iip1,jjp1),aireij2(iip1,jjp1),aireij3(iip1,jjp1)       , &
     & aireij4(iip1,jjp1),alpha1(iip1,jjp1),alpha2(iip1,jjp1)         , &
     & alpha3(iip1,jjp1),alpha4(iip1,jjp1),alpha1p2(iip1,jjp1)        , &
     & alpha1p4(iip1,jjp1),alpha2p3(iip1,jjp1),alpha3p4(iip1,jjp1)    , &
     & fext(iip1,jjm),constang(iip1,jjp1), rlatu(jjp1),rlatv(jjm),      &
     & rlonu(iip1),rlonv(iip1),cuvsurcv(iip1,jjm),cvsurcuv(iip1,jjm)  , &
     & cvusurcu(iip1,jjp1),cusurcvu(iip1,jjp1)                        , &
     & cuvscvgam1(iip1,jjm),cuvscvgam2(iip1,jjm),cvuscugam1(iip1,jjp1), &
     & cvuscugam2(iip1,jjp1),cvscuvgam(iip1,jjm),cuscvugam(iip1,jjp1) , &
     & unsapolnga1,unsapolnga2,unsapolsga1,unsapolsga2                , &
     & unsair_gam1(iip1,jjp1),unsair_gam2(iip1,jjp1)                  , &
     & unsairz_gam(iip1,jjm),aivscu2gam(iip1,jjm),aiuscv2gam(iip1,jjm)  &
     & , xprimu(iip1),xprimv(iip1)


      REAL                                                               &
     & cu,cv,unscu2,unscv2,aire,airesurg,aireu,airev,apoln,apols,unsaire &
     & ,unsairez,airuscv2,airvscu2,aireij1,aireij2,aireij3,aireij4     , &
     & alpha1,alpha2,alpha3,alpha4,alpha1p2,alpha1p4,alpha2p3,alpha3p4 , &
     & fext,constang,rlatu,rlatv,rlonu,rlonv,cuvscvgam1,cuvscvgam2     , &
     & cvuscugam1,cvuscugam2,cvscuvgam,cuscvugam,unsapolnga1           , &
     & unsapolnga2,unsapolsga1,unsapolsga2,unsair_gam1,unsair_gam2     , &
     & unsairz_gam,aivscu2gam,aiuscv2gam,cuvsurcv,cvsurcuv,cvusurcu    , &
     & cusurcvu,xprimu,xprimv
!
! $Header$
!
!
! gestion des impressions de sorties et de débogage
! lunout:    unité du fichier dans lequel se font les sorties 
!                           (par defaut 6, la sortie standard)
! prt_level: niveau d'impression souhaité (0 = minimum)
!
      INTEGER lunout, prt_level
      COMMON /comprint/ lunout, prt_level

! Arguments:
!------------
      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
      real,intent(inout) :: h(iip1,jjp1,llm) ! potential temperature
!      real,intent(in) :: pext(iip1,jjp1) ! extensive pressure
      real,intent(in) :: ps(iip1,jjp1) ! surface pressure (Pa)
      real,intent(in) :: dt   ! time step
      integer,intent(in) :: mode  ! sponge mode

!   Local:
!   ------

      real,save :: sig_s(llm) !sigma au milieu des couches
      REAL vm,um,hm,ptot(jjp1)
      real,save :: cst(llm)
      real :: pext(iip1,jjp1) ! extensive pressure

      INTEGER l,i,j
      integer :: jjb,jje
      integer,save :: l0 ! layer down to which sponge is applied

      real ssum

      real zkm
      logical,save :: firstcall=.true.



      if (firstcall) then
!$OMP BARRIER
!$OMP MASTER
       ! build approximative sigma levels at midlayer
        do l=1,llm
          sig_s(l)=((ap(l)+ap(l+1))/preff+bp(l)+bp(l+1))/2.
        enddo

        l0=llm-nsponge+1
 
        write(lunout,*)
        write(lunout,*)'sponge mode',mode
        write(lunout,*)'nsponge tetasponge ',nsponge,tetasponge
        write(lunout,*)'Coeffs for the sponge layer'
        write(lunout,*)'Z (km)         tau             cst'
        do l=llm,l0,-1
          ! double time scale with every level, starting from the top
          cst(l)=1.-exp(-dt/(tetasponge*2**(llm-l)))
        enddo

        do l=l0,llm
           zkm=-scaleheight*log(sig_s(l))
           print*,zkm,tetasponge*2**(llm-l),cst(l)
        enddo
        PRINT*

        firstcall=.false.
!$OMP END MASTER
!$OMP BARRIER
      endif ! of if (firstcall)

!-----------------------------------------------------------------------
!   compute sponge relaxation:
!   -------------------------
      jjb=jj_begin
      jje=jj_end+1 ! +1 because vcov(j) computations require pext(j+1) & ptot(j+1)
      IF (pole_sud)  jje=jj_end-1+1
!$OMP MASTER
      pext(1:iip1,jjb:jje)=ps(1:iip1,jjb:jje)*aire(1:iip1,jjb:jje)
      do j=jjb,jje
         ptot(j)=sum(pext(1:iim,j))
      enddo
!$OMP END MASTER
!$OMP BARRIER


!   potential temperature
      jjb=jj_begin
      jje=jj_end
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
      do l=l0,llm
         do j=jjb,jje
            ! compute zonal average
            hm=0.
            do i=1,iim
               hm=hm+h(i,j,l)*pext(i,j)
            enddo
            hm=hm/ptot(j)
            ! update h()
            do i=1,iim
               h(i,j,l)=h(i,j,l)-cst(l)*(h(i,j,l)-hm)
            enddo
            h(iip1,j,l)=h(1,j,l)
         enddo
      enddo
!$OMP END DO NOWAIT

!   zonal wind
      jjb=jj_begin
      jje=jj_end
      IF (pole_nord) jjb=jj_begin+1
      IF (pole_sud)  jje=jj_end-1
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
      do l=l0,llm
         do j=jjb,jje
            um=0.
            if(mode.ge.1) then ! compute zonal average
               do i=1,iim
                  um=um+0.5*ucov(i,j,l)*(pext(i,j)+pext(i+1,j))/cu(i,j)
               enddo
               um=um/ptot(j)
            endif
            ! update ucov()
            do i=1,iim
               ucov(i,j,l)=ucov(i,j,l)-cst(l)*(ucov(i,j,l)-um*cu(i,j))
            enddo
            ucov(iip1,j,l)=ucov(1,j,l)
         enddo
      enddo
!$OMP END DO NOWAIT

!  meridional wind
      jjb=jj_begin
      jje=jj_end
      IF (pole_sud) jje=jj_end-1
!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
      do l=l0,llm
         do j=jjb,jje
            vm=0.
            if(mode.ge.2) then ! compute zonal average
               do i=1,iim
                  vm=vm+vcov(i,j,l)*(pext(i,j)+pext(i,j+1))/cv(i,j)
               enddo
               vm=vm/(ptot(j)+ptot(j+1))
            endif
            ! update vcov
            do i=1,iim
               vcov(i,j,l)=vcov(i,j,l)-cst(l)*(vcov(i,j,l)-vm*cv(i,j))
            enddo
            vcov(iip1,j,l)=vcov(1,j,l)
         enddo
      enddo
!$OMP END DO

      end subroutine sponge_p
      
end module sponge_mod_p