Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |