Context Navigation

top_bound_p.F @ 3626

Last change on this file since 3626 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 Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 7.3 KB

Line
1	!
2	! $Id: top_bound_p.F 2600 2016-07-23 05:45:38Z lguez $
3	!
4	SUBROUTINE top_bound_p(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,jjp1,llm) ! covariant zonal wind
63	real,intent(inout) :: vcov(iip1,jjm,llm) ! covariant meridional wind
64	real,intent(inout) :: teta(iip1,jjp1,llm) ! potential temperature
65	real,intent(in) :: masse(iip1,jjp1,llm) ! mass of atmosphere
66	real,intent(in) :: dt ! time step (s) of sponge model
67
68	c Local:
69	c ------
70	REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm),zm
71	REAL uzon(jjp1,llm),vzon(jjm,llm),tzon(jjp1,llm)
72
73	integer i
74	REAL,SAVE :: rdamp(llm) ! quenching coefficient
75	real,save :: lambda(llm) ! inverse or quenching time scale (Hz)
76	LOGICAL,SAVE :: first=.true.
77	INTEGER j,l,jjb,jje
78
79
80	if (iflag_top_bound == 0) return
81
82	if (first) then
83	c$OMP BARRIER
84	c$OMP MASTER
85	if (iflag_top_bound == 1) then
86	! sponge quenching over the topmost 4 atmospheric layers
87	lambda(:)=0.
88	lambda(llm)=tau_top_bound
89	lambda(llm-1)=tau_top_bound/2.
90	lambda(llm-2)=tau_top_bound/4.
91	lambda(llm-3)=tau_top_bound/8.
92	else if (iflag_top_bound == 2) then
93	! sponge quenching over topmost layers down to pressures which are
94	! higher than 100 times the topmost layer pressure
95	lambda(:)=tau_top_bound
96	s *max(presnivs(llm)/presnivs(:)-0.01,0.)
97	endif
98
99	! quenching coefficient rdamp(:)
100	! rdamp(:)=dt*lambda(:) ! Explicit Euler approx.
101	rdamp(:)=1.-exp(-lambda(:)*dt)
102
103	write(lunout,*)'TOP_BOUND mode',mode_top_bound
104	write(lunout,*)'Sponge layer coefficients'
105	write(lunout,*)'p (Pa) z(km) tau(s) 1./tau (Hz)'
106	do l=1,llm
107	if (rdamp(l).ne.0.) then
108	write(lunout,'(6(1pe12.4,1x))')
109	& presnivs(l),log(preff/presnivs(l))*scaleheight,
110	& 1./lambda(l),lambda(l)
111	endif
112	enddo
113	first=.false.
114	c$OMP END MASTER
115	c$OMP BARRIER
116	endif ! of if (first)
117
118
119	CALL massbar_p(masse,massebx,masseby)
120
121	! compute zonal average of vcov (or set it to zero)
122	if (mode_top_bound.ge.2) then
123	jjb=jj_begin
124	jje=jj_end
125	IF (pole_sud) jje=jj_end-1
126	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
127	do l=1,llm
128	do j=jjb,jje
129	zm=0.
130	vzon(j,l)=0
131	do i=1,iim
132	! NB: we can work using vcov zonal mean rather than v since the
133	! cv coefficient (which relates the two) only varies with latitudes
134	vzon(j,l)=vzon(j,l)+vcov(i,j,l)*masseby(i,j,l)
135	zm=zm+masseby(i,j,l)
136	enddo
137	vzon(j,l)=vzon(j,l)/zm
138	enddo
139	enddo
140	c$OMP END DO NOWAIT
141	else
142	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
143	do l=1,llm
144	vzon(:,l)=0.
145	enddo
146	c$OMP END DO NOWAIT
147	endif ! of if (mode_top_bound.ge.2)
148
149	! compute zonal average of u (or set it to zero)
150	if (mode_top_bound.ge.2) then
151	jjb=jj_begin
152	jje=jj_end
153	IF (pole_nord) jjb=jj_begin+1
154	IF (pole_sud) jje=jj_end-1
155	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
156	do l=1,llm
157	do j=jjb,jje
158	uzon(j,l)=0.
159	zm=0.
160	do i=1,iim
161	uzon(j,l)=uzon(j,l)+massebx(i,j,l)*ucov(i,j,l)/cu(i,j)
162	zm=zm+massebx(i,j,l)
163	enddo
164	uzon(j,l)=uzon(j,l)/zm
165	enddo
166	enddo
167	c$OMP END DO NOWAIT
168	else
169	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
170	do l=1,llm
171	uzon(:,l)=0.
172	enddo
173	c$OMP END DO NOWAIT
174	endif ! of if (mode_top_bound.ge.2)
175
176	! compute zonal average of potential temperature, if necessary
177	if (mode_top_bound.ge.3) then
178	jjb=jj_begin
179	jje=jj_end
180	IF (pole_nord) jjb=jj_begin+1
181	IF (pole_sud) jje=jj_end-1
182	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
183	do l=1,llm
184	do j=jjb,jje
185	zm=0.
186	tzon(j,l)=0.
187	do i=1,iim
188	tzon(j,l)=tzon(j,l)+teta(i,j,l)*masse(i,j,l)
189	zm=zm+masse(i,j,l)
190	enddo
191	tzon(j,l)=tzon(j,l)/zm
192	enddo
193	enddo
194	c$OMP END DO NOWAIT
195	endif ! of if (mode_top_bound.ge.3)
196
197	if (mode_top_bound.ge.1) then
198	! Apply sponge quenching on vcov:
199	jjb=jj_begin
200	jje=jj_end
201	IF (pole_sud) jje=jj_end-1
202
203	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
204	do l=1,llm
205	do j=jjb,jje
206	do i=1,iip1
207	vcov(i,j,l)=vcov(i,j,l)
208	& -rdamp(l)*(vcov(i,j,l)-vzon(j,l))
209	enddo
210	enddo
211	enddo
212	c$OMP END DO NOWAIT
213
214	! Apply sponge quenching on ucov:
215	jjb=jj_begin
216	jje=jj_end
217	IF (pole_nord) jjb=jj_begin+1
218	IF (pole_sud) jje=jj_end-1
219
220	c$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
221	do l=1,llm
222	do j=jjb,jje
223	do i=1,iip1
224	ucov(i,j,l)=ucov(i,j,l)
225	& -rdamp(l)(ucov(i,j,l)-cu(i,j)uzon(j,l))
226	enddo
227	enddo
228	enddo
229	c$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
239	c$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
248	c$OMP END DO NOWAIT
249	endif ! of if (mode_top_bond.ge.3)
250
251	END

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

Download in other formats:

Original Format