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/dyn3dpar
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3dpar/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/dyn3dpar/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/dyn3dpar/conf_gcm.F

    r1734 r1793  
    334334       CALL getin('dissip_zref',dissip_zref )
    335335
     336! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0
     337!                   iflag_top_bound=0 for no sponge
     338!                   iflag_top_bound=1 for sponge over 4 topmost layers
     339!                   iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
    336340       iflag_top_bound=1
     341       CALL getin('iflag_top_bound',iflag_top_bound)
     342
     343! mode_top_bound : fields towards which sponge relaxation will be done:
     344!                  mode_top_bound=0: no relaxation
     345!                  mode_top_bound=1: u and v relax towards 0
     346!                  mode_top_bound=2: u and v relax towards their zonal mean
     347!                  mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
     348       mode_top_bound=3
     349       CALL getin('mode_top_bound',mode_top_bound)
     350
     351! top_bound sponge : inverse of charactericstic relaxation time scale for sponge
    337352       tau_top_bound=1.e-5
    338        CALL getin('iflag_top_bound',iflag_top_bound)
    339353       CALL getin('tau_top_bound',tau_top_bound)
    340354

  • LMDZ5/trunk/libf/dyn3dpar/leapfrog_p.F

    r1676 r1793  
    904904c      ajout des tendances physiques:
    905905c      ------------------------------
    906          IF (ok_strato) THEN
    907            CALL top_bound_p( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    908          ENDIF
    909        
    910906          CALL addfi_p( dtphys, leapf, forward   ,
    911907     $                  ucov, vcov, teta , q   ,ps ,
    912908     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
    913909
     910         IF (ok_strato) THEN
     911           CALL top_bound_p(vcov,ucov,teta,masse,dtphys)
     912         ENDIF
     913       
    914914c$OMP BARRIER
    915915c$OMP MASTER
    … …  
    10241024        ! Sponge layer (if any)
    10251025        IF (ok_strato) THEN
    1026           ! set dufi,dvfi,... to zero
    1027           ijb=ij_begin
    1028           ije=ij_end
    1029 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    1030           do l=1,llm
    1031             dufi(ijb:ije,l)=0
    1032             dtetafi(ijb:ije,l)=0
    1033             dqfi(ijb:ije,l,1:nqtot)=0
    1034           enddo
    1035 !$OMP END DO
    1036 !$OMP MASTER
    1037           dpfi(ijb:ije)=0
    1038 !$OMP END MASTER
    1039           ijb=ij_begin
    1040           ije=ij_end
    1041           if (pole_sud) ije=ije-iip1
    1042 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    1043           do l=1,llm
    1044             dvfi(ijb:ije,l)=0
    1045           enddo
    1046 !$OMP END DO
    1047 
    1048           CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    1049           CALL addfi_p( dtvr, leapf, forward   ,
    1050      $                  ucov, vcov, teta , q   ,ps ,
    1051      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     1026          CALL top_bound_p(vcov,ucov,teta,masse,dtvr)
    10521027!$OMP BARRIER
    10531028        ENDIF ! of IF (ok_strato)

  • 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.