Context Navigation

top_bound_loc.F @ 2598

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

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

Download in other formats:

Original Format