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/dyn3dmem
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/dyn3dmem/call_calfis_mod.F90

    r1676 r1793  
    321321#endif
    322322
    323     IF (ok_strato) THEN
    324       CALL top_bound_loc( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    325     ENDIF
    326 
    327323#ifdef DEBUG_IO           
    328324    CALL WriteField_u('ucovfi',ucov)
    … …  
    348344    ENDDO
    349345#endif
     346
     347    IF (ok_strato) THEN
     348!      CALL top_bound_loc( vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
     349      CALL top_bound_loc(vcov,ucov,teta,masse,dtphys)
     350    ENDIF
    350351
    351352  !$OMP BARRIER

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

    r1673 r1793  
    11!
    2 ! $Id$
     2! $Id: comconst.h 1671 2012-10-24 07:10:10Z emillour $
    33!
    44!-----------------------------------------------------------------------
    … …  
    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              &
    … …  
    2121      REAL dtdiss ! (s) time step for the dissipation
    2222      REAL rad ! (m) radius of the planet
    23       REAL r ! Reduced Gas constant r=R/mu 
    24              ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol) 
     23      REAL r ! Reduced Gas constant r=R/mu
     24             ! with R=8.31.. J.K-1.mol-1, mu: mol mass of atmosphere (kg/mol)
    2525      REAL cpp   ! Specific heat Cp (J.kg-1.K-1)
    2626      REAL kappa ! kappa=R/Cp
    … …  
    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/dyn3dmem/comvert.h

    r1673 r1793  
    11!
    2 ! $Id$
     2! $Id: comvert.h 1654 2012-09-24 15:07:18Z aslmd $
    33!
    44!-----------------------------------------------------------------------
    … …  
    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/dyn3dmem/conf_gcm.F

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

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

    r1705 r1793  
    11251125        ! Sponge layer (if any)
    11261126        IF (ok_strato) THEN
    1127           ! set dufi,dvfi,... to zero
    1128           ijb=ij_begin
    1129           ije=ij_end
    1130 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    1131           do l=1,llm
    1132             dufi(ijb:ije,l)=0
    1133             dtetafi(ijb:ije,l)=0
    1134             dqfi(ijb:ije,l,1:nqtot)=0
    1135           enddo
    1136 !$OMP END DO
    1137 !$OMP MASTER
    1138           dpfi(ijb:ije)=0
    1139 !$OMP END MASTER
    1140           ijb=ij_begin
    1141           ije=ij_end
    1142           if (pole_sud) ije=ije-iip1
    1143 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    1144           do l=1,llm
    1145             dvfi(ijb:ije,l)=0
    1146           enddo
    1147 !$OMP END DO
    1148 
    1149           CALL top_bound_loc(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
    1150           CALL addfi_loc( dtvr, leapf, forward   ,
    1151      $                  ucov, vcov, teta , q   ,ps ,
    1152      $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
     1127          CALL top_bound_loc(vcov,ucov,teta,masse,dtvr)
    11531128!$OMP BARRIER
    11541129        ENDIF ! of IF (ok_strato)

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