source: LMDZ5/branches/LF-private/libf/dyn3d/top_bound.F @ 4585

Last change on this file since 4585 was 1793, checked in by Ehouarn Millour, 11 years ago

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

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 5.9 KB
Line 
1!
2! $Id: top_bound.F 1793 2013-07-18 07:13:18Z fairhead $
3!
4      SUBROUTINE top_bound(vcov,ucov,teta,masse,dt)
5      IMPLICIT NONE
6c
7#include "dimensions.h"
8#include "paramet.h"
9#include "comconst.h"
10#include "comvert.h"
11#include "comgeom2.h"
12
13
14c ..  DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
15C     F. LOTT DEC. 2006
16c                                 (  10/12/06  )
17
18c=======================================================================
19c
20c   Auteur:  F. LOTT 
21c   -------
22c
23c   Objet:
24c   ------
25c
26c   Dissipation linéaire (ex top_bound de la physique)
27c
28c=======================================================================
29
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
53#include "comdissipn.h"
54#include "iniprint.h"
55
56c   Arguments:
57c   ----------
58
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
65c   Local:
66c   ------
67
68      REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
69      REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
70     
71      integer i
72      REAL,SAVE :: rdamp(llm) ! quenching coefficient
73      real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
74
75      LOGICAL,SAVE :: first=.true.
76
77      INTEGER j,l
78     
79      if (iflag_top_bound.eq.0) return
80
81      if (first) then
82         if (iflag_top_bound.eq.1) then
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.
89         else if (iflag_top_bound.eq.2) then
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
93     s       *max(presnivs(llm)/presnivs(:)-0.01,0.)
94         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
110         first=.false.
111      endif ! of if (first)
112
113      CALL massbar(masse,massebx,masseby)
114
115      ! compute zonal average of vcov and u
116      if (mode_top_bound.ge.2) then
117       do l=1,llm
118        do j=1,jjm
119          vzon(j,l)=0.
120          zm=0.
121          do i=1,iim
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
124            vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
125            zm=zm+masseby(i,j,l)
126          enddo
127          vzon(j,l)=vzon(j,l)/zm
128        enddo
129       enddo
130
131       do l=1,llm
132        do j=2,jjm ! excluding poles
133          uzon(j,l)=0.
134          zm=0.
135          do i=1,iim
136            uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
137            zm=zm+massebx(i,j,l)
138          enddo
139          uzon(j,l)=uzon(j,l)/zm
140        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)
146
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
151          zm=0.
152          tzon(j,l)=0.
153          do i=1,iim
154            tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
155            zm=zm+masse(i,j,l)
156          enddo
157          tzon(j,l)=tzon(j,l)/zm
158        enddo
159       enddo
160      endif ! of if (mode_top_bound.ge.3)
161
162      if (mode_top_bound.ge.1) then
163       ! Apply sponge quenching on vcov:
164       do l=1,llm
165        do i=1,iip1
166          do j=1,jjm
167            vcov(i,j,l)=vcov(i,j,l)
168     &                  -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
169          enddo
170        enddo
171       enddo
172
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   
196      END
Note: See TracBrowser for help on using the repository browser.