Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |