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/dyn3dpar/top_bound_p.F

    r1279 r1793  
    1       SUBROUTINE top_bound_p( vcov,ucov,teta,masse, du,dv,dh )
     1!
     2! $Id$
     3!
     4      SUBROUTINE top_bound_p(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,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)
     60      real,intent(inout) :: ucov(iip1,jjp1,llm) ! covariant zonal wind
     61      real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
     62      real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
     63      real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere
     64      real,intent(in) :: dt ! time step (s) of sponge model
    3965
    4066c   Local:
    … …  
    4369      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
    4470     
    45       INTEGER NDAMP
    46       PARAMETER (NDAMP=4)
    4771      integer i
    48       REAL,SAVE :: rdamp(llm)
    49 !     &   (/(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)
    5074      LOGICAL,SAVE :: first=.true.
    5175      INTEGER j,l,jjb,jje
    … …  
    5377
    5478      if (iflag_top_bound == 0) return
     79
    5580      if (first) then
    5681c$OMP BARRIER
    5782c$OMP MASTER
    5883         if (iflag_top_bound == 1) then
    59 ! couche eponge dans les 4 dernieres couches du modele
    60              rdamp(:)=0.
    61              rdamp(llm)=tau_top_bound
    62              rdamp(llm-1)=tau_top_bound/2.
    63              rdamp(llm-2)=tau_top_bound/4.
    64              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.
    6590         else if (iflag_top_bound == 2) then
    66 ! couce eponge dans toutes les couches de pression plus faible que
    67 ! 100 fois la pression de la derniere couche
    68              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
    6994     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    7095         endif
     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)'
     104         do l=1,llm
     105           if (rdamp(l).ne.0.) then
     106             write(lunout,'(6(1pe12.4,1x))')
     107     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     108     &           1./lambda(l),lambda(l)
     109           endif
     110         enddo
    71111         first=.false.
    72          print*,'TOP_BOUND rdamp=',rdamp
    73112c$OMP END MASTER
    74113c$OMP BARRIER
    75       endif
     114      endif ! of if (first)
    76115
    77116
    78117      CALL massbar_p(masse,massebx,masseby)
    79 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    80 
    81       jjb=jj_begin
    82       jje=jj_end
    83       IF (pole_sud) jje=jj_end-1
    84 
    85 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    86       do l=1,llm
     118
     119      ! compute zonal average of vcov (or set it to zero)
     120      if (mode_top_bound.ge.2) then
     121       jjb=jj_begin
     122       jje=jj_end
     123       IF (pole_sud) jje=jj_end-1
     124c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     125       do l=1,llm
    87126        do j=jjb,jje
    88127          zm=0.
    89128          vzon(j,l)=0
    90129          do i=1,iim
    91 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    92 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    93 ! ne varie qu'en latitude
     130! NB: we can work using vcov zonal mean rather than v since the
     131! cv coefficient (which relates the two) only varies with latitudes
    94132            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    95133            zm=zm+masseby(i,j,l)
    … …  
    97135          vzon(j,l)=vzon(j,l)/zm
    98136        enddo
    99       enddo
     137       enddo
    100138c$OMP END DO NOWAIT   
    101 
    102 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    103       do l=1,llm
    104         do j=jjb,jje
    105           do i=1,iip1
    106             dv(i,j,l)=dv(i,j,l)-rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    107           enddo
    108         enddo
    109       enddo
    110 c$OMP END DO NOWAIT
    111 
    112       jjb=jj_begin
    113       jje=jj_end
    114       IF (pole_nord) jjb=jj_begin+1
    115       IF (pole_sud)  jje=jj_end-1
    116 
    117 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    118       do l=1,llm
     139      else
     140c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     141       do l=1,llm
     142         vzon(:,l)=0.
     143       enddo
     144c$OMP END DO NOWAIT
     145      endif ! of if (mode_top_bound.ge.2)
     146
     147      ! compute zonal average of u (or set it to zero)
     148      if (mode_top_bound.ge.2) then
     149       jjb=jj_begin
     150       jje=jj_end
     151       IF (pole_nord) jjb=jj_begin+1
     152       IF (pole_sud)  jje=jj_end-1
     153c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     154       do l=1,llm
    119155        do j=jjb,jje
    120156          uzon(j,l)=0.
    … …  
    126162          uzon(j,l)=uzon(j,l)/zm
    127163        enddo
    128       enddo
    129 c$OMP END DO NOWAIT
    130 
     164       enddo
     165c$OMP END DO NOWAIT
     166      else
     167c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     168       do l=1,llm
     169         uzon(:,l)=0.
     170       enddo
     171c$OMP END DO NOWAIT
     172      endif ! of if (mode_top_bound.ge.2)
     173
     174      ! compute zonal average of potential temperature, if necessary
     175      if (mode_top_bound.ge.3) then
     176       jjb=jj_begin
     177       jje=jj_end
     178       IF (pole_nord) jjb=jj_begin+1
     179       IF (pole_sud)  jje=jj_end-1
    131180c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    132       do l=1,llm
     181       do l=1,llm
    133182        do j=jjb,jje
    134183          zm=0.
    … …  
    140189          tzon(j,l)=tzon(j,l)/zm
    141190        enddo
    142       enddo
    143 c$OMP END DO NOWAIT
    144 
    145 C   AMORTISSEMENTS LINEAIRES:
    146 
    147 c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    148       do l=1,llm
     191       enddo
     192c$OMP END DO NOWAIT
     193      endif ! of if (mode_top_bound.ge.3)
     194
     195      if (mode_top_bound.ge.1) then
     196       ! Apply sponge quenching on vcov:
     197       jjb=jj_begin
     198       jje=jj_end
     199       IF (pole_sud) jje=jj_end-1
     200
     201c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     202       do l=1,llm
    149203        do j=jjb,jje
    150204          do i=1,iip1
    151             du(i,j,l)=du(i,j,l)
    152      s               -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    153             dh(i,j,l)=dh(i,j,l)-rdamp(l)*(teta(i,j,l)-tzon(j,l))
    154           enddo
    155        enddo
    156       enddo
    157 c$OMP END DO NOWAIT
    158      
    159 
    160       RETURN
     205            vcov(i,j,l)=vcov(i,j,l)
     206     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     207          enddo
     208        enddo
     209       enddo
     210c$OMP END DO NOWAIT
     211
     212       ! Apply sponge quenching on ucov:
     213       jjb=jj_begin
     214       jje=jj_end
     215       IF (pole_nord) jjb=jj_begin+1
     216       IF (pole_sud)  jje=jj_end-1
     217
     218c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     219       do l=1,llm
     220        do j=jjb,jje
     221          do i=1,iip1
     222            ucov(i,j,l)=ucov(i,j,l)
     223     &                  -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     224          enddo
     225       enddo
     226       enddo
     227c$OMP END DO NOWAIT
     228      endif ! of if (mode_top_bound.ge.1)
     229
     230      if (mode_top_bound.ge.3) then   
     231       ! Apply sponge quenching on teta:
     232       jjb=jj_begin
     233       jje=jj_end
     234       IF (pole_nord) jjb=jj_begin+1
     235       IF (pole_sud)  jje=jj_end-1
     236
     237c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     238       do l=1,llm
     239        do j=jjb,jje
     240          do i=1,iip1
     241            teta(i,j,l)=teta(i,j,l)
     242     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     243          enddo
     244       enddo
     245       enddo
     246c$OMP END DO NOWAIT
     247      endif ! of if (mode_top_bond.ge.3)
     248
    161249      END
Note: See TracChangeset for help on using the changeset viewer.