Ignore:
Timestamp:
Jul 18, 2013, 9:13:18 AM (11 years ago)
Author:
Ehouarn Millour
Message:

Improved sponge layer:

  • Sponge tendencies are now computed analytically, instead than using a Forward Euler approximation.
  • Sponge tendencies are added within top_bound, and the sponge is applied after physics tendencies have been taken into account.

These changes imply that GCM results (when using sponge layer) will be differentwrt bench test cases using previous revisions.
EM

Location:
LMDZ5/trunk/libf/dyn3d
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3d/comconst.h

    r1671 r1793  
    66
    77      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
    8      &                 iflag_top_bound
     8     &                 iflag_top_bound,mode_top_bound
    99      COMMON/comconstr/dtvr,daysec,                                     &
    1010     & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
    … …  
    3030      REAL omeg ! (rad/s) rotation rate of the planet
    3131      REAL dissip_factz,dissip_deltaz,dissip_zref
    32       INTEGER iflag_top_bound
    33       REAL tau_top_bound
     32! top_bound sponge:
     33      INTEGER iflag_top_bound ! sponge type
     34      INTEGER mode_top_bound  ! sponge mode
     35      REAL tau_top_bound ! inverse of sponge characteristic time scale (Hz)
    3436      REAL daylen ! length of solar day, in 'standard' day length
    3537      REAL year_day ! Number of standard days in a year

  • LMDZ5/trunk/libf/dyn3d/comvert.h

    r1654 r1793  
    2323      real bps    ! hybrid sigma contribution at mid-layers
    2424      real scaleheight ! atmospheric (reference) scale height (km)
    25       real pseudoalt ! for planets
     25      real pseudoalt ! pseudo-altitude of model levels (km), based on presnivs(),
     26                     ! preff and scaleheight
    2627
    2728      integer disvert_type ! type of vertical discretization:

  • LMDZ5/trunk/libf/dyn3d/conf_gcm.F

    r1697 r1793  
    307307       CALL getin('dissip_zref',dissip_zref )
    308308
     309! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0
     310!                   iflag_top_bound=0 for no sponge
     311!                   iflag_top_bound=1 for sponge over 4 topmost layers
     312!                   iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
    309313       iflag_top_bound=1
     314       CALL getin('iflag_top_bound',iflag_top_bound)
     315
     316! mode_top_bound : fields towards which sponge relaxation will be done:
     317!                  mode_top_bound=0: no relaxation
     318!                  mode_top_bound=1: u and v relax towards 0
     319!                  mode_top_bound=2: u and v relax towards their zonal mean
     320!                  mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
     321       mode_top_bound=3
     322       CALL getin('mode_top_bound',mode_top_bound)
     323
     324! top_bound sponge : inverse of charactericstic relaxation time scale for sponge
    310325       tau_top_bound=1.e-5
    311        CALL getin('iflag_top_bound',iflag_top_bound)
    312326       CALL getin('tau_top_bound',tau_top_bound)
    313327

  • LMDZ5/trunk/libf/dyn3d/leapfrog.F

    r1676 r1793  
    436436     $               clesphy0, dufi,dvfi,dtetafi,dqfi,dpfi  )
    437437
    438          IF (ok_strato) THEN
    439            CALL top_bound( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    440          ENDIF
    441        
    442438c      ajout des tendances physiques:
    443439c      ------------------------------
    … …  
    445441     $                  ucov, vcov, teta , q   ,ps ,
    446442     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     443
     444         IF (ok_strato) THEN
     445           CALL top_bound( vcov,ucov,teta,masse,dtphys)
     446         ENDIF
     447       
    447448c
    448449c  Diagnostique de conservation de l'énergie : difference
    … …  
    476477        ! Sponge layer (if any)
    477478        IF (ok_strato) THEN
    478           dufi(:,:)=0.
    479           dvfi(:,:)=0.
    480           dtetafi(:,:)=0.
    481           dqfi(:,:,:)=0.
    482           dpfi(:)=0.
    483           CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    484           CALL addfi( dtvr, leapf, forward   ,
    485      $                  ucov, vcov, teta , q   ,ps ,
    486      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     479!          dufi(:,:)=0.
     480!          dvfi(:,:)=0.
     481!          dtetafi(:,:)=0.
     482!          dqfi(:,:,:)=0.
     483!          dpfi(:)=0.
     484!          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     485           CALL top_bound( vcov,ucov,teta,masse,dtvr)
     486!          CALL addfi( dtvr, leapf, forward   ,
     487!     $                  ucov, vcov, teta , q   ,ps ,
     488!     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    487489        ENDIF ! of IF (ok_strato)
    488490      ENDIF ! of IF (iflag_phys.EQ.2)

  • LMDZ5/trunk/libf/dyn3d/top_bound.F

    r1279 r1793  
    1       SUBROUTINE top_bound( vcov,ucov,teta,masse, du,dv,dh )
     1!
     2! $Id$
     3!
     4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt)
    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*t))
     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 masse(iip1,jjp1,llm)
    38       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
    3964
    4065c   Local:
    … …  
    4469      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
    4570     
    46       INTEGER NDAMP
    47       PARAMETER (NDAMP=4)
    4871      integer i
    49       REAL,SAVE :: rdamp(llm)
    50 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     72      REAL,SAVE :: rdamp(llm) ! quenching coefficient
     73      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5174
    5275      LOGICAL,SAVE :: first=.true.
    5376
    5477      INTEGER j,l
    55 
    56 
    57 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    5878     
    5979      if (iflag_top_bound.eq.0) return
    … …  
    6181      if (first) then
    6282         if (iflag_top_bound.eq.1) then
    63 ! couche eponge dans les 4 dernieres couches du modele
    64              rdamp(:)=0.
    65              rdamp(llm)=tau_top_bound
    66              rdamp(llm-1)=tau_top_bound/2.
    67              rdamp(llm-2)=tau_top_bound/4.
    68              rdamp(llm-3)=tau_top_bound/8.
     83! sponge quenching over the topmost 4 atmospheric layers
     84             lambda(:)=0.
     85             lambda(llm)=tau_top_bound
     86             lambda(llm-1)=tau_top_bound/2.
     87             lambda(llm-2)=tau_top_bound/4.
     88             lambda(llm-3)=tau_top_bound/8.
    6989         else if (iflag_top_bound.eq.2) then
    70 ! couce eponge dans toutes les couches de pression plus faible que
    71 ! 100 fois la pression de la derniere couche
    72              rdamp(:)=tau_top_bound
     90! sponge quenching over topmost layers down to pressures which are
     91! higher than 100 times the topmost layer pressure
     92             lambda(:)=tau_top_bound
    7393     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    7494         endif
     95
     96! quenching coefficient rdamp(:)
     97!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
     98         rdamp(:)=1.-exp(-lambda(:)*dt)
     99
     100         write(lunout,*)'TOP_BOUND mode',mode_top_bound
     101         write(lunout,*)'Sponge layer coefficients'
     102         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
     103         do l=1,llm
     104           if (rdamp(l).ne.0.) then
     105             write(lunout,'(6(1pe12.4,1x))')
     106     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     107     &           1./lambda(l),lambda(l)
     108           endif
     109         enddo
    75110         first=.false.
    76          print*,'TOP_BOUND rdamp=',rdamp
    77       endif
     111      endif ! of if (first)
    78112
    79113      CALL massbar(masse,massebx,masseby)
    80114
    81       do l=1,llm
     115      ! compute zonal average of vcov and u
     116      if (mode_top_bound.ge.2) then
     117       do l=1,llm
    82118        do j=1,jjm
    83119          vzon(j,l)=0.
    84120          zm=0.
    85121          do i=1,iim
    86 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    87 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    88 ! ne varie qu'en latitude
     122! NB: we can work using vcov zonal mean rather than v since the
     123! cv coefficient (which relates the two) only varies with latitudes
    89124            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    90125            zm=zm+masseby(i,j,l)
    … …  
    92127          vzon(j,l)=vzon(j,l)/zm
    93128        enddo
    94       enddo
     129       enddo
    95130
    96       do l=1,llm
    97         do i=1,iip1
    98           do j=1,jjm
    99             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    100           enddo
    101         enddo
    102       enddo
    103 
    104       do l=1,llm
    105         do j=2,jjm
     131       do l=1,llm
     132        do j=2,jjm ! excluding poles
    106133          uzon(j,l)=0.
    107134          zm=0.
    … …  
    112139          uzon(j,l)=uzon(j,l)/zm
    113140        enddo
    114       enddo
     141       enddo
     142      else ! ucov and vcov will relax towards 0
     143        vzon(:,:)=0.
     144        uzon(:,:)=0.
     145      endif ! of if (mode_top_bound.ge.2)
    115146
    116       do l=1,llm
    117         do j=2,jjm
     147      ! compute zonal average of potential temperature, if necessary
     148      if (mode_top_bound.ge.3) then
     149       do l=1,llm
     150        do j=2,jjm ! excluding poles
    118151          zm=0.
    119152          tzon(j,l)=0.
    … …  
    124157          tzon(j,l)=tzon(j,l)/zm
    125158        enddo
    126       enddo
     159       enddo
     160      endif ! of if (mode_top_bound.ge.3)
    127161
    128 C   AMORTISSEMENTS LINEAIRES:
    129 
    130       do l=1,llm
     162      if (mode_top_bound.ge.1) then
     163       ! Apply sponge quenching on vcov:
     164       do l=1,llm
    131165        do i=1,iip1
    132           do j=2,jjm
    133             du(i,j,l)=du(i,j,l)
    134      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    135             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
     166          do j=1,jjm
     167            vcov(i,j,l)=vcov(i,j,l)
     168     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    136169          enddo
    137170        enddo
    138       enddo
    139      
     171       enddo
    140172
    141       RETURN
     173       ! Apply sponge quenching on ucov:
     174       do l=1,llm
     175        do i=1,iip1
     176          do j=2,jjm ! excluding poles
     177            ucov(i,j,l)=ucov(i,j,l)
     178     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     179          enddo
     180        enddo
     181       enddo
     182      endif ! of if (mode_top_bound.ge.1)
     183
     184      if (mode_top_bound.ge.3) then
     185       ! Apply sponge quenching on teta:
     186       do l=1,llm
     187        do i=1,iip1
     188          do j=2,jjm ! excluding poles
     189            teta(i,j,l)=teta(i,j,l)
     190     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     191          enddo
     192        enddo
     193       enddo
     194      endif ! of if (mode_top_bound.ge.3)
     195   
    142196      END
Note: See TracChangeset for help on using the changeset viewer.