Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |