Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |