Home

Context Navigation

source: palm/trunk/SOURCE/nudging.f90 @ 1366

Last change on this file since 1366 was 1366, checked in by boeske, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 16.7 KB

Line
1	MODULE nudge_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: nudging.f90 1366 2014-04-22 15:06:33Z boeske $
27	!
28	! 1365 2014-04-22 15:03:56Z boeske
29	! Variable t renamed nt, variable currtnudge renamed tmp_tnudge,
30	! summation of nudging tendencies for data output added
31	! +sums_ls_l, tmp_tend
32	! Added new subroutine calc_tnudge, which calculates the current nudging time
33	! scale at each time step
34	!
35	! 1355 2014-04-10 10:21:29Z heinze
36	! Error message specified.
37	!
38	! 1353 2014-04-08 15:21:23Z heinze
39	! REAL constants provided with KIND-attribute
40	!
41	! 1320 2014-03-20 08:40:49Z raasch
42	! ONLY-attribute added to USE-statements,
43	! kind-parameters added to all INTEGER and REAL declaration statements,
44	! kinds are defined in new module kinds,
45	! old module precision_kind is removed,
46	! revision history before 2012 removed,
47	! comment fields (!:) to be used for variable explanations added to
48	! all variable declaration statements
49	!
50	! 1318 2014-03-17 13:35:16Z raasch
51	! module interfaces removed
52	!
53	! 1268 2013-12-12 09:47:53Z heinze
54	! bugfix: argument of calc_mean_profile corrected
55	!
56	! 1251 2013-11-07 08:14:30Z heinze
57	! bugfix: calculate dtm and dtp also in vector version
58	!
59	! 1249 2013-11-06 10:45:47Z heinze
60	! remove call of user module
61	! reformatting
62	!
63	! 1241 2013-10-30 11:36:58Z heinze
64	! Initial revision
65	!
66	! Description:
67	! ------------
68	! Nudges u, v, pt and q to given profiles on a relaxation timescale tnudge.
69	! Profiles are read in from NUDGIN_DATA. Code is based on Neggers et al. (2012)
70	! and also part of DALES and UCLA-LES.
71	!--------------------------------------------------------------------------------!
72
73	PRIVATE
74	PUBLIC init_nudge, calc_tnudge, nudge
75	SAVE
76
77	INTERFACE nudge
78	MODULE PROCEDURE nudge
79	MODULE PROCEDURE nudge_ij
80	END INTERFACE nudge
81
82	CONTAINS
83
84	SUBROUTINE init_nudge
85
86	USE arrays_3d, &
87	ONLY: ptnudge, qnudge, timenudge, tmp_tnudge, tnudge, unudge, &
88	vnudge, wnudge, zu
89
90	USE control_parameters, &
91	ONLY: dt_3d, lptnudge, lqnudge, lunudge, lvnudge, lwnudge, &
92	message_string, ntnudge
93
94	USE indices, &
95	ONLY: nzb, nzt
96
97	USE kinds
98
99	IMPLICIT NONE
100
101
102	INTEGER(iwp) :: finput = 90 !:
103	INTEGER(iwp) :: ierrn !:
104	INTEGER(iwp) :: k !:
105	INTEGER(iwp) :: nt !:
106
107	CHARACTER(1) :: hash !:
108
109	REAL(wp) :: highheight !:
110	REAL(wp) :: highqnudge !:
111	REAL(wp) :: highptnudge !:
112	REAL(wp) :: highunudge !:
113	REAL(wp) :: highvnudge !:
114	REAL(wp) :: highwnudge !:
115	REAL(wp) :: hightnudge !:
116
117	REAL(wp) :: lowheight !:
118	REAL(wp) :: lowqnudge !:
119	REAL(wp) :: lowptnudge !:
120	REAL(wp) :: lowunudge !:
121	REAL(wp) :: lowvnudge !:
122	REAL(wp) :: lowwnudge !:
123	REAL(wp) :: lowtnudge !:
124
125	REAL(wp) :: fac !:
126
127	ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
128	tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge), &
129	vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge) )
130
131	ALLOCATE( tmp_tnudge(nzb:nzt) )
132
133	ALLOCATE( timenudge(0:ntnudge) )
134
135	ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
136	vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
137	!
138	!-- Initialize array tmp_nudge with a current nudging time scale of 6 hours
139	tmp_tnudge = 21600.0_wp
140
141	nt = 0
142	OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
143	FORM='FORMATTED', IOSTAT=ierrn )
144
145	IF ( ierrn /= 0 ) THEN
146	message_string = 'file NUDGING_DATA does not exist'
147	CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
148	ENDIF
149
150	ierrn = 0
151
152	rloop:DO
153	nt = nt + 1
154	hash = "#"
155	ierrn = 1 ! not zero
156	!
157	!-- Search for the next line consisting of "# time",
158	!-- from there onwards the profiles will be read
159	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
160
161	READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
162	IF ( ierrn < 0 ) EXIT rloop
163
164	ENDDO
165
166	ierrn = 0
167	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge, &
168	lowvnudge, lowwnudge , lowptnudge, &
169	lowqnudge
170
171	IF ( ierrn /= 0 ) THEN
172	message_string = 'errors in file NUDGING_DATA'
173	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
174	ENDIF
175
176	ierrn = 0
177	READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge, &
178	highvnudge, highwnudge , highptnudge, &
179	highqnudge
180
181	IF ( ierrn /= 0 ) THEN
182	message_string = 'errors in file NUDGING_DATA'
183	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
184	ENDIF
185
186	DO k = nzb, nzt+1
187	IF ( highheight < zu(k) ) THEN
188	lowheight = highheight
189	lowtnudge = hightnudge
190	lowunudge = highunudge
191	lowvnudge = highvnudge
192	lowwnudge = highwnudge
193	lowptnudge = highptnudge
194	lowqnudge = highqnudge
195
196	ierrn = 0
197	READ ( finput, *, IOSTAT=ierrn ) highheight , hightnudge , &
198	highunudge , highvnudge , &
199	highwnudge , highptnudge, &
200	highqnudge
201	IF (ierrn /= 0 ) THEN
202	WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm is ',&
203	'higher than the maximum height in NUDING_DATA which ', &
204	'is ', lowheight, 'm. Interpolation on PALM ', &
205	'grid is not possible.'
206	CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
207	ENDIF
208	ENDIF
209
210	!
211	!-- Interpolation of prescribed profiles in space
212
213	fac = ( highheight - zu(k) ) / ( highheight - lowheight )
214
215	tnudge(k,nt) = fac * lowtnudge + ( 1.0_wp - fac ) * hightnudge
216	unudge(k,nt) = fac * lowunudge + ( 1.0_wp - fac ) * highunudge
217	vnudge(k,nt) = fac * lowvnudge + ( 1.0_wp - fac ) * highvnudge
218	wnudge(k,nt) = fac * lowwnudge + ( 1.0_wp - fac ) * highwnudge
219	ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
220	qnudge(k,nt) = fac * lowqnudge + ( 1.0_wp - fac ) * highqnudge
221	ENDDO
222
223	ENDDO rloop
224
225	CLOSE ( finput )
226
227	!
228	!-- Prevent nudging if nudging profiles exhibt too small values
229	!-- not used so far
230	lptnudge = ANY( ABS( ptnudge ) > 1.0e-8_wp )
231	lqnudge = ANY( ABS( qnudge ) > 1.0e-8_wp )
232	lunudge = ANY( ABS( unudge ) > 1.0e-8_wp )
233	lvnudge = ANY( ABS( vnudge ) > 1.0e-8_wp )
234	lwnudge = ANY( ABS( wnudge ) > 1.0e-8_wp )
235
236	END SUBROUTINE init_nudge
237
238
239	SUBROUTINE calc_tnudge ( time )
240
241	USE arrays_3d, &
242	ONLY: timenudge, tmp_tnudge, tnudge
243
244	USE control_parameters, &
245	ONLY: dt_3d
246
247	USE indices, &
248	ONLY: nzb, nzt
249
250	USE kinds
251
252	IMPLICIT NONE
253
254
255	REAL(wp) :: dtm !:
256	REAL(wp) :: dtp !:
257	REAL(wp) :: time !:
258
259	INTEGER(iwp) :: k !:
260	INTEGER(iwp) :: nt !:
261
262	nt = 1
263	DO WHILE ( time > timenudge(nt) )
264	nt = nt+1
265	ENDDO
266	IF ( time /= timenudge(1) ) THEN
267	nt = nt-1
268	ENDIF
269
270	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
271	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
272
273	DO k = nzb, nzt
274	tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
275	ENDDO
276
277	END SUBROUTINE calc_tnudge
278
279	!------------------------------------------------------------------------------!
280	! Call for all grid points
281	!------------------------------------------------------------------------------!
282	SUBROUTINE nudge ( time, prog_var )
283
284	USE arrays_3d, &
285	ONLY: pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
286	u, unudge, v, vnudge
287
288	USE control_parameters, &
289	ONLY: dt_3d, intermediate_timestep_count, message_string
290
291	USE indices, &
292	ONLY: nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
293
294	USE kinds
295
296	USE statistics, &
297	ONLY: hom, sums_ls_l, weight_pres
298
299	IMPLICIT NONE
300
301	CHARACTER (LEN=*) :: prog_var !:
302
303	REAL(wp) :: tmp_tend !:
304	REAL(wp) :: dtm !:
305	REAL(wp) :: dtp !:
306	REAL(wp) :: time !:
307
308	INTEGER(iwp) :: i !:
309	INTEGER(iwp) :: j !:
310	INTEGER(iwp) :: k !:
311	INTEGER(iwp) :: nt !:
312
313
314	nt = 1
315	DO WHILE ( time > timenudge(nt) )
316	nt = nt+1
317	ENDDO
318	IF ( time /= timenudge(1) ) THEN
319	nt = nt-1
320	ENDIF
321
322	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
323	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
324
325	SELECT CASE ( prog_var )
326
327	CASE ( 'u' )
328
329	DO i = nxl, nxr
330	DO j = nys, nyn
331	DO k = nzb_u_inner(j,i)+1, nzt
332
333	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
334	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
335
336	tend(k,j,i) = tend(k,j,i) + tmp_tend
337
338	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend * &
339	weight_pres(intermediate_timestep_count)
340	ENDDO
341	ENDDO
342	ENDDO
343
344	CASE ( 'v' )
345
346	DO i = nxl, nxr
347	DO j = nys, nyn
348	DO k = nzb_u_inner(j,i)+1, nzt
349
350	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
351	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
352
353	tend(k,j,i) = tend(k,j,i) + tmp_tend
354
355	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend * &
356	weight_pres(intermediate_timestep_count)
357	ENDDO
358	ENDDO
359	ENDDO
360
361	CASE ( 'pt' )
362
363	DO i = nxl, nxr
364	DO j = nys, nyn
365	DO k = nzb_u_inner(j,i)+1, nzt
366
367	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
368	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
369
370	tend(k,j,i) = tend(k,j,i) + tmp_tend
371
372	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend * &
373	weight_pres(intermediate_timestep_count)
374	ENDDO
375	ENDDO
376	ENDDO
377
378	CASE ( 'q' )
379
380	DO i = nxl, nxr
381	DO j = nys, nyn
382	DO k = nzb_u_inner(j,i)+1, nzt
383
384	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
385	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
386
387	tend(k,j,i) = tend(k,j,i) + tmp_tend
388
389	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend * &
390	weight_pres(intermediate_timestep_count)
391	ENDDO
392	ENDDO
393	ENDDO
394
395	CASE DEFAULT
396	message_string = 'unknown prognostic variable "' // prog_var // '"'
397	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
398
399	END SELECT
400
401	END SUBROUTINE nudge
402
403
404	!------------------------------------------------------------------------------!
405	! Call for grid point i,j
406	!------------------------------------------------------------------------------!
407
408	SUBROUTINE nudge_ij( i, j, time, prog_var )
409
410	USE arrays_3d, &
411	ONLY: pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
412	u, unudge, v, vnudge
413
414	USE control_parameters, &
415	ONLY: dt_3d, intermediate_timestep_count, message_string
416
417	USE indices, &
418	ONLY: nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
419
420	USE kinds
421
422	USE statistics, &
423	ONLY: hom, sums_ls_l, weight_pres
424
425	IMPLICIT NONE
426
427
428	CHARACTER (LEN=*) :: prog_var !:
429
430	REAL(wp) :: tmp_tend !:
431	REAL(wp) :: dtm !:
432	REAL(wp) :: dtp !:
433	REAL(wp) :: time !:
434
435	INTEGER(iwp) :: i !:
436	INTEGER(iwp) :: j !:
437	INTEGER(iwp) :: k !:
438	INTEGER(iwp) :: nt !:
439
440
441	nt = 1
442	DO WHILE ( time > timenudge(nt) )
443	nt = nt+1
444	ENDDO
445	IF ( time /= timenudge(1) ) THEN
446	nt = nt-1
447	ENDIF
448
449	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
450	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
451
452	SELECT CASE ( prog_var )
453
454	CASE ( 'u' )
455
456	DO k = nzb_u_inner(j,i)+1, nzt
457
458	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
459	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
460
461	tend(k,j,i) = tend(k,j,i) + tmp_tend
462
463	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend &
464	* weight_pres(intermediate_timestep_count)
465	ENDDO
466
467	CASE ( 'v' )
468
469	DO k = nzb_u_inner(j,i)+1, nzt
470
471	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
472	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
473
474	tend(k,j,i) = tend(k,j,i) + tmp_tend
475
476	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend &
477	* weight_pres(intermediate_timestep_count)
478	ENDDO
479
480
481	CASE ( 'pt' )
482
483	DO k = nzb_u_inner(j,i)+1, nzt
484
485	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
486	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
487
488	tend(k,j,i) = tend(k,j,i) + tmp_tend
489
490	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend &
491	* weight_pres(intermediate_timestep_count)
492	ENDDO
493
494
495	CASE ( 'q' )
496
497	DO k = nzb_u_inner(j,i)+1, nzt
498
499	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
500	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
501
502	tend(k,j,i) = tend(k,j,i) + tmp_tend
503
504	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend &
505	* weight_pres(intermediate_timestep_count)
506	ENDDO
507
508	CASE DEFAULT
509	message_string = 'unknown prognostic variable "' // prog_var // '"'
510	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
511
512	END SELECT
513
514
515	END SUBROUTINE nudge_ij
516
517	END MODULE nudge_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |