Home

Context Navigation

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

Last change on this file since 1365 was 1365, checked in by boeske, 10 years ago
large scale forcing enabled
Property svn:keywords set to `Id`
File size: 16.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |