Home

Context Navigation

source: palm/trunk/SOURCE/ls_forcing_mod.f90 @ 2310

Last change on this file since 2310 was 2271, checked in by sward, 7 years ago
error messages and numbers updated
Property svn:keywords set to `Id`
File size: 25.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |