Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |