Context Navigation

top_bound_loc.F @ 5456

Last change on this file since 5456 was 1823, checked in by Ehouarn Millour, 11 years ago
Remplacement de parallel.F90 (en conflit avec orchidée) par parallel_lmdz.F90. UG ......................................... Renaming parallel.F90 (conflicting with orchidée) into parallel_lmdz.F90. UG
File size: 7.3 KB

Line
1	!
2	! $Id: $
3	!
4	SUBROUTINE top_bound_loc(vcov,ucov,teta,masse,dt)
5	USE parallel_lmdz
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