source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/radsponge.F90 @ 310

subroutine radsponge(ngrid,nlayer,pplay,pplev,temp,pdt,dtphys,dtradsponge)
! add a (vertical) sponge effect on temperature to damp oscillations
! near the model top (topmost 
! 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 ioipsl_getincom_p, only: getin_p
use mod_phys_lmdz_para, only : is_master

implicit none

integer,intent(in) :: ngrid ! number of atmospheric columns
integer,intent(in) :: nlayer ! number of atmospheric layers
real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
real,intent(in) :: pplev(ngrid,nlayer+1) ! interlayer pressure (Pa)
real,intent(in) :: temp(ngrid,nlayer) ! temperature (K)
real,intent(in) :: pdt(ngrid,nlayer) ! tendecies on temperature from previous processes (K/s)
real,intent(in) :: dtphys ! physics time step (s)
real,intent(out) :: dtradsponge(ngrid,nlayer) ! tendencies on temperature from sponge

! local variables:
logical,save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
integer,save :: iflag_sponge   !0 --> for no sponge
                               !1 --> for sponge over 4 topmost layers
                               !2 --> for sponge from top to ~1% of top layer pressure
                               !3 --> for sponge over last nb_sponge_layers
real,save :: tau_sponge !inverse of charactericstic relaxation time scale at the topmost layer (Hz)
real,save :: nb_sponge_layers ! number of layers over which the sponge extends
!$OMP THREADPRIVATE(iflag_sponge,tau_sponge,nb_sponge_layers)
integer,save :: nlayer0 ! layer down to which sponge is applied
real,allocatable,save :: rdamp(:) ! quenching coefficient
real,allocatable,save :: lambda(:) ! inverse or quenching time scale (Hz)
real,allocatable,save :: mass(:) ! g*mass (per m2) of atmospheric layer
real,save :: spongemass ! (pseudo-) mass of atmosphere over the entire sponge layer

real :: temp_cur(ngrid,nlayer) ! temperature (K)
real :: mean_temp(ngrid) ! mean temperature towards which to relax

integer :: i,k

if (firstcall) then
  ! load sponge parameters
  iflag_sponge=0 ! default: no sponge
  call getin_p("iflag_sponge",iflag_sponge)
  tau_sponge=1.E-4 ! default value
  call getin_p("tau_sponge",tau_sponge)
  if (iflag_sponge==3) then
    call getin_p("nb_sponge_layers",nb_sponge_layers)
  endif

!$OMP MASTER
  allocate(rdamp(nlayer))
  allocate(lambda(nlayer))
  allocate(mass(nlayer))
  
  if (iflag_sponge == 1) then
! sponge quenching over the topmost 4 atmospheric layers
    nb_sponge_layers=4
    nlayer0=nlayer-3
    lambda(:)=0.
    lambda(nlayer)=tau_sponge
    lambda(nlayer-1)=tau_sponge/2.
    lambda(nlayer-2)=tau_sponge/4.
    lambda(nlayer-3)=tau_sponge/8.
  else if (iflag_sponge == 2) then
! sponge quenching over topmost layers down to pressures which are
! higher than 100 times the topmost layer pressure
    do k=nlayer,1,-1
      lambda(k)=tau_sponge*max(pplay(1,nlayer)/pplay(1,k)-0.01,0.)
      if (lambda(k).ne.0.) nlayer0=k
    enddo
  else if (iflag_sponge == 3) then
    nlayer0=nlayer-nb_sponge_layers+1
    lambda(:)=0
    do k=nlayer,nlayer0,-1
      lambda(k)=tau_sponge/(2.**(nlayer-k))
    enddo
  else
    write(*,*) "radsponge: bad value for iflag_sponge:",iflag_sponge
    stop
  endif

  rdamp(:)=1.-exp(-lambda(:)*dtphys)
  
  ! compute (pseudo-) atmospheric mass of each layer
  ! (here we assume that they won't vay with time,
  ! which is always the case near model top)
  do k=1,nlayer
    mass(k)=pplev(1,k)-pplev(1,k+1)
  enddo
  spongemass=0
  do k=nlayer0,nlayer
    spongemass=spongemass+mass(k)
  enddo
  
  if (is_master) then
    write(*,*) "radsponge: nlayer0=",nlayer0," spongemass=",spongemass
    write(*,*) " layer    lambda       rdamp        mass"
    do k=nlayer0,nlayer
      write(*,'(i3,3x,3(1pe14.7,1x))') k,lambda(k),rdamp(k),mass(k)
    enddo
  endif
!$OMP END MASTER
!$OMP BARRIER
  
  firstcall=.false.
endif ! of if (firstcall)


! 1. compute current temperature and average teperature toward which to relax

temp_cur(:,:)=temp(:,:)+dtphys*pdt(:,:)

! compute mean temperature as a mass-weighed average
mean_temp(1:ngrid)=0
do k=nlayer0,nlayer
  mean_temp(1:ngrid)=mean_temp(1:ngrid)+temp_cur(1:ngrid,k)*mass(k)
enddo
mean_temp(1:ngrid)=mean_temp(1:ngrid)/spongemass

! 2. Temperature tendency due to sponge
dtradsponge(:,:)=0
do k=nlayer0,nlayer
  ! A(t+dt)-A(t)=-(A(t)-Am)*(1-exp(-lambda*dt))
  dtradsponge(1:ngrid,k)=-(temp_cur(1:ngrid,k)-mean_temp(1:ngrid))*rdamp(k)
enddo
! Ehouarn debug:
!write(*,*)"radsponge: mean_temp(:)=",mean_temp(:)
!write(*,*)" layer    temp     temp_cur     dtradsponge"
!do k=nlayer0,nlayer
!  write(*,'(i3,3x,3(1pe14.7,1x))') k,temp(1,k),temp_cur(1,k),dtradsponge(1,k)
!enddo

! return dtradsponge as a tendency in  K/s
dtradsponge(1:ngrid,nlayer0:nlayer)=dtradsponge(1:ngrid,nlayer0:nlayer)/dtphys

end subroutine radsponge