Context Navigation

top_bound_loc.F @ 4379

Last change on this file since 4379 was 2600, checked in by Ehouarn Millour, 8 years ago
Cleanup in the dynamics: turn comvert.h into module comvert_mod.F90 EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 7.4 KB

Line
1	!
2	! $Id: $
3	!
4	SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt)
5	USE parallel_lmdz
6	USE comconst_mod, ONLY: iflag_top_bound, mode_top_bound,
7	& tau_top_bound
8	USE comvert_mod, ONLY: presnivs, preff, scaleheight
9
10	IMPLICIT NONE
11	c
12	include "dimensions.h"
13	include "paramet.h"
14	include "comgeom2.h"
15
16
17	c .. DISSIPATION LINEAIRE A HAUT NIVEAU, RUN MESO,
18	C F. LOTT DEC. 2006
19	c ( 10/12/06 )
20
21	c=======================================================================
22	c
23	c Auteur: F. LOTT
24	c -------
25	c
26	c Objet:
27	c ------
28	c
29	c Dissipation linéaire (ex top_bound de la physique)
30	c
31	c=======================================================================
32
33	! top_bound sponge layer model:
34	! Quenching is modeled as: A(t)=Am+A0exp(-lambdat)
35	! where Am is the zonal average of the field (or zero), and lambda the inverse
36	! of the characteristic quenching/relaxation time scale
37	! Thus, assuming Am to be time-independent, field at time t+dt is given by:
38	! A(t+dt)=A(t)-(A(t)-Am)(1-exp(-lambdat))
39	! Moreover lambda can be a function of model level (see below), and relaxation
40	! can be toward the average zonal field or just zero (see below).
41
42	! NB: top_bound sponge is only called from leapfrog if ok_strato=.true.
43
44	! sponge parameters: (loaded/set in conf_gcm.F ; stored in comconst_mod)
45	! iflag_top_bound=0 for no sponge
46	! iflag_top_bound=1 for sponge over 4 topmost layers
47	! iflag_top_bound=2 for sponge from top to ~1% of top layer pressure
48	! mode_top_bound=0: no relaxation
49	! mode_top_bound=1: u and v relax towards 0
50	! mode_top_bound=2: u and v relax towards their zonal mean
51	! mode_top_bound=3: u,v and pot. temp. relax towards their zonal mean
52	! tau_top_bound : inverse of charactericstic relaxation time scale at
53	! the topmost layer (Hz)
54
55
56	#include "comdissipn.h"
57	#include "iniprint.h"
58
59	c Arguments:
60	c ----------
61
62	real,intent(inout) :: ucov(iip1,jjb_u:jje_u,llm) ! covariant zonal wind
63	real,intent(inout) :: vcov(iip1,jjb_v:jje_v,llm) ! covariant meridional wind
64	real,intent(inout) :: teta(iip1,jjb_u:jje_u,llm) ! potential temperature
65	real,intent(in) :: masse(iip1,jjb_u:jje_u,llm) ! mass of atmosphere
66	real,intent(in) :: dt ! time step (s) of sponge model
67
68	! REAL dv(iip1,jjb_v:jje_v,llm),du(iip1,jjb_u:jje_u,llm)
69	! REAL dh(iip1,jjb_u:jje_u,llm)
70
71	c Local:
72	c ------
73	REAL massebx(iip1,jjb_u:jje_u,llm),masseby(iip1,jjb_v:jje_v,llm)
74	REAL zm
75	REAL uzon(jjb_u:jje_u,llm),vzon(jjb_v:jje_v,llm)
76	REAL tzon(jjb_u:jje_u,llm)
77
78	integer i
79	REAL,SAVE :: rdamp(llm)
80	real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
81	LOGICAL,SAVE :: first=.true.
82	INTEGER j,l,jjb,jje
83
84
85	if (iflag_top_bound == 0) return
86
87	if (first) then
88	c$OMP BARRIER
89	c$OMP MASTER
90	if (iflag_top_bound == 1) then
91	! sponge quenching over the topmost 4 atmospheric layers
92	lambda(:)=0.
93	lambda(llm)=tau_top_bound
94	lambda(llm-1)=tau_top_bound/2.
95	lambda(llm-2)=tau_top_bound/4.
96	lambda(llm-3)=tau_top_bound/8.
97	else if (iflag_top_bound == 2) then
98	! sponge quenching over topmost layers down to pressures which are
99	! higher than 100 times the topmost layer pressure
100	lambda(:)=tau_top_bound
101	s *max(presnivs(llm)/presnivs(:)-0.01,0.)
102	endif
103
104	! quenching coefficient rdamp(:)
105	! rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
106	rdamp(:)=1.-exp(-lambda(:)*dt)
107
108	write(lunout,*)'TOP_BOUND mode',mode_top_bound
109	write(lunout,*)'Sponge layer coefficients'
110	write(lunout,*)'p (Pa) z(km) tau(s) 1./tau (Hz)'
111	do l=1,llm
112	if (rdamp(l).ne.0.) then
113	write(lunout,'(6(1pe12.4,1x))')
114	& presnivs(l),log(preff/presnivs(l))*scaleheight,
115	& 1./lambda(l),lambda(l)
116	endif
117	enddo
118	first=.false.
119	c$OMP END MASTER
120	c$OMP BARRIER
121	endif ! of if (first)
122
123
124	CALL massbar_loc(masse,massebx,masseby)
125
126	! compute zonal average of vcov (or set it to zero)
127	if (mode_top_bound.ge.2) then
128	jjb=jj_begin
129	jje=jj_end
130	IF (pole_sud) jje=jj_end-1
131	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
132	do l=1,llm
133	do j=jjb,jje
134	zm=0.
135	vzon(j,l)=0
136	do i=1,iim
137	! NB: we can work using vcov zonal mean rather than v since the
138	! cv coefficient (which relates the two) only varies with latitudes
139	vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
140	zm=zm+masseby(i,j,l)
141	enddo
142	vzon(j,l)=vzon(j,l)/zm
143	enddo
144	enddo
145	c$OMP END DO NOWAIT
146	else
147	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
148	do l=1,llm
149	vzon(:,l)=0.
150	enddo
151	c$OMP END DO NOWAIT
152	endif ! of if (mode_top_bound.ge.2)
153
154	! compute zonal average of u (or set it to zero)
155	if (mode_top_bound.ge.2) then
156	jjb=jj_begin
157	jje=jj_end
158	IF (pole_nord) jjb=jj_begin+1
159	IF (pole_sud) jje=jj_end-1
160	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
161	do l=1,llm
162	do j=jjb,jje
163	uzon(j,l)=0.
164	zm=0.
165	do i=1,iim
166	uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
167	zm=zm+massebx(i,j,l)
168	enddo
169	uzon(j,l)=uzon(j,l)/zm
170	enddo
171	enddo
172	c$OMP END DO NOWAIT
173	else
174	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
175	do l=1,llm
176	uzon(:,l)=0.
177	enddo
178	c$OMP END DO NOWAIT
179	endif ! of if (mode_top_bound.ge.2)
180
181	! compute zonal average of potential temperature, if necessary
182	if (mode_top_bound.ge.3) then
183	jjb=jj_begin
184	jje=jj_end
185	IF (pole_nord) jjb=jj_begin+1
186	IF (pole_sud) jje=jj_end-1
187	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
188	do l=1,llm
189	do j=jjb,jje
190	zm=0.
191	tzon(j,l)=0.
192	do i=1,iim
193	tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
194	zm=zm+masse(i,j,l)
195	enddo
196	tzon(j,l)=tzon(j,l)/zm
197	enddo
198	enddo
199	c$OMP END DO NOWAIT
200	endif ! of if (mode_top_bound.ge.3)
201
202	if (mode_top_bound.ge.1) then
203	! Apply sponge quenching on vcov:
204	jjb=jj_begin
205	jje=jj_end
206	IF (pole_sud) jje=jj_end-1
207
208	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
209	do l=1,llm
210	do j=jjb,jje
211	do i=1,iip1
212	vcov(i,j,l)=vcov(i,j,l)
213	& -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
214	enddo
215	enddo
216	enddo
217	c$OMP END DO NOWAIT
218
219	! Apply sponge quenching on ucov:
220	jjb=jj_begin
221	jje=jj_end
222	IF (pole_nord) jjb=jj_begin+1
223	IF (pole_sud) jje=jj_end-1
224
225	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
226	do l=1,llm
227	do j=jjb,jje
228	do i=1,iip1
229	ucov(i,j,l)=ucov(i,j,l)
230	& -rdamp(l)(ucov(i,j,l)-cu(i,j)uzon(j,l))
231	enddo
232	enddo
233	enddo
234	c$OMP END DO NOWAIT
235	endif ! of if (mode_top_bound.ge.1)
236
237	if (mode_top_bound.ge.3) then
238	! Apply sponge quenching on teta:
239	jjb=jj_begin
240	jje=jj_end
241	IF (pole_nord) jjb=jj_begin+1
242	IF (pole_sud) jje=jj_end-1
243
244	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
245	do l=1,llm
246	do j=jjb,jje
247	do i=1,iip1
248	teta(i,j,l)=teta(i,j,l)
249	& -rdamp(l)*(teta(i,j,l)-tzon(j,l))
250	enddo
251	enddo
252	enddo
253	c$OMP END DO NOWAIT
254	endif ! of if (mode_top_bond.ge.3)
255
256	END

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

Download in other formats:

Original Format