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

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3dmem/top_bound_loc.F

    r1632 r1793  
    1       SUBROUTINE top_bound_loc( vcov,ucov,teta,masse, du,dv,dh )
     1!
     2! $Id: $
     3!
     4      SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt)
    25      USE parallel
    36      IMPLICIT NONE
    … …  
    2528c
    2629c=======================================================================
    27 c-----------------------------------------------------------------------
    28 c   Declarations:
    29 c   -------------
     30
     31! top_bound sponge layer model:
     32! Quenching is modeled as: A(t)=Am+A0*exp(-lambda*t)
     33! where Am is the zonal average of the field (or zero), and lambda the inverse
     34! of the characteristic quenching/relaxation time scale
     35! Thus, assuming Am to be time-independent, field at time t+dt is given by:
     36! A(t+dt)=A(t)-(A(t)-Am)*(1-exp(-lambda*t))
     37! Moreover lambda can be a function of model level (see below), and relaxation
     38! can be toward the average zonal field or just zero (see below).
     39
     40! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.
     41
     42! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst.h)
     43!    iflag_top_bound=0 for no sponge
     44!    iflag_top_bound=1 for sponge over 4 topmost layers
     45!    iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
     46!    mode_top_bound=0: no relaxation
     47!    mode_top_bound=1: u and v relax towards 0
     48!    mode_top_bound=2: u and v relax towards their zonal mean
     49!    mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
     50!    tau_top_bound : inverse of charactericstic relaxation time scale at
     51!                       the topmost layer (Hz)
     52
    3053
    3154#include "comdissipn.h"
     55#include "iniprint.h"
    3256
    3357c   Arguments:
    3458c   ----------
    3559
    36       REAL ucov(iip1,jjb_u:jje_u,llm),vcov(iip1,jjb_v:jje_v,llm)
    37       REAL teta(iip1,jjb_u:jje_u,llm)
    38       REAL masse(iip1,jjb_u:jje_u,llm)
    39       REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)
    40       REAL dh(iip1,jjb_u:jje_u,llm)
     60      real,intent(inout) :: ucov(iip1,jjb_u:jje_u,llm) ! covariant zonal wind
     61      real,intent(inout) :: vcov(iip1,jjb_v:jje_v,llm) ! covariant meridional wind
     62      real,intent(inout) :: teta(iip1,jjb_u:jje_u,llm) ! potential temperature
     63      real,intent(in) :: masse(iip1,jjb_u:jje_u,llm) ! mass of atmosphere
     64      real,intent(in) :: dt ! time step (s) of sponge model
     65
     66!      REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)
     67!      REAL dh(iip1,jjb_u:jje_u,llm)
    4168
    4269c   Local:
    … …  
    4774      REAL tzon(jjb_u:jje_u,llm)
    4875     
    49       INTEGER NDAMP
    50       PARAMETER (NDAMP=4)
    5176      integer i
    5277      REAL,SAVE :: rdamp(llm)
    53 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     78      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5479      LOGICAL,SAVE :: first=.true.
    5580      INTEGER j,l,jjb,jje
    … …  
    5782
    5883      if (iflag_top_bound == 0) return
     84
    5985      if (first) then
    6086c$OMP BARRIER
    6187c$OMP MASTER
    6288         if (iflag_top_bound == 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.
     89! sponge quenching over the topmost 4 atmospheric layers
     90             lambda(:)=0.
     91             lambda(llm)=tau_top_bound
     92             lambda(llm-1)=tau_top_bound/2.
     93             lambda(llm-2)=tau_top_bound/4.
     94             lambda(llm-3)=tau_top_bound/8.
    6995         else if (iflag_top_bound == 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
     96! sponge quenching over topmost layers down to pressures which are
     97! higher than 100 times the topmost layer pressure
     98             lambda(:)=tau_top_bound
    7399     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    74100         endif
     101
     102! quenching coefficient rdamp(:)
     103!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
     104         rdamp(:)=1.-exp(-lambda(:)*dt)
     105
     106         write(lunout,*)'TOP_BOUND mode',mode_top_bound
     107         write(lunout,*)'Sponge layer coefficients'
     108         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
     109         do l=1,llm
     110           if (rdamp(l).ne.0.) then
     111             write(lunout,'(6(1pe12.4,1x))')
     112     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     113     &           1./lambda(l),lambda(l)
     114           endif
     115         enddo
    75116         first=.false.
    76          print*,'TOP_BOUND rdamp=',rdamp
    77117c$OMP END MASTER
    78118c$OMP BARRIER
    79       endif
     119      endif ! of if (first)
    80120
    81121
    82122      CALL massbar_loc(masse,massebx,masseby)
    83 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    84 
    85       jjb=jj_begin
    86       jje=jj_end
    87       IF (pole_sud) jje=jj_end-1
    88 
    89 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    90       do l=1,llm
     123
     124      ! compute zonal average of vcov (or set it to zero)
     125      if (mode_top_bound.ge.2) then
     126       jjb=jj_begin
     127       jje=jj_end
     128       IF (pole_sud) jje=jj_end-1
     129c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     130       do l=1,llm
    91131        do j=jjb,jje
    92132          zm=0.
    93133          vzon(j,l)=0
    94134          do i=1,iim
    95 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    96 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    97 ! ne varie qu'en latitude
     135! NB: we can work using vcov zonal mean rather than v since the
     136! cv coefficient (which relates the two) only varies with latitudes
    98137            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    99138            zm=zm+masseby(i,j,l)
    … …  
    101140          vzon(j,l)=vzon(j,l)/zm
    102141        enddo
    103       enddo
     142       enddo
    104143c$OMP END DO NOWAIT   
    105 
    106 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    107       do l=1,llm
    108         do j=jjb,jje
    109           do i=1,iip1
    110             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    111           enddo
    112         enddo
    113       enddo
    114 c$OMP END DO NOWAIT
    115 
    116       jjb=jj_begin
    117       jje=jj_end
    118       IF (pole_nord) jjb=jj_begin+1
    119       IF (pole_sud)  jje=jj_end-1
    120 
    121 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    122       do l=1,llm
     144      else
     145c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     146       do l=1,llm
     147         vzon(:,l)=0.
     148       enddo
     149c$OMP END DO NOWAIT
     150      endif ! of if (mode_top_bound.ge.2)
     151
     152      ! compute zonal average of u (or set it to zero)
     153      if (mode_top_bound.ge.2) then
     154       jjb=jj_begin
     155       jje=jj_end
     156       IF (pole_nord) jjb=jj_begin+1
     157       IF (pole_sud)  jje=jj_end-1
     158c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     159       do l=1,llm
    123160        do j=jjb,jje
    124161          uzon(j,l)=0.
    … …  
    130167          uzon(j,l)=uzon(j,l)/zm
    131168        enddo
    132       enddo
    133 c$OMP END DO NOWAIT
    134 
     169       enddo
     170c$OMP END DO NOWAIT
     171      else
     172c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     173       do l=1,llm
     174         uzon(:,l)=0.
     175       enddo
     176c$OMP END DO NOWAIT
     177      endif ! of if (mode_top_bound.ge.2)
     178
     179      ! compute zonal average of potential temperature, if necessary
     180      if (mode_top_bound.ge.3) then
     181       jjb=jj_begin
     182       jje=jj_end
     183       IF (pole_nord) jjb=jj_begin+1
     184       IF (pole_sud)  jje=jj_end-1
    135185c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    136       do l=1,llm
     186       do l=1,llm
    137187        do j=jjb,jje
    138188          zm=0.
    … …  
    144194          tzon(j,l)=tzon(j,l)/zm
    145195        enddo
    146       enddo
    147 c$OMP END DO NOWAIT
    148 
    149 C   AMORTISSEMENTS LINEAIRES:
    150 
    151 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    152       do l=1,llm
     196       enddo
     197c$OMP END DO NOWAIT
     198      endif ! of if (mode_top_bound.ge.3)
     199
     200      if (mode_top_bound.ge.1) then
     201       ! Apply sponge quenching on vcov:
     202       jjb=jj_begin
     203       jje=jj_end
     204       IF (pole_sud) jje=jj_end-1
     205
     206c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     207       do l=1,llm
    153208        do j=jjb,jje
    154209          do i=1,iip1
    155             du(i,j,l)=du(i,j,l)
    156      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    157             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
    158           enddo
    159        enddo
    160       enddo
    161 c$OMP END DO NOWAIT
    162      
    163 
    164       RETURN
     210            vcov(i,j,l)=vcov(i,j,l)
     211     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     212          enddo
     213        enddo
     214       enddo
     215c$OMP END DO NOWAIT
     216
     217       ! Apply sponge quenching on ucov:
     218       jjb=jj_begin
     219       jje=jj_end
     220       IF (pole_nord) jjb=jj_begin+1
     221       IF (pole_sud)  jje=jj_end-1
     222
     223c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     224       do l=1,llm
     225        do j=jjb,jje
     226          do i=1,iip1
     227            ucov(i,j,l)=ucov(i,j,l)
     228     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     229          enddo
     230       enddo
     231       enddo
     232c$OMP END DO NOWAIT
     233      endif ! of if (mode_top_bound.ge.1)
     234
     235      if (mode_top_bound.ge.3) then   
     236       ! Apply sponge quenching on teta:
     237       jjb=jj_begin
     238       jje=jj_end
     239       IF (pole_nord) jjb=jj_begin+1
     240       IF (pole_sud)  jje=jj_end-1
     241
     242c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     243       do l=1,llm
     244        do j=jjb,jje
     245          do i=1,iip1
     246            teta(i,j,l)=teta(i,j,l)
     247     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     248          enddo
     249       enddo
     250       enddo
     251c$OMP END DO NOWAIT
     252      endif ! of if (mode_top_bond.ge.3)
     253
    165254      END
Note: See TracChangeset for help on using the changeset viewer.