Ignore:
Timestamp:
Jul 24, 2013, 1:36:44 PM (11 years ago)
Author:
emillour
Message:

Common dynamics: Improved sponge layer scheme (top_bound):

  • Sponge tendencies are now computed analytically, instead of using a Forward Euler approximation.
  • Sponge tendencies are now added within top_bound.

EM

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.COMMON/libf/dyn3d/top_bound.F

    r110 r1010  
    1       SUBROUTINE top_bound( vcov,ucov,teta,phi,masse, du,dv,dh )
     1!
     2! $Id: top_bound.F 1793 2013-07-18 07:13:18Z emillour $
     3!
     4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt,ducov)
    25      IMPLICIT NONE
    36c
    … …  
    2427c
    2528c=======================================================================
    26 c-----------------------------------------------------------------------
    27 c   Declarations:
    28 c   -------------
    2929
    30 ! #include "comgeom.h"
     30! top_bound sponge layer model:
     31! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t)
     32! where Am is the zonal average of the field (or zero), and lambda the inverse
     33! of the characteristic quenching/relaxation time scale
     34! Thus, assuming Am to be time-independent, field at time t+dt is given by:
     35! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*dt))
     36! Moreover lambda can be a function of model level (see below), and relaxation
     37! can be toward the average zonal field or just zero (see below).
     38
     39! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.
     40
     41! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst.h)
     42!    iflag_top_bound=0 for no sponge
     43!    iflag_top_bound=1 for sponge over 4 topmost layers
     44!    iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
     45!    mode_top_bound=0: no relaxation
     46!    mode_top_bound=1: u and v relax towards 0
     47!    mode_top_bound=2: u and v relax towards their zonal mean
     48!    mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
     49!    tau_top_bound : inverse of charactericstic relaxation time scale at
     50!                       the topmost layer (Hz)
     51
     52
    3153#include "comdissipn.h"
     54#include "iniprint.h"
    3255
    3356c   Arguments:
    3457c   ----------
    3558
    36       REAL ucov(iip1,jjp1,llm),vcov(iip1,jjm,llm),teta(iip1,jjp1,llm)
    37       REAL phi(iip1,jjp1,llm)                  ! geopotentiel
    38       REAL masse(iip1,jjp1,llm)
    39       REAL dv(iip1,jjm,llm),du(iip1,jjp1,llm),dh(iip1,jjp1,llm)
     59      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
     60      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
     61      real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
     62      real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere
     63      real,intent(in) :: dt ! time step (s) of sponge model
     64      real,intent(out) :: ducov(iip1,jjp1,llm) ! increment on ucov due to sponge
    4065
    4166c   Local:
    … …  
    4570      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
    4671     
    47       INTEGER NDAMP
    48       PARAMETER (NDAMP=4)
    4972      integer i
    50       REAL,SAVE :: rdamp(llm)
    51 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     73      REAL,SAVE :: rdamp(llm) ! quenching coefficient
     74      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5275
    5376      LOGICAL,SAVE :: first=.true.
    5477
    55       REAL zkm
    5678      INTEGER j,l
    57 
    58 
    59 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    6079     
    6180      if (iflag_top_bound.eq.0) return
    … …  
    6382      if (first) then
    6483         if (iflag_top_bound.eq.1) then
    65 ! couche eponge dans les 4 dernieres couches du modele
    66              rdamp(:)=0.
    67              rdamp(llm)=tau_top_bound
    68              rdamp(llm-1)=tau_top_bound/2.
    69              rdamp(llm-2)=tau_top_bound/4.
    70              rdamp(llm-3)=tau_top_bound/8.
     84! sponge quenching over the topmost 4 atmospheric layers
     85             lambda(:)=0.
     86             lambda(llm)=tau_top_bound
     87             lambda(llm-1)=tau_top_bound/2.
     88             lambda(llm-2)=tau_top_bound/4.
     89             lambda(llm-3)=tau_top_bound/8.
    7190         else if (iflag_top_bound.eq.2) then
    72 ! couche eponge dans toutes les couches de pression plus faible que
    73 ! 100 fois la pression de la derniere couche
    74              rdamp(:)=tau_top_bound
     91! sponge quenching over topmost layers down to pressures which are
     92! higher than 100 times the topmost layer pressure
     93             lambda(:)=tau_top_bound
    7594     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    7695         endif
    77          first=.false.
    78          print*,'TOP_BOUND mode',mode_top_bound
    79          print*,'Coeffs pour la couche eponge a l equateur'
    80          print*,'p (Pa)  z(km)  tau (s)'
     96
     97! quenching coefficient rdamp(:)
     98!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
     99         rdamp(:)=1.-exp(-lambda(:)*dt)
     100
     101         write(lunout,*)'TOP_BOUND mode',mode_top_bound
     102         write(lunout,*)'Sponge layer coefficients'
     103         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
    81104         do l=1,llm
    82105           if (rdamp(l).ne.0.) then
    83             zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
    84           print*,presnivs(l),zkm,1./rdamp(l)
     106             write(lunout,'(6(1pe12.4,1x))')
     107     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     108     &           1./lambda(l),lambda(l)
    85109           endif
    86110         enddo
    87       endif
     111         first=.false.
     112      endif ! of if (first)
    88113
    89114      CALL massbar(masse,massebx,masseby)
    90115
    91 c   mode = 0 : pas de sponge
    92 c   mode = 1 : u et v -> 0
    93 c   mode = 2 : u et v -> moyenne zonale
    94 c   mode = 3 : u, v et h -> moyenne zonale
    95 
     116      ! compute zonal average of vcov and u
    96117      if (mode_top_bound.ge.2) then
    97118       do l=1,llm
    … …  
    100121          zm=0.
    101122          do i=1,iim
    102 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    103 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    104 ! ne varie qu'en latitude
     123! NB: we can work using vcov zonal mean rather than v since the
     124! cv coefficient (which relates the two) only varies with latitudes
    105125            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    106126            zm=zm+masseby(i,j,l)
    … …  
    111131
    112132       do l=1,llm
    113         do j=2,jjm
     133        do j=2,jjm ! excluding poles
    114134          uzon(j,l)=0.
    115135          zm=0.
    … …  
    121141        enddo
    122142       enddo
    123       else
    124           vzon(:,:)=0.
    125           uzon(:,:)=0.
    126       endif
     143      else ! ucov and vcov will relax towards 0
     144        vzon(:,:)=0.
     145        uzon(:,:)=0.
     146      endif ! of if (mode_top_bound.ge.2)
    127147
     148      ! compute zonal average of potential temperature, if necessary
    128149      if (mode_top_bound.ge.3) then
    129150       do l=1,llm
    130         do j=2,jjm
     151        do j=2,jjm ! excluding poles
    131152          zm=0.
    132153          tzon(j,l)=0.
    … …  
    138159        enddo
    139160       enddo
    140       endif
    141 
    142 C   AMORTISSEMENTS LINEAIRES:
     161      endif ! of if (mode_top_bound.ge.3)
    143162
    144163      if (mode_top_bound.ge.1) then
     164       ! Apply sponge quenching on vcov:
    145165       do l=1,llm
    146166        do i=1,iip1
    147167          do j=1,jjm
    148             dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     168            vcov(i,j,l)=vcov(i,j,l)
     169     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    149170          enddo
    150171        enddo
    151172       enddo
    152173
     174       ! Apply sponge quenching on ucov:
    153175       do l=1,llm
    154176        do i=1,iip1
    155           do j=2,jjm
    156             du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     177          do j=2,jjm ! excluding poles
     178            ducov(i,j,l)=-rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     179            ucov(i,j,l)=ucov(i,j,l)
     180     &                  +ducov(i,j,l)
    157181          enddo
    158182        enddo
    159183       enddo
    160       endif
     184      endif ! of if (mode_top_bound.ge.1)
    161185
    162186      if (mode_top_bound.ge.3) then
     187       ! Apply sponge quenching on teta:
    163188       do l=1,llm
    164189        do i=1,iip1
    165           do j=2,jjm
    166             dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     190          do j=2,jjm ! excluding poles
     191            teta(i,j,l)=teta(i,j,l)
     192     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
    167193          enddo
    168194        enddo
    169195       enddo
    170       endif
    171      
    172       RETURN
     196      endif ! of if (mode_top_bound.ge.3)
     197   
    173198      END
Note: See TracChangeset for help on using the changeset viewer.