Home

Context Navigation

source: palm/trunk/SOURCE/ls_forcing.f90 @ 1382

Last change on this file since 1382 was 1382, checked in by boeske, 10 years ago
minor changes in profile data output of lsf tendencies, variables renamed
Property svn:keywords set to `Id`
File size: 20.8 KB

Line
1	MODULE ls_forcing_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	! Renamed variables which store large scale forcing tendencies
23	! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
24	! q_lsa -> td_lsa_q, q_subs -> td_sub_q,
25	! high\|lowpt_lsa -> high\|low_td_lsa_lpt, ...
26	!
27	! Former revisions:
28	! -----------------
29	! $Id: ls_forcing.f90 1382 2014-04-30 12:15:41Z boeske $
30	!
31	! 1365 2014-04-22 15:03:56Z boeske
32	! Usage of large scale forcing for pt and q enabled:
33	! Added new subroutine ls_advec for horizontal large scale advection and large
34	! scale subsidence,
35	! error message in init_ls_forcing specified,
36	! variable t renamed nt
37	!
38	! 1353 2014-04-08 15:21:23Z heinze
39	! REAL constants provided with KIND-attribute
40	!
41	! 1320 2014-03-20 08:40:49Z raasch
42	! ONLY-attribute added to USE-statements,
43	! kind-parameters added to all INTEGER and REAL declaration statements,
44	! kinds are defined in new module kinds,
45	! comment fields (!:) to be used for variable explanations added to
46	! all variable declaration statements
47	!
48	! 1318 2014-03-17 13:35:16Z raasch
49	! module interfaces removed
50	!
51	! 1299 2014-03-06 13:15:21Z heinze
52	! Ensure a zero large scale vertical velocity at the surface
53	! Bugfix: typo in case of boundary condition in if-clause
54	!
55	! 1276 2014-01-15 13:40:41Z heinze
56	! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
57	!
58	! 1249 2013-11-06 10:45:47Z heinze
59	! remove call of user module
60	! reformatting
61	!
62	! 1241 2013-10-30 11:36:58Z heinze
63	! Initial revision
64	!
65	! Description:
66	! ------------
67	! Calculates large scale forcings (geostrophic wind and subsidence velocity) as
68	! well as surfaces fluxes dependent on time given in an external file (LSF_DATA).
69	! Code is based in parts on DALES and UCLA-LES.
70	!--------------------------------------------------------------------------------!
71
72	PRIVATE
73	PUBLIC init_ls_forcing, ls_forcing_surf, ls_forcing_vert, ls_advec
74	SAVE
75
76	INTERFACE ls_advec
77	MODULE PROCEDURE ls_advec
78	MODULE PROCEDURE ls_advec_ij
79	END INTERFACE ls_advec
80
81	CONTAINS
82
83	SUBROUTINE init_ls_forcing
84
85	USE arrays_3d, &
86	ONLY: p_surf, pt_surf, q_surf, qsws_surf, shf_surf, td_lsa_lpt, &
87	td_lsa_q, td_sub_lpt, td_sub_q, time_surf, time_vert, &
88	ug_vert, vg_vert, wsubs_vert, zu
89
90	USE control_parameters, &
91	ONLY: end_time, lsf_surf, lsf_vert, message_string, nlsf
92
93	USE indices, &
94	ONLY: ngp_sums_ls, nzb, nz, nzt
95
96	USE kinds
97
98	USE statistics, &
99	ONLY: sums_ls_l
100
101
102	IMPLICIT NONE
103
104	CHARACTER(100) :: chmess !:
105	CHARACTER(1) :: hash !:
106
107	INTEGER(iwp) :: ierrn !:
108	INTEGER(iwp) :: finput = 90 !:
109	INTEGER(iwp) :: k !:
110	INTEGER(iwp) :: nt !:
111
112	REAL(wp) :: fac !:
113	REAL(wp) :: highheight !:
114	REAL(wp) :: highug_vert !:
115	REAL(wp) :: highvg_vert !:
116	REAL(wp) :: highwsubs_vert !:
117	REAL(wp) :: lowheight !:
118	REAL(wp) :: lowug_vert !:
119	REAL(wp) :: lowvg_vert !:
120	REAL(wp) :: lowwsubs_vert !:
121	REAL(wp) :: high_td_lsa_lpt !:
122	REAL(wp) :: low_td_lsa_lpt !:
123	REAL(wp) :: high_td_lsa_q !:
124	REAL(wp) :: low_td_lsa_q !:
125	REAL(wp) :: high_td_sub_lpt !:
126	REAL(wp) :: low_td_sub_lpt !:
127	REAL(wp) :: high_td_sub_q !:
128	REAL(wp) :: low_td_sub_q !:
129	REAL(wp) :: r_dummy !:
130
131	ALLOCATE( p_surf(0:nlsf), pt_surf(0:nlsf), q_surf(0:nlsf), &
132	qsws_surf(0:nlsf), shf_surf(0:nlsf), &
133	td_lsa_lpt(nzb:nzt+1,0:nlsf), td_lsa_q(nzb:nzt+1,0:nlsf), &
134	td_sub_lpt(nzb:nzt+1,0:nlsf), td_sub_q(nzb:nzt+1,0:nlsf), &
135	time_vert(0:nlsf), time_surf(0:nlsf), &
136	ug_vert(nzb:nzt+1,0:nlsf), vg_vert(nzb:nzt+1,0:nlsf), &
137	wsubs_vert(nzb:nzt+1,0:nlsf) )
138
139	p_surf = 0.0_wp; pt_surf = 0.0_wp; q_surf = 0.0_wp; qsws_surf = 0.0_wp
140	shf_surf = 0.0_wp; time_vert = 0.0_wp; td_lsa_lpt = 0.0_wp
141	td_lsa_q = 0.0_wp; td_sub_lpt = 0.0_wp; td_sub_q = 0.0_wp
142	time_surf = 0.0_wp; ug_vert = 0.0_wp; vg_vert = 0.0_wp
143	wsubs_vert = 0.0_wp
144
145	!
146	!-- Array for storing large scale forcing and nudging tendencies at each
147	!-- timestep for data output
148	ALLOCATE( sums_ls_l(nzb:nzt+1,0:7) )
149	sums_ls_l = 0.0_wp
150
151	ngp_sums_ls = (nz+2)*6
152
153	OPEN ( finput, FILE='LSF_DATA', STATUS='OLD', &
154	FORM='FORMATTED', IOSTAT=ierrn )
155
156	IF ( ierrn /= 0 ) THEN
157	message_string = 'file LSF_DATA does not exist'
158	CALL message( 'ls_forcing', 'PA0368', 1, 2, 0, 6, 0 )
159	ENDIF
160
161	ierrn = 0
162	!
163	!-- First three lines of LSF_DATA contain header
164	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
165	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
166	READ ( finput, FMT='(a100)', IOSTAT=ierrn ) chmess
167
168	IF ( ierrn /= 0 ) THEN
169	message_string = 'errors in file LSF_DATA'
170	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
171	ENDIF
172
173	!
174	!-- Surface values are read in
175	nt = 0
176	ierrn = 0
177
178	DO WHILE ( time_surf(nt) < end_time )
179	nt = nt + 1
180	READ ( finput, *, IOSTAT = ierrn ) time_surf(nt), shf_surf(nt), &
181	qsws_surf(nt), pt_surf(nt), &
182	q_surf(nt), p_surf(nt)
183
184	IF ( ierrn < 0 ) THEN
185	WRITE ( message_string, * ) 'No time dependent surface variables ',&
186	'in&LSF_DATA for end of run found'
187
188	CALL message( 'ls_forcing', 'PA0370', 1, 2, 0, 6, 0 )
189	ENDIF
190	ENDDO
191
192
193	IF ( time_surf(1) > end_time ) THEN
194	WRITE ( message_string, * ) 'No time dependent surface variables in ',&
195	'&LSF_DATA for end of run found - ', &
196	'lsf_surf is set to FALSE'
197	CALL message( 'ls_forcing', 'PA0371', 0, 0, 0, 6, 0 )
198	lsf_surf = .FALSE.
199	ENDIF
200
201	!
202	!-- Go to the end of the list with surface variables
203	DO WHILE ( ierrn == 0 )
204	READ ( finput, *, IOSTAT = ierrn ) r_dummy
205	ENDDO
206
207	!
208	!-- Profiles of ug, vg and w_subs are read in (large scale forcing)
209
210	nt = 0
211	DO WHILE ( time_vert(nt) < end_time )
212	nt = nt + 1
213	hash = "#"
214	ierrn = 1 ! not zero
215	!
216	!-- Search for the next line consisting of "# time",
217	!-- from there onwards the profiles will be read
218	DO WHILE ( .NOT. ( hash == "#" .AND. ierrn == 0 ) )
219	READ ( finput, *, IOSTAT=ierrn ) hash, time_vert(nt)
220	IF ( ierrn < 0 ) THEN
221	WRITE( message_string, * ) 'No time dependent vertical profiles',&
222	' in&LSF_DATA for end of run found'
223	CALL message( 'ls_forcing', 'PA0372', 1, 2, 0, 6, 0 )
224	ENDIF
225	ENDDO
226
227	IF ( nt == 1 .AND. time_vert(nt) > end_time ) EXIT
228
229	READ ( finput, *, IOSTAT=ierrn ) lowheight, lowug_vert, lowvg_vert, &
230	lowwsubs_vert, low_td_lsa_lpt, &
231	low_td_lsa_q, low_td_sub_lpt, &
232	low_td_sub_q
233	IF ( ierrn /= 0 ) THEN
234	message_string = 'errors in file LSF_DATA'
235	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
236	ENDIF
237
238	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
239	highvg_vert, highwsubs_vert, &
240	high_td_lsa_lpt, high_td_lsa_q, &
241	high_td_sub_lpt, high_td_sub_q
242
243	IF ( ierrn /= 0 ) THEN
244	message_string = 'errors in file LSF_DATA'
245	CALL message( 'ls_forcing', 'PA0369', 1, 2, 0, 6, 0 )
246	ENDIF
247
248
249	DO k = nzb, nzt+1
250	IF ( highheight < zu(k) ) THEN
251	lowheight = highheight
252	lowug_vert = highug_vert
253	lowvg_vert = highvg_vert
254	lowwsubs_vert = highwsubs_vert
255	low_td_lsa_lpt = high_td_lsa_lpt
256	low_td_lsa_q = high_td_lsa_q
257	low_td_sub_lpt = high_td_sub_lpt
258	low_td_sub_q = high_td_sub_q
259
260	ierrn = 0
261	READ ( finput, *, IOSTAT=ierrn ) highheight, highug_vert, &
262	highvg_vert, highwsubs_vert, &
263	high_td_lsa_lpt, &
264	high_td_lsa_q, &
265	high_td_sub_lpt, high_td_sub_q
266
267	IF ( ierrn /= 0 ) THEN
268	WRITE( message_string, * ) 'zu(nzt+1) = ', zu(nzt+1), 'm ', &
269	'is higher than the maximum height in LSF_DATA which ',&
270	'is ', lowheight, 'm. Interpolation on PALM ', &
271	'grid is not possible.'
272	CALL message( 'ls_forcing', 'PA0395', 1, 2, 0, 6, 0 )
273	ENDIF
274
275	ENDIF
276
277	!
278	!-- Interpolation of prescribed profiles in space
279	fac = (highheight-zu(k))/(highheight - lowheight)
280
281	ug_vert(k,nt) = fac * lowug_vert &
282	+ ( 1.0_wp - fac ) * highug_vert
283	vg_vert(k,nt) = fac * lowvg_vert &
284	+ ( 1.0_wp - fac ) * highvg_vert
285	wsubs_vert(k,nt) = fac * lowwsubs_vert &
286	+ ( 1.0_wp - fac ) * highwsubs_vert
287
288	td_lsa_lpt(k,nt) = fac * low_td_lsa_lpt &
289	+ ( 1.0_wp - fac ) * high_td_lsa_lpt
290	td_lsa_q(k,nt) = fac * low_td_lsa_q &
291	+ ( 1.0_wp - fac ) * high_td_lsa_q
292	td_sub_lpt(k,nt) = fac * low_td_sub_lpt &
293	+ ( 1.0_wp - fac ) * high_td_sub_lpt
294	td_sub_q(k,nt) = fac * low_td_sub_q &
295	+ ( 1.0_wp - fac ) * high_td_sub_q
296
297	ENDDO
298
299	ENDDO
300
301	!
302	!-- Large scale vertical velocity has to be zero at the surface
303	wsubs_vert(nzb,:) = 0.0_wp
304
305	IF ( time_vert(1) > end_time ) THEN
306	WRITE ( message_string, * ) 'Time dependent large scale profile ', &
307	'forcing from&LSF_DATA sets in after end of ' , &
308	'simulation - lsf_vert is set to FALSE'
309	CALL message( 'ls_forcing', 'PA0373', 0, 0, 0, 6, 0 )
310	lsf_vert = .FALSE.
311	ENDIF
312
313	CLOSE( finput )
314
315
316	END SUBROUTINE init_ls_forcing
317
318
319	SUBROUTINE ls_forcing_surf ( time )
320
321	USE arrays_3d, &
322	ONLY: p_surf, pt_surf, q_surf, qsws, qsws_surf, shf, shf_surf, &
323	time_surf, time_vert, ug, ug_vert, vg, vg_vert
324
325	USE control_parameters, &
326	ONLY: bc_q_b, ibc_pt_b, ibc_q_b, pt_surface, q_surface, &
327	surface_pressure
328
329	USE kinds
330
331	IMPLICIT NONE
332
333	INTEGER(iwp) :: nt !:
334
335	REAL(wp) :: fac !:
336	REAL(wp), INTENT(in) :: time !:
337
338	!
339	!-- Interpolation in time of LSF_DATA at the surface
340	nt = 1
341	DO WHILE ( time > time_surf(nt) )
342	nt = nt + 1
343	ENDDO
344	IF ( time /= time_surf(nt) ) THEN
345	nt = nt - 1
346	ENDIF
347
348	fac = ( time -time_surf(nt) ) / ( time_surf(nt+1) - time_surf(nt) )
349
350	IF ( ibc_pt_b == 0 ) THEN
351	!
352	!-- In case of Dirichlet boundary condition shf must not
353	!-- be set - it is calculated via MOST in prandtl_fluxes
354	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
355
356	ELSEIF ( ibc_pt_b == 1 ) THEN
357	!
358	!-- In case of Neumann boundary condition pt_surface is needed for
359	!-- calculation of reference density
360	shf = shf_surf(nt) + fac * ( shf_surf(nt+1) - shf_surf(nt) )
361	pt_surface = pt_surf(nt) + fac * ( pt_surf(nt+1) - pt_surf(nt) )
362
363	ENDIF
364
365	IF ( ibc_q_b == 0 ) THEN
366	!
367	!-- In case of Dirichlet boundary condition qsws must not
368	!-- be set - it is calculated via MOST in prandtl_fluxes
369	q_surface = q_surf(nt) + fac * ( q_surf(nt+1) - q_surf(nt) )
370
371	ELSEIF ( ibc_q_b == 1 ) THEN
372
373	qsws = qsws_surf(nt) + fac * ( qsws_surf(nt+1) - qsws_surf(nt) )
374
375	ENDIF
376
377	surface_pressure = p_surf(nt) + fac * ( p_surf(nt+1) - p_surf(nt) )
378
379	END SUBROUTINE ls_forcing_surf
380
381
382	SUBROUTINE ls_forcing_vert ( time )
383
384	USE arrays_3d, &
385	ONLY: time_vert, ug, ug_vert, vg, vg_vert, w_subs, wsubs_vert
386
387	USE control_parameters, &
388	ONLY: large_scale_subsidence
389
390	USE kinds
391
392	IMPLICIT NONE
393
394	INTEGER(iwp) :: nt !:
395
396	REAL(wp) :: fac !:
397	REAL(wp), INTENT(in) :: time !:
398
399	!
400	!-- Interpolation in time of LSF_DATA for ug, vg and w_subs
401	nt = 1
402	DO WHILE ( time > time_vert(nt) )
403	nt = nt + 1
404	ENDDO
405	IF ( time /= time_vert(nt) ) THEN
406	nt = nt - 1
407	ENDIF
408
409	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
410
411	ug = ug_vert(:,nt) + fac * ( ug_vert(:,nt+1) - ug_vert(:,nt) )
412	vg = vg_vert(:,nt) + fac * ( vg_vert(:,nt+1) - vg_vert(:,nt) )
413
414	IF ( large_scale_subsidence ) THEN
415	w_subs = wsubs_vert(:,nt) &
416	+ fac * ( wsubs_vert(:,nt+1) - wsubs_vert(:,nt) )
417	ENDIF
418
419	END SUBROUTINE ls_forcing_vert
420
421
422	!------------------------------------------------------------------------------!
423	! Call for all grid points
424	!------------------------------------------------------------------------------!
425	SUBROUTINE ls_advec ( time, prog_var )
426
427	USE arrays_3d, &
428	ONLY: td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, tend, time_vert
429
430	USE control_parameters, &
431	ONLY: large_scale_subsidence, use_subsidence_tendencies
432
433	USE indices
434
435	USE kinds
436
437	IMPLICIT NONE
438
439	CHARACTER (LEN=*) :: prog_var !:
440
441	REAL(wp), INTENT(in) :: time !:
442	REAL(wp) :: fac !:
443
444	INTEGER(iwp) :: i !:
445	INTEGER(iwp) :: j !:
446	INTEGER(iwp) :: k !:
447	INTEGER(iwp) :: nt !:
448
449	!
450	!-- Interpolation in time of LSF_DATA
451	nt = 1
452	DO WHILE ( time > time_vert(nt) )
453	nt = nt + 1
454	ENDDO
455	IF ( time /= time_vert(nt) ) THEN
456	nt = nt - 1
457	ENDIF
458
459	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
460
461	!
462	!-- Add horizontal large scale advection tendencies of pt and q
463	SELECT CASE ( prog_var )
464
465	CASE ( 'pt' )
466
467	DO i = nxl, nxr
468	DO j = nys, nyn
469	DO k = nzb_u_inner(j,i)+1, nzt
470	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) + fac * &
471	( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
472	ENDDO
473	ENDDO
474	ENDDO
475
476	CASE ( 'q' )
477
478	DO i = nxl, nxr
479	DO j = nys, nyn
480	DO k = nzb_u_inner(j,i)+1, nzt
481	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) + fac * &
482	( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
483	ENDDO
484	ENDDO
485	ENDDO
486
487	END SELECT
488
489	!
490	!-- Subsidence of pt and q with prescribed subsidence tendencies
491	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
492
493	SELECT CASE ( prog_var )
494
495	CASE ( 'pt' )
496
497	DO i = nxl, nxr
498	DO j = nys, nyn
499	DO k = nzb_u_inner(j,i)+1, nzt
500	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
501	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
502	ENDDO
503	ENDDO
504	ENDDO
505
506	CASE ( 'q' )
507
508	DO i = nxl, nxr
509	DO j = nys, nyn
510	DO k = nzb_u_inner(j,i)+1, nzt
511	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
512	( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
513	ENDDO
514	ENDDO
515	ENDDO
516
517	END SELECT
518
519	ENDIF
520
521	END SUBROUTINE ls_advec
522
523
524	!------------------------------------------------------------------------------!
525	! Call for grid point i,j
526	!------------------------------------------------------------------------------!
527	SUBROUTINE ls_advec_ij ( i, j, time, prog_var )
528
529	USE arrays_3d, &
530	ONLY: td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, tend, time_vert
531
532	USE control_parameters, &
533	ONLY: large_scale_subsidence, use_subsidence_tendencies
534
535	USE indices
536
537	USE kinds
538
539	IMPLICIT NONE
540
541	CHARACTER (LEN=*) :: prog_var !:
542
543	REAL(wp), INTENT(in) :: time !:
544	REAL(wp) :: fac !:
545
546	INTEGER(iwp) :: i !:
547	INTEGER(iwp) :: j !:
548	INTEGER(iwp) :: k !:
549	INTEGER(iwp) :: nt !:
550
551	!
552	!-- Interpolation in time of LSF_DATA
553	nt = 1
554	DO WHILE ( time > time_vert(nt) )
555	nt = nt + 1
556	ENDDO
557	IF ( time /= time_vert(nt) ) THEN
558	nt = nt - 1
559	ENDIF
560
561	fac = ( time-time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
562
563	!
564	!-- Add horizontal large scale advection tendencies of pt and q
565	SELECT CASE ( prog_var )
566
567	CASE ( 'pt' )
568
569	DO k = nzb_u_inner(j,i)+1, nzt
570	tend(k,j,i) = tend(k,j,i) + td_lsa_lpt(k,nt) &
571	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
572	ENDDO
573
574	CASE ( 'q' )
575
576	DO k = nzb_u_inner(j,i)+1, nzt
577	tend(k,j,i) = tend(k,j,i) + td_lsa_q(k,nt) &
578	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
579	ENDDO
580
581	END SELECT
582
583	!
584	!-- Subsidence of pt and q with prescribed profiles
585	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
586
587	SELECT CASE ( prog_var )
588
589	CASE ( 'pt' )
590
591	DO k = nzb_u_inner(j,i)+1, nzt
592	tend(k,j,i) = tend(k,j,i) + td_sub_lpt(k,nt) + fac * &
593	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
594	ENDDO
595
596	CASE ( 'q' )
597
598	DO k = nzb_u_inner(j,i)+1, nzt
599	tend(k,j,i) = tend(k,j,i) + td_sub_q(k,nt) + fac * &
600	( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
601	ENDDO
602
603	END SELECT
604
605	ENDIF
606
607	END SUBROUTINE ls_advec_ij
608
609
610	END MODULE ls_forcing_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |