Changeset 1010


Ignore:
Timestamp:
Jul 24, 2013, 1:36:44 PM (11 years ago)
Author:
emillour
Message:

Common dynamics: Improved sponge layer scheme (top_bound):

  • Sponge tendencies are now computed analytically, instead of using a Forward Euler approximation.
  • Sponge tendencies are now added within top_bound.

EM

Location:
trunk
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/DOC/chantiers/commit_importants.log

    r1000 r1010  
    12251225Ehouarn: Some cleanup in the dynamics: ismax.F,ismin.F and cray.F (scopy/ssum)
    12261226moved to 'bibio' (rather than being in dyn3d and dyn3dpar).
     1227
     1228**********************
     1229**** commit_v1010 ****
     1230**********************
     1231Ehouarn: Improved sponge layer scheme (top_bound):
     1232- Sponge tendencies are now computed analytically, instead than using
     1233  a Forward Euler approximation.
     1234- Sponge tendencies are now added within top_bound.

  • trunk/DOC/documentation/top_bound.tex

    r108 r1010  
    3737
    3838\vspace{1cm}
    39 S\'ebastien Lebonnois
     39S\'ebastien Lebonnois , Ehouarn Millour
    4040
    4141\vspace{1cm}
    … …  
    4343\end{center}
    4444
    45 %\section{Theoretical aspects}
     45\section{Theoretical aspects}
     46Because of the inevitable numerical boundary at the top of the model,
     47upward travelling waves are found to non-physically reflect down into the
     48atmosphere.
     49A common remedy to this unwanted behaviour is to apply a sponge layer near
     50the top of the model in order to quench these waves and avoid significant
     51reflection thereof.
     52In practice such quenching is done by adding a dissipative term which forces
     53a relaxation of potential temperature and/or winds of the form:
     54\[
     55 A(t)=A_m+A_0 \exp(-\lambda t )
     56\]
     57Where $A_m$ is the value towards which $A$ is to asymptotically relax, and
     58$\lambda$ is the inverse of the characteristic relaxation time scale.
     59As there is no obvious value of $A_m$ towards which to relax, in practice
     60it is often chosen to be either the zonal average of $A$ (evaluated at time $t$,
     61i.e. conveniently ignoring that $A_m$ then is in fact not time-independent),
     62or zero (at least for winds, since this would have little physical meaning for
     63potential temperature).
    4664
    4765\section{Pratical aspects in the code}
    … …  
    5068flag is set to {\em True} in \textsf{gcm.def}
    5169(this parameter also controls the application of a second step in the
    52 horizontal dissipation).
     70determination of vertical variation of coefficients for
     71the horizontal dissipation, see \textsf{inidissip.F} and
     72\textsf{disspi\_horiz.pdf} document).
    5373
    5474The tendencies for the upper boundary sponge layer are computed separately in
    55 the \textsf{top\_bound.F} routine, called from \textsf{leapfrog.F}.
    56 These tendencies are \textsf{dutop}, \textsf{dvtop} and \textsf{dhtop}, in
    57 unit/s.
     75the \textsf{top\_bound.F} routine (called from \textsf{leapfrog.F}) and
     76added in place.
     77The resulting sponge tendency \textsf{dutop}, in m/s, is also given as an output for
     78diagnostics.
    5879
    5980Three parameters may be adjusted in the \textsf{gcm.def} file:
    … …  
    7495  averaged value.
    7596  \end{itemize}
    76 \item \textsf{tau\_top\_bound}: damping rate (in /s) in the top layer.
     97\item \textsf{tau\_top\_bound}: damping rate ($\lambda$ in equation above,
     98expressed in Hz) in the topmost layer.
    7799\end{itemize}
    78100

  • trunk/LMDZ.COMMON/libf/dyn3d/conf_gcm.F

    r979 r1010  
    343343       CALL getin('dissip_pupstart',dissip_pupstart )
    344344
    345 ! Parametres controlant la couche eponge
    346 ! Actifs uniquement avec ok_strato=y
    347 !   mode = 0 : pas de sponge
    348 !   mode = 1 : u et v -> 0
    349 !   mode = 2 : u et v -> moyenne zonale
    350 !   mode = 3 : u, v et h -> moyenne zonale
     345! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0
     346!                   iflag_top_bound=0 for no sponge
     347!                   iflag_top_bound=1 for sponge over 4 topmost layers
     348!                   iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
    351349       iflag_top_bound=1
     350       CALL getin('iflag_top_bound',iflag_top_bound)
     351
     352! mode_top_bound : fields towards which sponge relaxation will be done:
     353!                  mode_top_bound=0: no relaxation
     354!                  mode_top_bound=1: u and v relax towards 0
     355!                  mode_top_bound=2: u and v relax towards their zonal mean
     356!                  mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
    352357       mode_top_bound=3
     358       CALL getin('mode_top_bound',mode_top_bound)
     359
     360! top_bound sponge : inverse of charactericstic relaxation time scale for sponge
    353361       tau_top_bound=1.e-5
    354        CALL getin('iflag_top_bound',iflag_top_bound)
    355        CALL getin('mode_top_bound',mode_top_bound)
    356362       CALL getin('tau_top_bound',tau_top_bound)
    357363

  • trunk/LMDZ.COMMON/libf/dyn3d/leapfrog.F

    r841 r1010  
    103103
    104104c   tendances de la couche superieure */s
    105       REAL dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
    106       REAL dtetatop(ip1jmp1,llm)
    107       REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
     105c      REAL dvtop(ip1jm,llm)
     106      REAL dutop(ip1jmp1,llm)
     107c      REAL dtetatop(ip1jmp1,llm)
     108c      REAL dqtop(ip1jmp1,llm,nqtot),dptop(ip1jmp1)
    108109
    109110c   TITAN : tendances due au forces de marees */s
    … …  
    216217        dtetadis(:,:)=0.
    217218        dutop(:,:)   =0.
    218         dvtop(:,:)   =0.
    219         dtetatop(:,:)=0.
    220         dqtop(:,:,:) =0.
    221         dptop(:)     =0.
     219c        dvtop(:,:)   =0.
     220c        dtetatop(:,:)=0.
     221c        dqtop(:,:,:) =0.
     222c        dptop(:)     =0.
    222223        dufi(:,:)   =0.
    223224        dvfi(:,:)   =0.
    … …  
    504505c      -------------------
    505506         IF (ok_strato) THEN
    506            CALL top_bound( vcov,ucov,teta,phi,masse,
    507      $                     dutop,dvtop,dtetatop)
    508 c dqtop=0, dptop=0
    509            CALL addfi( dtphys, leapf, forward   ,
    510      $                  ucov, vcov, teta , q   ,ps ,
    511      $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
     507           CALL top_bound(vcov,ucov,teta,masse,dtphys,dutop)
     508           dutop(:,:)=dutop(:,:)/dtphys   ! convert to a tendency in (m/s)/s
    512509         ENDIF
    513510
    … …  
    543540        ! Sponge layer (if any)
    544541        IF (ok_strato) THEN
    545           CALL top_bound(vcov,ucov,teta,phi,
    546      $                   masse,dutop,dvtop,dtetatop)
    547 c dqtop=0, dptop=0
    548           CALL addfi( dtvr, leapf, forward   ,
    549      $                  ucov, vcov, teta , q   ,ps ,
    550      $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
     542           CALL top_bound(vcov,ucov,teta,masse,dtvr,dutop)
     543           dutop(:,:)=dutop(:,:)/dtvr   ! convert to a tendency in (m/s)/s
    551544        ENDIF ! of IF (ok_strato)
    552545      ENDIF ! of IF (iflag_phys.EQ.2)

  • trunk/LMDZ.COMMON/libf/dyn3d/top_bound.F

    r110 r1010  
    1       SUBROUTINE top_bound( vcov,ucov,teta,phi,masse, du,dv,dh )
     1!
     2! $Id: top_bound.F 1793 2013-07-18 07:13:18Z emillour $
     3!
     4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt,ducov)
    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*dt))
     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 phi(iip1,jjp1,llm)                  ! geopotentiel
    38       REAL masse(iip1,jjp1,llm)
    39       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
     64      real,intent(out) :: ducov(iip1,jjp1,llm) ! increment on ucov due to sponge
    4065
    4166c   Local:
    … …  
    4570      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
    4671     
    47       INTEGER NDAMP
    48       PARAMETER (NDAMP=4)
    4972      integer i
    50       REAL,SAVE :: rdamp(llm)
    51 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     73      REAL,SAVE :: rdamp(llm) ! quenching coefficient
     74      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5275
    5376      LOGICAL,SAVE :: first=.true.
    5477
    55       REAL zkm
    5678      INTEGER j,l
    57 
    58 
    59 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    6079     
    6180      if (iflag_top_bound.eq.0) return
    … …  
    6382      if (first) then
    6483         if (iflag_top_bound.eq.1) then
    65 ! couche eponge dans les 4 dernieres couches du modele
    66              rdamp(:)=0.
    67              rdamp(llm)=tau_top_bound
    68              rdamp(llm-1)=tau_top_bound/2.
    69              rdamp(llm-2)=tau_top_bound/4.
    70              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.
    7190         else if (iflag_top_bound.eq.2) then
    72 ! couche eponge dans toutes les couches de pression plus faible que
    73 ! 100 fois la pression de la derniere couche
    74              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
    7594     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    7695         endif
    77          first=.false.
    78          print*,'TOP_BOUND mode',mode_top_bound
    79          print*,'Coeffs pour la couche eponge a l equateur'
    80          print*,'p (Pa)  z(km)  tau (s)'
     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)'
    81104         do l=1,llm
    82105           if (rdamp(l).ne.0.) then
    83             zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
    84           print*,presnivs(l),zkm,1./rdamp(l)
     106             write(lunout,'(6(1pe12.4,1x))')
     107     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     108     &           1./lambda(l),lambda(l)
    85109           endif
    86110         enddo
    87       endif
     111         first=.false.
     112      endif ! of if (first)
    88113
    89114      CALL massbar(masse,massebx,masseby)
    90115
    91 c   mode = 0 : pas de sponge
    92 c   mode = 1 : u et v -> 0
    93 c   mode = 2 : u et v -> moyenne zonale
    94 c   mode = 3 : u, v et h -> moyenne zonale
    95 
     116      ! compute zonal average of vcov and u
    96117      if (mode_top_bound.ge.2) then
    97118       do l=1,llm
    … …  
    100121          zm=0.
    101122          do i=1,iim
    102 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    103 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    104 ! ne varie qu'en latitude
     123! NB: we can work using vcov zonal mean rather than v since the
     124! cv coefficient (which relates the two) only varies with latitudes
    105125            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    106126            zm=zm+masseby(i,j,l)
    … …  
    111131
    112132       do l=1,llm
    113         do j=2,jjm
     133        do j=2,jjm ! excluding poles
    114134          uzon(j,l)=0.
    115135          zm=0.
    … …  
    121141        enddo
    122142       enddo
    123       else
    124           vzon(:,:)=0.
    125           uzon(:,:)=0.
    126       endif
     143      else ! ucov and vcov will relax towards 0
     144        vzon(:,:)=0.
     145        uzon(:,:)=0.
     146      endif ! of if (mode_top_bound.ge.2)
    127147
     148      ! compute zonal average of potential temperature, if necessary
    128149      if (mode_top_bound.ge.3) then
    129150       do l=1,llm
    130         do j=2,jjm
     151        do j=2,jjm ! excluding poles
    131152          zm=0.
    132153          tzon(j,l)=0.
    … …  
    138159        enddo
    139160       enddo
    140       endif
    141 
    142 C   AMORTISSEMENTS LINEAIRES:
     161      endif ! of if (mode_top_bound.ge.3)
    143162
    144163      if (mode_top_bound.ge.1) then
     164       ! Apply sponge quenching on vcov:
    145165       do l=1,llm
    146166        do i=1,iip1
    147167          do j=1,jjm
    148             dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     168            vcov(i,j,l)=vcov(i,j,l)
     169     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    149170          enddo
    150171        enddo
    151172       enddo
    152173
     174       ! Apply sponge quenching on ucov:
    153175       do l=1,llm
    154176        do i=1,iip1
    155           do j=2,jjm
    156             du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     177          do j=2,jjm ! excluding poles
     178            ducov(i,j,l)=-rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     179            ucov(i,j,l)=ucov(i,j,l)
     180     &                  +ducov(i,j,l)
    157181          enddo
    158182        enddo
    159183       enddo
    160       endif
     184      endif ! of if (mode_top_bound.ge.1)
    161185
    162186      if (mode_top_bound.ge.3) then
     187       ! Apply sponge quenching on teta:
    163188       do l=1,llm
    164189        do i=1,iip1
    165           do j=2,jjm
    166             dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     190          do j=2,jjm ! excluding poles
     191            teta(i,j,l)=teta(i,j,l)
     192     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
    167193          enddo
    168194        enddo
    169195       enddo
    170       endif
    171      
    172       RETURN
     196      endif ! of if (mode_top_bound.ge.3)
     197   
    173198      END

  • trunk/LMDZ.COMMON/libf/dyn3dpar/conf_gcm.F

    r979 r1010  
    370370       CALL getin('dissip_pupstart',dissip_pupstart )
    371371
    372 ! Parametres controlant la couche eponge
    373 ! Actifs uniquement avec ok_strato=y
    374 !   mode = 0 : pas de sponge
    375 !   mode = 1 : u et v -> 0
    376 !   mode = 2 : u et v -> moyenne zonale
    377 !   mode = 3 : u, v et h -> moyenne zonale
     372! top_bound sponge: only active if ok_strato=.true. and iflag_top_bound!=0
     373!                   iflag_top_bound=0 for no sponge
     374!                   iflag_top_bound=1 for sponge over 4 topmost layers
     375!                   iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
    378376       iflag_top_bound=1
     377       CALL getin('iflag_top_bound',iflag_top_bound)
     378
     379! mode_top_bound : fields towards which sponge relaxation will be done:
     380!                  mode_top_bound=0: no relaxation
     381!                  mode_top_bound=1: u and v relax towards 0
     382!                  mode_top_bound=2: u and v relax towards their zonal mean
     383!                  mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
    379384       mode_top_bound=3
     385       CALL getin('mode_top_bound',mode_top_bound)
     386
     387! top_bound sponge : inverse of charactericstic relaxation time scale for sponge
    380388       tau_top_bound=1.e-5
    381        CALL getin('iflag_top_bound',iflag_top_bound)
    382        CALL getin('mode_top_bound',mode_top_bound)
    383389       CALL getin('tau_top_bound',tau_top_bound)
    384390

  • trunk/LMDZ.COMMON/libf/dyn3dpar/leapfrog_p.F

    r953 r1010  
    113113
    114114c   tendances top_bound (sponge layer)
    115       REAL,SAVE :: dvtop(ip1jm,llm),dutop(ip1jmp1,llm)
    116       REAL,SAVE :: dtetatop(ip1jmp1,llm)
    117       REAL,SAVE :: dptop(ip1jmp1)
    118       REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
     115c      REAL,SAVE :: dvtop(ip1jm,llm)
     116      REAL,SAVE :: dutop(ip1jmp1,llm)
     117c      REAL,SAVE :: dtetatop(ip1jmp1,llm)
     118c      REAL,SAVE :: dptop(ip1jmp1)
     119c      REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqtop
    119120
    120121c   TITAN : tendances due au forces de marees */s
    … …  
    247248            ALLOCATE(dqfi(ip1jmp1,llm,nqtot))
    248249            ALLOCATE(dqfi_tmp(iip1,llm,nqtot))
    249             ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
     250c            ALLOCATE(dqtop(ip1jmp1,llm,nqtot))
    250251         END IF
    251252c$OMP END MASTER     
    … …  
    262263        dtetadis(:,:)=0.
    263264        dutop(:,:)   =0.
    264         dvtop(:,:)   =0.
    265         dtetatop(:,:)=0.
    266         dqtop(:,:,:) =0.
    267         dptop(:)     =0.
     265c        dvtop(:,:)   =0.
     266c        dtetatop(:,:)=0.
     267c        dqtop(:,:,:) =0.
     268c        dptop(:)     =0.
    268269        dufi(:,:)   =0.
    269270        dvfi(:,:)   =0.
    … …  
    10011002c      -------------------
    10021003         IF (ok_strato) THEN
    1003            CALL top_bound_p( vcov,ucov,teta,phi,masse,
    1004      $                       dutop,dvtop,dtetatop)
    1005            CALL addfi_p( dtphys, leapf, forward   ,
    1006      $                  ucov, vcov, teta , q   ,ps ,
    1007      $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
    1008 
    1009          ENDIF
     1004           CALL top_bound_p(vcov,ucov,teta,masse,dtphys,
     1005     $                       dutop)
     1006        ijb=ij_begin
     1007        ije=ij_end
     1008c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
     1009        DO l=1,llm
     1010          dutop(ijb:ije,l)=dutop(ijb:ije,l)/dtphys ! convert to a tendency in (m/s)/s
     1011        ENDDO
     1012c$OMP END DO NOWAIT       
     1013         ENDIF ! of IF (ok_strato)
    10101014       
    10111015c$OMP BARRIER
    … …  
    11241128        ! Sponge layer (if any)
    11251129        IF (ok_strato) THEN
    1126           ! set dutop,dvtop,... to zero
     1130          CALL top_bound_p(vcov,ucov,teta,masse,dtvr,
     1131     $                     dutop)
    11271132          ijb=ij_begin
    11281133          ije=ij_end
    1129 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    1130           do l=1,llm
    1131             dutop(ijb:ije,l)=0
    1132             dtetatop(ijb:ije,l)=0
    1133             dqtop(ijb:ije,l,1:nqtot)=0
    1134           enddo
    1135 !$OMP END DO
    1136 !$OMP SINGLE
    1137           dptop(ijb:ije)=0
    1138 !$OMP END SINGLE
    1139           ijb=ij_begin
    1140           ije=ij_end
    1141           if (pole_sud) ije=ije-iip1
    1142 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    1143           do l=1,llm
    1144             dvtop(ijb:ije,l)=0
    1145           enddo
    1146 !$OMP END DO
    1147 
    1148           CALL top_bound_p(vcov,ucov,teta,phi,masse,
    1149      $                     dutop,dvtop,dtetatop)
    1150           CALL addfi_p( dtvr, leapf, forward   ,
    1151      $                  ucov, vcov, teta , q   ,ps ,
    1152      $                 dutop, dvtop, dtetatop , dqtop ,dptop  )
     1134c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)       
     1135          DO l=1,llm
     1136            dutop(ijb:ije,l)=dutop(ijb:ije,l)/dtvr ! convert to a tendency in (m/s)/s
     1137          ENDDO
     1138c$OMP END DO NOWAIT       
    11531139!$OMP BARRIER
    11541140        ENDIF ! of IF (ok_strato)

  • trunk/LMDZ.COMMON/libf/dyn3dpar/top_bound_p.F

    r124 r1010  
    1       SUBROUTINE top_bound_p( vcov,ucov,teta,phi,masse, du,dv,dh )
     1!
     2! $Id: top_bound_p.F 1793 2013-07-18 07:13:18Z emillour $
     3!
     4      SUBROUTINE top_bound_p(vcov,ucov,teta,masse,dt,ducov)
    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*dt)
     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 phi(iip1,jjp1,llm)                  ! geopotentiel
    38       REAL masse(iip1,jjp1,llm)
    39       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
     65      real,intent(out) :: ducov(iip1,jjp1,llm) ! increment on ucov due to sponge
    4066
    4167c   Local:
    … …  
    4470      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
    4571     
    46       INTEGER NDAMP
    47       PARAMETER (NDAMP=4)
    4872      integer i
    49       REAL,SAVE :: rdamp(llm)
    50 !     &   (/(0., i =1,llm-NDAMP),0.125E-5,.25E-5,.5E-5,1.E-5/)
     73      REAL,SAVE :: rdamp(llm) ! quenching coefficient
     74      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
    5175      LOGICAL,SAVE :: first=.true.
    52      
    53       REAL zkm
    5476      INTEGER j,l,jjb,jje
    5577
    5678
    5779      if (iflag_top_bound == 0) return
     80
    5881      if (first) then
    5982c$OMP BARRIER
    6083c$OMP MASTER
    6184         if (iflag_top_bound == 1) then
    62 ! couche eponge dans les 4 dernieres couches du modele
    63              rdamp(:)=0.
    64              rdamp(llm)=tau_top_bound
    65              rdamp(llm-1)=tau_top_bound/2.
    66              rdamp(llm-2)=tau_top_bound/4.
    67              rdamp(llm-3)=tau_top_bound/8.
     85! sponge quenching over the topmost 4 atmospheric layers
     86             lambda(:)=0.
     87             lambda(llm)=tau_top_bound
     88             lambda(llm-1)=tau_top_bound/2.
     89             lambda(llm-2)=tau_top_bound/4.
     90             lambda(llm-3)=tau_top_bound/8.
    6891         else if (iflag_top_bound == 2) then
    69 ! couche eponge dans toutes les couches de pression plus faible que
    70 ! 100 fois la pression de la derniere couche
    71              rdamp(:)=tau_top_bound
     92! sponge quenching over topmost layers down to pressures which are
     93! higher than 100 times the topmost layer pressure
     94             lambda(:)=tau_top_bound
    7295     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
    7396         endif
    74          first=.false.
    75          print*,'TOP_BOUND mode',mode_top_bound
    76          print*,'Coeffs pour la couche eponge a l equateur'
    77          print*,'p (Pa)  z(km)  tau (s)'
     97
     98! quenching coefficient rdamp(:)
     99!         rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
     100         rdamp(:)=1.-exp(-lambda(:)*dt)
     101
     102         write(lunout,*)'TOP_BOUND mode',mode_top_bound
     103         write(lunout,*)'Sponge layer coefficients'
     104         write(lunout,*)'p (Pa)  z(km)  tau(s)   1./tau (Hz)'
    78105         do l=1,llm
    79106           if (rdamp(l).ne.0.) then
    80             zkm        = phi(iip1/2,jjp1/2,l)/(1000*g)
    81           print*,presnivs(l),zkm,1./rdamp(l)
     107             write(lunout,'(6(1pe12.4,1x))')
     108     &        presnivs(l),log(preff/presnivs(l))*scaleheight,
     109     &           1./lambda(l),lambda(l)
    82110           endif
    83111         enddo
     112         first=.false.
    84113c$OMP END MASTER
    85114c$OMP BARRIER
    86       endif
     115      endif ! of if (first)
    87116
    88117
    89118      CALL massbar_p(masse,massebx,masseby)
    90119
    91 c   mode = 0 : pas de sponge
    92 c   mode = 1 : u et v -> 0
    93 c   mode = 2 : u et v -> moyenne zonale
    94 c   mode = 3 : u, v et h -> moyenne zonale
    95 
    96 C POUR V
    97 
    98 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    99 
    100       jjb=jj_begin
    101       jje=jj_end
    102       IF (pole_sud) jje=jj_end-1
    103 
     120      ! compute zonal average of vcov (or set it to zero)
    104121      if (mode_top_bound.ge.2) then
    105 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     122       jjb=jj_begin
     123       jje=jj_end
     124       IF (pole_sud) jje=jj_end-1
     125c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    106126       do l=1,llm
    107127        do j=jjb,jje
    … …  
    109129          vzon(j,l)=0
    110130          do i=1,iim
    111 ! Rm: on peut travailler directement avec la moyenne zonale de vcov
    112 ! plutot qu'avec celle de v car le coefficient cv qui relie les deux
    113 ! ne varie qu'en latitude
     131! NB: we can work using vcov zonal mean rather than v since the
     132! cv coefficient (which relates the two) only varies with latitudes
    114133            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
    115134            zm=zm+masseby(i,j,l)
    … …  
    118137        enddo
    119138       enddo
    120 !$OMP END DO NOWAIT   
     139c$OMP END DO NOWAIT   
    121140      else
    122 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
    123        do l=1,llm
    124         do j=jjb,jje
    125           vzon(j,l)=0.
    126         enddo
    127        enddo
    128 !$OMP END DO NOWAIT
    129       endif
    130 
    131 C   AMORTISSEMENTS LINEAIRES:
    132 
    133       if (mode_top_bound.ge.1) then
    134 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    135        do l=1,llm
    136         do j=jjb,jje
    137           do i=1,iip1
    138             dv(i,j,l)= -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
    139           enddo
    140         enddo
    141        enddo
    142 !$OMP END DO NOWAIT
    143       endif
    144 
    145 C POUR U ET H
    146 
    147 C  CALCUL DES CHAMPS EN MOYENNE ZONALE:
    148 
    149       jjb=jj_begin
    150       jje=jj_end
    151       IF (pole_nord) jjb=jj_begin+1
    152       IF (pole_sud)  jje=jj_end-1
    153 
     141c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     142       do l=1,llm
     143         vzon(:,l)=0.
     144       enddo
     145c$OMP END DO NOWAIT
     146      endif ! of if (mode_top_bound.ge.2)
     147
     148      ! compute zonal average of u (or set it to zero)
    154149      if (mode_top_bound.ge.2) then
    155 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     150       jjb=jj_begin
     151       jje=jj_end
     152       IF (pole_nord) jjb=jj_begin+1
     153       IF (pole_sud)  jje=jj_end-1
     154c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    156155       do l=1,llm
    157156        do j=jjb,jje
    … …  
    165164        enddo
    166165       enddo
    167 !$OMP END DO NOWAIT
     166c$OMP END DO NOWAIT
    168167      else
    169 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    170        do l=1,llm
    171         do j=jjb,jje
    172           uzon(j,l)=0.
    173         enddo
    174        enddo
    175 !$OMP END DO NOWAIT
    176       endif
    177 
     168c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     169       do l=1,llm
     170         uzon(:,l)=0.
     171       enddo
     172c$OMP END DO NOWAIT
     173      endif ! of if (mode_top_bound.ge.2)
     174
     175      ! compute zonal average of potential temperature, if necessary
    178176      if (mode_top_bound.ge.3) then
    179 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     177       jjb=jj_begin
     178       jje=jj_end
     179       IF (pole_nord) jjb=jj_begin+1
     180       IF (pole_sud)  jje=jj_end-1
     181c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)   
    180182       do l=1,llm
    181183        do j=jjb,jje
    … …  
    189191        enddo
    190192       enddo
    191 !$OMP END DO NOWAIT
    192       endif
    193 
    194 C   AMORTISSEMENTS LINEAIRES:
     193c$OMP END DO NOWAIT
     194      endif ! of if (mode_top_bound.ge.3)
    195195
    196196      if (mode_top_bound.ge.1) then
    197 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     197       ! Apply sponge quenching on vcov:
     198       jjb=jj_begin
     199       jje=jj_end
     200       IF (pole_sud) jje=jj_end-1
     201
     202c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    198203       do l=1,llm
    199204        do j=jjb,jje
    200205          do i=1,iip1
    201             du(i,j,l)= -rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
    202           enddo
    203         enddo
    204        enddo
    205 !$OMP END DO NOWAIT
    206       endif
    207      
    208       if (mode_top_bound.ge.3) then
    209 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
     206            vcov(i,j,l)=vcov(i,j,l)
     207     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
     208          enddo
     209        enddo
     210       enddo
     211c$OMP END DO NOWAIT
     212
     213       ! Apply sponge quenching on ucov:
     214       jjb=jj_begin
     215       jje=jj_end
     216       IF (pole_nord) jjb=jj_begin+1
     217       IF (pole_sud)  jje=jj_end-1
     218
     219c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
    210220       do l=1,llm
    211221        do j=jjb,jje
    212222          do i=1,iip1
    213             dh(i,j,l)= -rdamp(l)*(teta(i,j,l)-tzon(j,l))
    214           enddo
    215         enddo
    216        enddo
    217 !$OMP END DO NOWAIT
    218       endif     
    219 
    220 !$OMP BARRIER
    221       RETURN
     223            ducov(i,j,l)=-rdamp(l)*(ucov(i,j,l)-cu(i,j)*uzon(j,l))
     224            ucov(i,j,l)=ucov(i,j,l)
     225     &                  +ducov(i,j,l)
     226          enddo
     227       enddo
     228       enddo
     229c$OMP END DO NOWAIT
     230      endif ! of if (mode_top_bound.ge.1)
     231
     232      if (mode_top_bound.ge.3) then   
     233       ! Apply sponge quenching on teta:
     234       jjb=jj_begin
     235       jje=jj_end
     236       IF (pole_nord) jjb=jj_begin+1
     237       IF (pole_sud)  jje=jj_end-1
     238
     239c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)     
     240       do l=1,llm
     241        do j=jjb,jje
     242          do i=1,iip1
     243            teta(i,j,l)=teta(i,j,l)
     244     &                  -rdamp(l)*(teta(i,j,l)-tzon(j,l))
     245          enddo
     246       enddo
     247       enddo
     248c$OMP END DO NOWAIT
     249      endif ! of if (mode_top_bond.ge.3)
     250
    222251      END
Note: See TracChangeset for help on using the changeset viewer.