Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |