Home

Context Navigation

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

Last change on this file since 1587 was 1354, checked in by heinze, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 12.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |