Home

Context Navigation

source: palm/trunk/SOURCE/calc_radiation.f90 @ 3186

Last change on this file since 3186 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 14.1 KB

Rev	Line
[1873]	1	!> @file calc_radiation.f90
[2000]	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[1036]	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.
[1036]	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	!
[2718]	17	! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1354]	22	!
[2233]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: calc_radiation.f90 2718 2018-01-02 08:49:38Z suehring $
[2716]	27	! Corrected "Former revisions" section
	28	!
	29	! 2696 2017-12-14 17:12:51Z kanani
	30	! Change in file header (GPL part)
[1321]	31	!
[2716]	32	! 2233 2017-05-30 18:08:54Z suehring
	33	!
[2233]	34	! 2232 2017-05-30 17:47:52Z suehring
	35	! Adjustments to new topography concept
	36	!
[2001]	37	! 2000 2016-08-20 18:09:15Z knoop
	38	! Forced header and separation lines into 80 columns
	39	!
[1874]	40	! 1873 2016-04-18 14:50:06Z maronga
	41	! Module renamed (removed _mod)
	42	!
	43	!
[1851]	44	! 1850 2016-04-08 13:29:27Z maronga
	45	! Module renamed
	46	!
	47	!
[1683]	48	! 1682 2015-10-07 23:56:08Z knoop
	49	! Code annotations made doxygen readable
	50	!
[1354]	51	! 1353 2014-04-08 15:21:23Z heinze
	52	! REAL constants provided with KIND-attribute
	53	!
[1323]	54	! 1322 2014-03-20 16:38:49Z raasch
	55	! exponent 4.0 changed to integer
	56	!
[1321]	57	! 1320 2014-03-20 08:40:49Z raasch
[1320]	58	! ONLY-attribute added to USE-statements,
	59	! kind-parameters added to all INTEGER and REAL declaration statements,
	60	! kinds are defined in new module kinds,
	61	! revision history before 2012 removed,
	62	! comment fields (!:) to be used for variable explanations added to
	63	! all variable declaration statements
[1]	64	!
[1037]	65	! 1036 2012-10-22 13:43:42Z raasch
	66	! code put under GPL (PALM 3.9)
	67	!
[1]	68	! Revision 1.1 2000/04/13 14:42:45 schroeter
	69	! Initial revision
	70	!
	71	!
	72	! Description:
	73	! -------------
[1682]	74	!> Calculation of the vertical divergences of the long-wave radiation-fluxes
	75	!> based on the parameterization of the cloud effective emissivity
[1]	76	!------------------------------------------------------------------------------!
[1682]	77	MODULE calc_radiation_mod
	78
[1320]	79	USE kinds
	80
[1]	81	PRIVATE
	82	PUBLIC calc_radiation
	83
[1682]	84	LOGICAL, SAVE :: first_call = .TRUE. !<
	85	REAL(wp), SAVE :: sigma = 5.67E-08_wp !<
[1]	86
[1682]	87	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: lwp_ground !<
	88	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: lwp_top !<
	89	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: blackbody_emission !<
[1]	90
	91	INTERFACE calc_radiation
	92	MODULE PROCEDURE calc_radiation
	93	MODULE PROCEDURE calc_radiation_ij
	94	END INTERFACE calc_radiation
	95
	96	CONTAINS
	97
	98
	99	!------------------------------------------------------------------------------!
[1682]	100	! Description:
	101	! ------------
	102	!> Call for all grid points
[1]	103	!------------------------------------------------------------------------------!
	104	SUBROUTINE calc_radiation
	105
[1320]	106	USE arrays_3d, &
	107	ONLY: dzw, pt, ql, tend
	108
	109	USE cloud_parameters, &
	110	ONLY: cp, l_d_cp, pt_d_t, t_d_pt
	111
	112	USE control_parameters, &
	113	ONLY: rho_surface
	114
	115	USE indices, &
[2232]	116	ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
[1320]	117
	118	USE kinds
	119
[1]	120	USE pegrid
	121
[1320]	122
[1]	123	IMPLICIT NONE
	124
[1682]	125	INTEGER(iwp) :: i !<
	126	INTEGER(iwp) :: j !<
	127	INTEGER(iwp) :: k !<
	128	INTEGER(iwp) :: k_help !<
[1]	129
[1682]	130	REAL(wp) :: df_p !<
	131	REAL(wp) :: df_m !<
	132	REAL(wp) :: effective_emission_up_m !<
	133	REAL(wp) :: effective_emission_up_p !<
	134	REAL(wp) :: effective_emission_down_m !<
	135	REAL(wp) :: effective_emission_down_p !<
	136	REAL(wp) :: f_up_m !<
	137	REAL(wp) :: f_up_p !<
	138	REAL(wp) :: f_down_m !<
	139	REAL(wp) :: f_down_p !<
	140	REAL(wp) :: impinging_flux_at_top !<
	141	REAL(wp) :: temperature !<
[1]	142
	143
	144	!
	145	!-- On first call, allocate temporary arrays
	146	IF ( first_call ) THEN
[1320]	147	ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1), &
[1]	148	lwp_top(nzb:nzt+1) )
	149	first_call = .FALSE.
	150	ENDIF
	151
	152
	153	DO i = nxl, nxr
	154	DO j = nys, nyn
	155	!
	156	!-- Compute the liquid water path (LWP) and blackbody_emission
	157	!-- at all vertical levels
[1353]	158	lwp_ground(nzb) = 0.0_wp
[1]	159	lwp_top(nzt+1) = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
	160
	161	temperature = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
[1322]	162	blackbody_emission(nzb) = sigma * temperature**4
[1]	163
[2232]	164	DO k = nzb+1, nzt
[1]	165
	166	k_help = ( nzt+nzb+1 ) - k
[1320]	167	lwp_ground(k) = lwp_ground(k-1) + rho_surface * ql(k,j,i) * &
[1]	168	dzw(k)
	169
[1320]	170	lwp_top(k_help) = lwp_top(k_help+1) + &
[1]	171	rho_surface * ql(k_help,j,i) * dzw(k_help)
	172
	173	temperature = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
[2232]	174	blackbody_emission(k) = sigma * temperature**4 &
	175	* MERGE( 1.0_wp, 0.0_wp, &
	176	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	177
	178	ENDDO
	179
[1320]	180	lwp_ground(nzt+1) = lwp_ground(nzt) + &
[1]	181	rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
	182	lwp_top(nzb) = lwp_top(nzb+1)
	183
[1320]	184	temperature = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp * &
[1]	185	ql(nzt+1,j,i)
[1322]	186	blackbody_emission(nzt+1) = sigma * temperature**4
[1]	187
	188	!
	189	!-- See Chlond '92, this is just a first guess
[1353]	190	impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
[1]	191
[2232]	192	DO k = nzb+1, nzt
[1]	193	!
	194	!-- Save some computational time, but this may cause load
	195	!-- imbalances if ql is not distributed uniformly
[1353]	196	IF ( ql(k,j,i) /= 0.0_wp ) THEN
[1]	197	!
	198	!-- Compute effective emissivities
[1353]	199	effective_emission_up_p = 1.0_wp - &
	200	EXP( -130.0_wp * lwp_ground(k+1) )
	201	effective_emission_up_m = 1.0_wp - &
	202	EXP( -130.0_wp * lwp_ground(k-1) )
	203	effective_emission_down_p = 1.0_wp - &
	204	EXP( -158.0_wp * lwp_top(k+1) )
	205	effective_emission_down_m = 1.0_wp - &
	206	EXP( -158.0_wp * lwp_top(k-1) )
[1]	207
	208	!
	209	!-- Compute vertical long wave radiation fluxes
[1320]	210	f_up_p = blackbody_emission(nzb) + &
	211	effective_emission_up_p * &
[1]	212	( blackbody_emission(k) - blackbody_emission(nzb) )
	213
[1320]	214	f_up_m = blackbody_emission(nzb) + &
	215	effective_emission_up_m * &
[1]	216	( blackbody_emission(k-1) - blackbody_emission(nzb) )
	217
[1320]	218	f_down_p = impinging_flux_at_top + &
	219	effective_emission_down_p * &
[1]	220	( blackbody_emission(k) - impinging_flux_at_top )
	221
[1320]	222	f_down_m = impinging_flux_at_top + &
	223	effective_emission_down_m * &
[1]	224	( blackbody_emission(k-1) - impinging_flux_at_top )
	225
	226	!
	227	!-- Divergence of vertical long wave radiation fluxes
	228	df_p = f_up_p - f_down_p
	229	df_m = f_up_m - f_down_m
	230
	231	!
	232	!-- Compute tendency term
[1320]	233	tend(k,j,i) = tend(k,j,i) - &
	234	( pt_d_t(k) / ( rho_surface * cp ) * &
[2232]	235	( df_p - df_m ) / dzw(k) ) &
	236	* MERGE( 1.0_wp, 0.0_wp, &
	237	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	238
	239	ENDIF
	240
	241	ENDDO
	242	ENDDO
	243	ENDDO
	244
	245	END SUBROUTINE calc_radiation
	246
	247
	248	!------------------------------------------------------------------------------!
[1682]	249	! Description:
	250	! ------------
	251	!> Call for grid point i,j
[1]	252	!------------------------------------------------------------------------------!
	253	SUBROUTINE calc_radiation_ij( i, j )
	254
[1320]	255	USE arrays_3d, &
	256	ONLY: dzw, pt, ql, tend
	257
	258	USE cloud_parameters, &
	259	ONLY: cp, l_d_cp, pt_d_t, t_d_pt
	260
	261	USE control_parameters, &
	262	ONLY: rho_surface
	263
	264	USE indices, &
[2232]	265	ONLY: nzb, nzt, wall_flags_0
[1320]	266
	267	USE kinds
	268
[1]	269	USE pegrid
[1320]	270
[1]	271
	272	IMPLICIT NONE
	273
[1682]	274	INTEGER(iwp) :: i !<
	275	INTEGER(iwp) :: j !<
	276	INTEGER(iwp) :: k !<
	277	INTEGER(iwp) :: k_help !<
[1]	278
[1682]	279	REAL(wp) :: df_p !<
	280	REAL(wp) :: df_m !<
	281	REAL(wp) :: effective_emission_up_m !<
	282	REAL(wp) :: effective_emission_up_p !<
	283	REAL(wp) :: effective_emission_down_m !<
	284	REAL(wp) :: effective_emission_down_p !<
	285	REAL(wp) :: f_up_m !<
	286	REAL(wp) :: f_up_p !<
	287	REAL(wp) :: f_down_m !<
	288	REAL(wp) :: f_down_p !<
	289	REAL(wp) :: impinging_flux_at_top !<
	290	REAL(wp) :: temperature !<
[1]	291
[1320]	292
[1]	293	!
	294	!-- On first call, allocate temporary arrays
	295	IF ( first_call ) THEN
[1320]	296	ALLOCATE( blackbody_emission(nzb:nzt+1), lwp_ground(nzb:nzt+1), &
[1]	297	lwp_top(nzb:nzt+1) )
	298	first_call = .FALSE.
	299	ENDIF
	300
	301	!
	302	!-- Compute the liquid water path (LWP) and blackbody_emission
	303	!-- at all vertical levels
[1353]	304	lwp_ground(nzb) = 0.0_wp
[1]	305	lwp_top(nzt+1) = rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
	306
	307	temperature = pt(nzb,j,i) * t_d_pt(nzb) + l_d_cp * ql(nzb,j,i)
[1322]	308	blackbody_emission(nzb) = sigma * temperature**4
[1]	309
[2232]	310	DO k = nzb+1, nzt
[1]	311	k_help = ( nzt+nzb+1 ) - k
	312	lwp_ground(k) = lwp_ground(k-1) + rho_surface * ql(k,j,i) * dzw(k)
	313
[1320]	314	lwp_top(k_help) = lwp_top(k_help+1) + &
[1]	315	rho_surface * ql(k_help,j,i) * dzw(k_help)
	316
	317	temperature = pt(k,j,i) * t_d_pt(k) + l_d_cp * ql(k,j,i)
[2232]	318	blackbody_emission(k) = sigma * temperature**4 &
	319	* MERGE( 1.0_wp, 0.0_wp, &
	320	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	321
	322	ENDDO
[1320]	323	lwp_ground(nzt+1) = lwp_ground(nzt) + &
[1]	324	rho_surface * ql(nzt+1,j,i) * dzw(nzt+1)
	325	lwp_top(nzb) = lwp_top(nzb+1)
	326
[1320]	327	temperature = pt(nzt+1,j,i) * t_d_pt(nzt+1) + l_d_cp * &
[1]	328	ql(nzt+1,j,i)
[1322]	329	blackbody_emission(nzt+1) = sigma * temperature**4
[1]	330
	331	!
	332	!-- See Chlond '92, this is just a first guess
[1353]	333	impinging_flux_at_top = blackbody_emission(nzb) - 100.0_wp
[1]	334
[2232]	335	DO k = nzb+1, nzt
[1]	336	!
	337	!-- Store some computational time,
	338	!-- this may cause load imbalances if ql is not distributed uniformly
[1353]	339	IF ( ql(k,j,i) /= 0.0_wp ) THEN
[1]	340	!
	341	!-- Compute effective emissivities
[1353]	342	effective_emission_up_p = 1.0_wp - &
	343	EXP( -130.0_wp * lwp_ground(k+1) )
	344	effective_emission_up_m = 1.0_wp - &
	345	EXP( -130.0_wp * lwp_ground(k-1) )
	346	effective_emission_down_p = 1.0_wp - &
	347	EXP( -158.0_wp * lwp_top(k+1) )
	348	effective_emission_down_m = 1.0_wp - &
	349	EXP( -158.0_wp * lwp_top(k-1) )
[1]	350
	351	!
	352	!-- Compute vertical long wave radiation fluxes
[1320]	353	f_up_p = blackbody_emission(nzb) + effective_emission_up_p * &
[1]	354	( blackbody_emission(k) - blackbody_emission(nzb) )
	355
[1320]	356	f_up_m = blackbody_emission(nzb) + effective_emission_up_m * &
[1]	357	( blackbody_emission(k-1) - blackbody_emission(nzb) )
	358
[1320]	359	f_down_p = impinging_flux_at_top + effective_emission_down_p * &
[1]	360	( blackbody_emission(k) - impinging_flux_at_top )
	361
[1320]	362	f_down_m = impinging_flux_at_top + effective_emission_down_m * &
[1]	363	( blackbody_emission(k-1) - impinging_flux_at_top )
	364
	365	!
	366	!- Divergence of vertical long wave radiation fluxes
	367	df_p = f_up_p - f_down_p
	368	df_m = f_up_m - f_down_m
	369
	370	!
	371	!-- Compute tendency term
[1320]	372	tend(k,j,i) = tend(k,j,i) - ( pt_d_t(k) / ( rho_surface * cp ) * &
[2232]	373	( df_p - df_m ) / dzw(k) ) &
	374	* MERGE( 1.0_wp, 0.0_wp, &
	375	BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	376
	377	ENDIF
	378
	379	ENDDO
	380
	381	END SUBROUTINE calc_radiation_ij
	382
	383	END MODULE calc_radiation_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |