Changeset 1793 for LMDZ5/trunk


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
Files:
16 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

  • 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

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