Context Navigation

top_bound_p.F @ 3537

Last change on this file since 3537 was 1442, checked in by slebonnois, 10 years ago
SL: update of the Venus GCM, + corrections on routines used for newstart/start2archive for Titan and Venus, + some modifications on tools
File size: 7.4 KB

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

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

Download in other formats:

Original Format