Context Navigation

top_bound_loc.F @ 1793

Last change on this file since 1793 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

File size: 7.3 KB

Line
1	!
2	! $Id: $
3	!
4	SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt)
5	USE parallel
6	IMPLICIT NONE
7	c
8	#include "dimensions.h"
9	#include "paramet.h"
10	#include "comconst.h"
11	#include "comvert.h"
12	#include "comgeom2.h"
13
14
15	c .. DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
16	C F. LOTT DEC. 2006
17	c ( 10/12/06 )
18
19	c=======================================================================
20	c
21	c Auteur: F. LOTT
22	c -------
23	c
24	c Objet:
25	c ------
26	c
27	c Dissipation linéaire (ex top_bound de la physique)
28	c
29	c=======================================================================
30
31	! top_bound sponge layer model:
32	! Quenching is modeled as: A(t)=Am+A0exp(-lambdat)
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(-lambdat))
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
53
54	#include "comdissipn.h"
55	#include "iniprint.h"
56
57	c Arguments:
58	c ----------
59
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)
68
69	c Local:
70	c ------
71	REAL massebx(iip1,jjb_u:jje_u,llm),masseby(iip1,jjb_v:jje_v,llm)
72	REAL zm
73	REAL uzon(jjb_u:jje_u,llm),vzon(jjb_v:jje_v,llm)
74	REAL tzon(jjb_u:jje_u,llm)
75
76	integer i
77	REAL,SAVE :: rdamp(llm)
78	real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
79	LOGICAL,SAVE :: first=.true.
80	INTEGER j,l,jjb,jje
81
82
83	if (iflag_top_bound == 0) return
84
85	if (first) then
86	c$OMP BARRIER
87	c$OMP MASTER
88	if (iflag_top_bound == 1) then
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.
95	else if (iflag_top_bound == 2) then
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
99	s *max(presnivs(llm)/presnivs(:)-0.01,0.)
100	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
116	first=.false.
117	c$OMP END MASTER
118	c$OMP BARRIER
119	endif ! of if (first)
120
121
122	CALL massbar_loc(masse,massebx,masseby)
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
129	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
130	do l=1,llm
131	do j=jjb,jje
132	zm=0.
133	vzon(j,l)=0
134	do i=1,iim
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
137	vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
138	zm=zm+masseby(i,j,l)
139	enddo
140	vzon(j,l)=vzon(j,l)/zm
141	enddo
142	enddo
143	c$OMP END DO NOWAIT
144	else
145	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
146	do l=1,llm
147	vzon(:,l)=0.
148	enddo
149	c$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
158	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
159	do l=1,llm
160	do j=jjb,jje
161	uzon(j,l)=0.
162	zm=0.
163	do i=1,iim
164	uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
165	zm=zm+massebx(i,j,l)
166	enddo
167	uzon(j,l)=uzon(j,l)/zm
168	enddo
169	enddo
170	c$OMP END DO NOWAIT
171	else
172	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
173	do l=1,llm
174	uzon(:,l)=0.
175	enddo
176	c$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
185	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
186	do l=1,llm
187	do j=jjb,jje
188	zm=0.
189	tzon(j,l)=0.
190	do i=1,iim
191	tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
192	zm=zm+masse(i,j,l)
193	enddo
194	tzon(j,l)=tzon(j,l)/zm
195	enddo
196	enddo
197	c$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
206	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
207	do l=1,llm
208	do j=jjb,jje
209	do i=1,iip1
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
215	c$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
223	c$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
232	c$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
242	c$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
251	c$OMP END DO NOWAIT
252	endif ! of if (mode_top_bond.ge.3)
253
254	END

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format