Home

Context Navigation

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

Last change on this file since 1754 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 19.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |