Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |