Home

Context Navigation

source: palm/trunk/SOURCE/nudging_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 9 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 19.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |