Home

Context Navigation

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

Last change on this file since 1394 was 1383, checked in by boeske, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 18.6 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 1383 2014-04-30 12:17:54Z keck $
27	!
28	! 1382 2014-04-30 12:15:41Z boeske
29	! Changed the weighting factor that is used in the summation of nudging
30	! tendencies for profile data output from weight_pres to weight_substep,
31	! added Neumann boundary conditions for profile data output of nudging terms at
32	! nzt+1
33	!
34	! 1380 2014-04-28 12:40:45Z heinze
35	! Subroutine nudge_ref added to account for proper upper scalar boundary
36	! conditions in case of nudging
37	!
38	! 1365 2014-04-22 15:03:56Z boeske
39	! Variable t renamed nt, variable currtnudge renamed tmp_tnudge,
40	! summation of nudging tendencies for data output added
41	! +sums_ls_l, tmp_tend
42	! Added new subroutine calc_tnudge, which calculates the current nudging time
43	! scale at each time step
44	!
45	! 1355 2014-04-10 10:21:29Z heinze
46	! Error message specified.
47	!
48	! 1353 2014-04-08 15:21:23Z heinze
49	! REAL constants provided with KIND-attribute
50	!
51	! 1320 2014-03-20 08:40:49Z raasch
52	! ONLY-attribute added to USE-statements,
53	! kind-parameters added to all INTEGER and REAL declaration statements,
54	! kinds are defined in new module kinds,
55	! old module precision_kind is removed,
56	! revision history before 2012 removed,
57	! comment fields (!:) to be used for variable explanations added to
58	! all variable declaration statements
59	!
60	! 1318 2014-03-17 13:35:16Z raasch
61	! module interfaces removed
62	!
63	! 1268 2013-12-12 09:47:53Z heinze
64	! bugfix: argument of calc_mean_profile corrected
65	!
66	! 1251 2013-11-07 08:14:30Z heinze
67	! bugfix: calculate dtm and dtp also in vector version
68	!
69	! 1249 2013-11-06 10:45:47Z heinze
70	! remove call of user module
71	! reformatting
72	!
73	! 1241 2013-10-30 11:36:58Z heinze
74	! Initial revision
75	!
76	! Description:
77	! ------------
78	! Nudges u, v, pt and q to given profiles on a relaxation timescale tnudge.
79	! Profiles are read in from NUDGIN_DATA. Code is based on Neggers et al. (2012)
80	! and also part of DALES and UCLA-LES.
81	!--------------------------------------------------------------------------------!
82
83	PRIVATE
84	PUBLIC init_nudge, calc_tnudge, nudge, nudge_ref
85	SAVE
86
87	INTERFACE nudge
88	MODULE PROCEDURE nudge
89	MODULE PROCEDURE nudge_ij
90	END INTERFACE nudge
91
92	CONTAINS
93
94	SUBROUTINE init_nudge
95
96	USE arrays_3d, &
97	ONLY: ptnudge, qnudge, timenudge, tmp_tnudge, tnudge, unudge, &
98	vnudge, wnudge, zu
99
100	USE control_parameters, &
101	ONLY: dt_3d, lptnudge, lqnudge, lunudge, lvnudge, lwnudge, &
102	message_string, ntnudge
103
104	USE indices, &
105	ONLY: nzb, nzt
106
107	USE kinds
108
109	IMPLICIT NONE
110
111
112	INTEGER(iwp) :: finput = 90 !:
113	INTEGER(iwp) :: ierrn !:
114	INTEGER(iwp) :: k !:
115	INTEGER(iwp) :: nt !:
116
117	CHARACTER(1) :: hash !:
118
119	REAL(wp) :: highheight !:
120	REAL(wp) :: highqnudge !:
121	REAL(wp) :: highptnudge !:
122	REAL(wp) :: highunudge !:
123	REAL(wp) :: highvnudge !:
124	REAL(wp) :: highwnudge !:
125	REAL(wp) :: hightnudge !:
126
127	REAL(wp) :: lowheight !:
128	REAL(wp) :: lowqnudge !:
129	REAL(wp) :: lowptnudge !:
130	REAL(wp) :: lowunudge !:
131	REAL(wp) :: lowvnudge !:
132	REAL(wp) :: lowwnudge !:
133	REAL(wp) :: lowtnudge !:
134
135	REAL(wp) :: fac !:
136
137	ALLOCATE( ptnudge(nzb:nzt+1,1:ntnudge), qnudge(nzb:nzt+1,1:ntnudge), &
138	tnudge(nzb:nzt+1,1:ntnudge), unudge(nzb:nzt+1,1:ntnudge), &
139	vnudge(nzb:nzt+1,1:ntnudge), wnudge(nzb:nzt+1,1:ntnudge) )
140
141	ALLOCATE( tmp_tnudge(nzb:nzt) )
142
143	ALLOCATE( timenudge(0:ntnudge) )
144
145	ptnudge = 0.0_wp; qnudge = 0.0_wp; tnudge = 0.0_wp; unudge = 0.0_wp
146	vnudge = 0.0_wp; wnudge = 0.0_wp; timenudge = 0.0_wp
147	!
148	!-- Initialize array tmp_nudge with a current nudging time scale of 6 hours
149	tmp_tnudge = 21600.0_wp
150
151	nt = 0
152	OPEN ( finput, FILE='NUDGING_DATA', STATUS='OLD', &
153	FORM='FORMATTED', IOSTAT=ierrn )
154
155	IF ( ierrn /= 0 ) THEN
156	message_string = 'file NUDGING_DATA does not exist'
157	CALL message( 'nudging', 'PA0365', 1, 2, 0, 6, 0 )
158	ENDIF
159
160	ierrn = 0
161
162	rloop:DO
163	nt = nt + 1
164	hash = "#"
165	ierrn = 1 ! not zero
166	!
167	!-- Search for the next line consisting of "# time",
168	!-- from there onwards the profiles will be read
169	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
170
171	READ ( finput, *, IOSTAT=ierrn ) hash, timenudge(nt)
172	IF ( ierrn < 0 ) EXIT rloop
173
174	ENDDO
175
176	ierrn = 0
177	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowtnudge, lowunudge, &
178	lowvnudge, lowwnudge , lowptnudge, &
179	lowqnudge
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	ierrn = 0
187	READ ( finput, *, IOSTAT=ierrn ) highheight, hightnudge, highunudge, &
188	highvnudge, highwnudge , highptnudge, &
189	highqnudge
190
191	IF ( ierrn /= 0 ) THEN
192	message_string = 'errors in file NUDGING_DATA'
193	CALL message( 'nudging', 'PA0366', 1, 2, 0, 6, 0 )
194	ENDIF
195
196	DO k = nzb, nzt+1
197	IF ( highheight < zu(k) ) THEN
198	lowheight = highheight
199	lowtnudge = hightnudge
200	lowunudge = highunudge
201	lowvnudge = highvnudge
202	lowwnudge = highwnudge
203	lowptnudge = highptnudge
204	lowqnudge = highqnudge
205
206	ierrn = 0
207	READ ( finput, *, IOSTAT=ierrn ) highheight , hightnudge , &
208	highunudge , highvnudge , &
209	highwnudge , highptnudge, &
210	highqnudge
211	IF (ierrn /= 0 ) THEN
212	WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm is ',&
213	'higher than the maximum height in NUDING_DATA which ', &
214	'is ', lowheight, 'm. Interpolation on PALM ', &
215	'grid is not possible.'
216	CALL message( 'nudging', 'PA0364', 1, 2, 0, 6, 0 )
217	ENDIF
218	ENDIF
219
220	!
221	!-- Interpolation of prescribed profiles in space
222
223	fac = ( highheight - zu(k) ) / ( highheight - lowheight )
224
225	tnudge(k,nt) = fac * lowtnudge + ( 1.0_wp - fac ) * hightnudge
226	unudge(k,nt) = fac * lowunudge + ( 1.0_wp - fac ) * highunudge
227	vnudge(k,nt) = fac * lowvnudge + ( 1.0_wp - fac ) * highvnudge
228	wnudge(k,nt) = fac * lowwnudge + ( 1.0_wp - fac ) * highwnudge
229	ptnudge(k,nt) = fac * lowptnudge + ( 1.0_wp - fac ) * highptnudge
230	qnudge(k,nt) = fac * lowqnudge + ( 1.0_wp - fac ) * highqnudge
231	ENDDO
232
233	ENDDO rloop
234
235	CLOSE ( finput )
236
237	!
238	!-- Prevent nudging if nudging profiles exhibt too small values
239	!-- not used so far
240	lptnudge = ANY( ABS( ptnudge ) > 1.0e-8_wp )
241	lqnudge = ANY( ABS( qnudge ) > 1.0e-8_wp )
242	lunudge = ANY( ABS( unudge ) > 1.0e-8_wp )
243	lvnudge = ANY( ABS( vnudge ) > 1.0e-8_wp )
244	lwnudge = ANY( ABS( wnudge ) > 1.0e-8_wp )
245
246	END SUBROUTINE init_nudge
247
248
249	SUBROUTINE calc_tnudge ( time )
250
251	USE arrays_3d, &
252	ONLY: timenudge, tmp_tnudge, tnudge
253
254	USE control_parameters, &
255	ONLY: dt_3d
256
257	USE indices, &
258	ONLY: nzb, nzt
259
260	USE kinds
261
262	IMPLICIT NONE
263
264
265	REAL(wp) :: dtm !:
266	REAL(wp) :: dtp !:
267	REAL(wp) :: time !:
268
269	INTEGER(iwp) :: k !:
270	INTEGER(iwp) :: nt !:
271
272	nt = 1
273	DO WHILE ( time > timenudge(nt) )
274	nt = nt+1
275	ENDDO
276	IF ( time /= timenudge(1) ) THEN
277	nt = nt-1
278	ENDIF
279
280	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
281	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
282
283	DO k = nzb, nzt
284	tmp_tnudge(k) = MAX( dt_3d, tnudge(k,nt) * dtp + tnudge(k,nt+1) * dtm )
285	ENDDO
286
287	END SUBROUTINE calc_tnudge
288
289	!------------------------------------------------------------------------------!
290	! Call for all grid points
291	!------------------------------------------------------------------------------!
292	SUBROUTINE nudge ( time, prog_var )
293
294	USE arrays_3d, &
295	ONLY: pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
296	u, unudge, v, vnudge
297
298	USE control_parameters, &
299	ONLY: dt_3d, intermediate_timestep_count, message_string
300
301	USE indices, &
302	ONLY: nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
303
304	USE kinds
305
306	USE statistics, &
307	ONLY: hom, sums_ls_l, weight_substep
308
309	IMPLICIT NONE
310
311	CHARACTER (LEN=*) :: prog_var !:
312
313	REAL(wp) :: tmp_tend !:
314	REAL(wp) :: dtm !:
315	REAL(wp) :: dtp !:
316	REAL(wp) :: time !:
317
318	INTEGER(iwp) :: i !:
319	INTEGER(iwp) :: j !:
320	INTEGER(iwp) :: k !:
321	INTEGER(iwp) :: nt !:
322
323
324	nt = 1
325	DO WHILE ( time > timenudge(nt) )
326	nt = nt+1
327	ENDDO
328	IF ( time /= timenudge(1) ) THEN
329	nt = nt-1
330	ENDIF
331
332	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
333	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
334
335	SELECT CASE ( prog_var )
336
337	CASE ( 'u' )
338
339	DO i = nxl, nxr
340	DO j = nys, nyn
341
342	DO k = nzb_u_inner(j,i)+1, nzt
343
344	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
345	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
346
347	tend(k,j,i) = tend(k,j,i) + tmp_tend
348
349	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend * &
350	weight_substep(intermediate_timestep_count)
351	ENDDO
352
353	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
354
355	ENDDO
356	ENDDO
357
358	CASE ( 'v' )
359
360	DO i = nxl, nxr
361	DO j = nys, nyn
362
363	DO k = nzb_u_inner(j,i)+1, nzt
364
365	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
366	vnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
367
368	tend(k,j,i) = tend(k,j,i) + tmp_tend
369
370	sums_ls_l(k,7) = sums_ls_l(k,7) + tmp_tend * &
371	weight_substep(intermediate_timestep_count)
372	ENDDO
373
374	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
375
376	ENDDO
377	ENDDO
378
379	CASE ( 'pt' )
380
381	DO i = nxl, nxr
382	DO j = nys, nyn
383
384	DO k = nzb_u_inner(j,i)+1, nzt
385
386	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
387	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
388
389	tend(k,j,i) = tend(k,j,i) + tmp_tend
390
391	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend * &
392	weight_substep(intermediate_timestep_count)
393	ENDDO
394
395	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
396
397	ENDDO
398	ENDDO
399
400	CASE ( 'q' )
401
402	DO i = nxl, nxr
403	DO j = nys, nyn
404
405	DO k = nzb_u_inner(j,i)+1, nzt
406
407	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
408	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
409
410	tend(k,j,i) = tend(k,j,i) + tmp_tend
411
412	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend * &
413	weight_substep(intermediate_timestep_count)
414	ENDDO
415
416	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
417
418	ENDDO
419	ENDDO
420
421	CASE DEFAULT
422	message_string = 'unknown prognostic variable "' // prog_var // '"'
423	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
424
425	END SELECT
426
427	END SUBROUTINE nudge
428
429
430	!------------------------------------------------------------------------------!
431	! Call for grid point i,j
432	!------------------------------------------------------------------------------!
433
434	SUBROUTINE nudge_ij( i, j, time, prog_var )
435
436	USE arrays_3d, &
437	ONLY: pt, ptnudge, q, qnudge, tend, timenudge, tmp_tnudge, tnudge, &
438	u, unudge, v, vnudge
439
440	USE control_parameters, &
441	ONLY: dt_3d, intermediate_timestep_count, message_string
442
443	USE indices, &
444	ONLY: nxl, nxr, nys, nyn, nzb, nzb_u_inner, nzt
445
446	USE kinds
447
448	USE statistics, &
449	ONLY: hom, sums_ls_l, weight_substep
450
451	IMPLICIT NONE
452
453
454	CHARACTER (LEN=*) :: prog_var !:
455
456	REAL(wp) :: tmp_tend !:
457	REAL(wp) :: dtm !:
458	REAL(wp) :: dtp !:
459	REAL(wp) :: time !:
460
461	INTEGER(iwp) :: i !:
462	INTEGER(iwp) :: j !:
463	INTEGER(iwp) :: k !:
464	INTEGER(iwp) :: nt !:
465
466
467	nt = 1
468	DO WHILE ( time > timenudge(nt) )
469	nt = nt+1
470	ENDDO
471	IF ( time /= timenudge(1) ) THEN
472	nt = nt-1
473	ENDIF
474
475	dtm = ( time - timenudge(nt) ) / ( timenudge(nt+1) - timenudge(nt) )
476	dtp = ( timenudge(nt+1) - time ) / ( timenudge(nt+1) - timenudge(nt) )
477
478	SELECT CASE ( prog_var )
479
480	CASE ( 'u' )
481
482	DO k = nzb_u_inner(j,i)+1, nzt
483
484	tmp_tend = - ( hom(k,1,1,0) - ( unudge(k,nt) * dtp + &
485	unudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
486
487	tend(k,j,i) = tend(k,j,i) + tmp_tend
488
489	sums_ls_l(k,6) = sums_ls_l(k,6) + tmp_tend &
490	* weight_substep(intermediate_timestep_count)
491	ENDDO
492
493	sums_ls_l(nzt+1,6) = sums_ls_l(nzt,6)
494
495	CASE ( 'v' )
496
497	DO k = nzb_u_inner(j,i)+1, nzt
498
499	tmp_tend = - ( hom(k,1,2,0) - ( vnudge(k,nt) * dtp + &
500	vnudge(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,7) = sums_ls_l(k,7) + tmp_tend &
505	* weight_substep(intermediate_timestep_count)
506	ENDDO
507
508	sums_ls_l(nzt+1,7) = sums_ls_l(nzt,7)
509
510	CASE ( 'pt' )
511
512	DO k = nzb_u_inner(j,i)+1, nzt
513
514	tmp_tend = - ( hom(k,1,4,0) - ( ptnudge(k,nt) * dtp + &
515	ptnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
516
517	tend(k,j,i) = tend(k,j,i) + tmp_tend
518
519	sums_ls_l(k,4) = sums_ls_l(k,4) + tmp_tend &
520	* weight_substep(intermediate_timestep_count)
521	ENDDO
522
523	sums_ls_l(nzt+1,4) = sums_ls_l(nzt,4)
524
525
526	CASE ( 'q' )
527
528	DO k = nzb_u_inner(j,i)+1, nzt
529
530	tmp_tend = - ( hom(k,1,41,0) - ( qnudge(k,nt) * dtp + &
531	qnudge(k,nt+1) * dtm ) ) / tmp_tnudge(k)
532
533	tend(k,j,i) = tend(k,j,i) + tmp_tend
534
535	sums_ls_l(k,5) = sums_ls_l(k,5) + tmp_tend &
536	* weight_substep(intermediate_timestep_count)
537	ENDDO
538
539	sums_ls_l(nzt+1,5) = sums_ls_l(nzt,5)
540
541	CASE DEFAULT
542	message_string = 'unknown prognostic variable "' // prog_var // '"'
543	CALL message( 'nudge', 'PA0367', 1, 2, 0, 6, 0 )
544
545	END SELECT
546
547
548	END SUBROUTINE nudge_ij
549
550
551	SUBROUTINE nudge_ref ( time )
552
553	USE arrays_3d, &
554	ONLY: time_vert, ptnudge, pt_init, qnudge, q_init
555
556	USE kinds
557
558
559	IMPLICIT NONE
560
561	INTEGER(iwp) :: nt !:
562
563	REAL(wp) :: fac !:
564	REAL(wp), INTENT(in) :: time !:
565
566	!
567	!-- Interpolation in time of NUDGING_DATA for pt_init and q_init. This is
568	!-- needed for correct upper boundary conditions for pt and q and in case that
569	! large scale subsidence as well as scalar Rayleigh-damping are used
570	nt = 1
571	DO WHILE ( time > time_vert(nt) )
572	nt = nt + 1
573	ENDDO
574	IF ( time /= time_vert(nt) ) THEN
575	nt = nt - 1
576	ENDIF
577
578	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
579
580	pt_init = ptnudge(:,nt) + fac * ( ptnudge(:,nt+1) - ptnudge(:,nt) )
581	q_init = qnudge(:,nt) + fac * ( qnudge(:,nt+1) - qnudge(:,nt) )
582
583	END SUBROUTINE nudge_ref
584
585
586	END MODULE nudge_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |