Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

radiation_tripleclouds_lw.F90

Last change on this file was 5185, checked in by abarral, 2 months ago
Replace REPROBUS CPP KEY by logical using handmade wonky wrapper
File size: 25.0 KB

Rev	Line
[3908]	1	! radiation_tripleclouds_lw.F90 - Longwave "Tripleclouds" solver
[5159]	2
[3908]	3	! (C) Copyright 2016- ECMWF.
[5159]	4
[3908]	5	! This software is licensed under the terms of the Apache Licence Version 2.0
	6	! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
[5159]	7
[3908]	8	! In applying this licence, ECMWF does not waive the privileges and immunities
	9	! granted to it by virtue of its status as an intergovernmental organisation
	10	! nor does it submit to any jurisdiction.
[5159]	11
[3908]	12	! Author: Robin Hogan
	13	! Email: r.j.hogan@ecmwf.int
[5159]	14
[3908]	15	! Modifications
	16	! 2017-04-28 R. Hogan Receive emission/albedo rather than planck/emissivity
	17	! 2017-04-22 R. Hogan Store surface fluxes at all g-points
	18	! 2017-10-23 R. Hogan Renamed single-character variables
	19	! 2018-10-08 R. Hogan Call calc_region_properties
[4489]	20	! 2020-09-18 R. Hogan Replaced some array expressions with loops
	21	! 2020-09-19 R. Hogan Implement the cloud-only-scattering optimization
[3908]	22
	23	module radiation_tripleclouds_lw
	24
	25	public
	26
	27	contains
	28	! Small routine for scaling cloud optical depth in the cloudy
	29	! regions
	30	#include "radiation_optical_depth_scaling.h"
	31
[4489]	32	!---------------------------------------------------------------------
[3908]	33	! This module contains just one subroutine, the longwave
	34	! "Tripleclouds" solver in which cloud inhomogeneity is treated by
	35	! dividing each model level into three regions, one clear and two
	36	! cloudy (with differing optical depth). This approach was described
	37	! by Shonk and Hogan (2008).
	38	subroutine solver_tripleclouds_lw(nlev,istartcol,iendcol, &
	39	& config, cloud, &
	40	& od, ssa, g, od_cloud, ssa_cloud, g_cloud, planck_hl, &
	41	& emission, albedo, &
	42	& flux)
	43
	44	use parkind1, only : jprb
	45	use yomhook, only : lhook, dr_hook
	46
	47	! use radiation_io, only : nulout
	48	use radiation_config, only : config_type, IPdfShapeGamma
	49	use radiation_cloud, only : cloud_type
	50	use radiation_regions, only : calc_region_properties
	51	use radiation_overlap, only : calc_overlap_matrices
[4489]	52	use radiation_flux, only : flux_type, indexed_sum
[3908]	53	use radiation_matrix, only : singlemat_x_vec
	54	use radiation_two_stream, only : calc_two_stream_gammas_lw, &
	55	& calc_reflectance_transmittance_lw, &
	56	& calc_no_scattering_transmittance_lw
[4489]	57	use radiation_adding_ica_lw, only : adding_ica_lw, calc_fluxes_no_scattering_lw
[3908]	58	use radiation_lw_derivatives, only : calc_lw_derivatives_region
	59
	60	implicit none
	61
	62	! Inputs
	63	integer, intent(in) :: nlev ! number of model levels
	64	integer, intent(in) :: istartcol, iendcol ! range of columns to process
	65	type(config_type), intent(in) :: config
	66	type(cloud_type), intent(in) :: cloud
	67
	68	! Gas and aerosol optical depth of each layer at each longwave
	69	! g-point
	70	real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od
	71
	72	! Gas and aerosol single-scattering albedo and asymmetry factor,
	73	! only if longwave scattering by aerosols is to be represented
	74	real(jprb), intent(in), &
	75	& dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: ssa, g
	76
	77	! Cloud and precipitation optical depth of each layer in each
	78	! longwave band
	79	real(jprb), intent(in) :: od_cloud(config%n_bands_lw,nlev,istartcol:iendcol)
	80
	81	! Cloud and precipitation single-scattering albedo and asymmetry
	82	! factor, only if longwave scattering by clouds is to be
	83	! represented
	84	real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering, &
	85	& nlev,istartcol:iendcol) :: ssa_cloud, g_cloud
	86
	87	! Planck function (emitted flux from a black body) at half levels
	88	! and at the surface at each longwave g-point
	89	real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: planck_hl
	90
	91	! Emission (Planck*emissivity) and albedo (1-emissivity) at the
	92	! surface at each longwave g-point
	93	real(jprb), intent(in), dimension(config%n_g_lw, istartcol:iendcol) :: emission, albedo
	94
	95	! Optical depth, single scattering albedo and asymmetry factor in
	96	! each g-point including gas, aerosol and clouds
	97	real(jprb), dimension(config%n_g_lw) :: od_total, ssa_total, g_total
	98
	99	! Modified optical depth after Tripleclouds scaling to represent
	100	! cloud inhomogeneity
	101	real(jprb), dimension(config%n_g_lw) :: od_cloud_new
	102
	103	! Output
	104	type(flux_type), intent(inout):: flux
	105
	106	! Local constants
	107	integer, parameter :: nregions = 3
	108
	109	! In a clear-sky layer this will be 1, otherwise equal to nregions
	110	integer :: nreg
	111
	112	! Local variables
	113
	114	! The area fractions of each region
	115	real(jprb) :: region_fracs(1:nregions,nlev,istartcol:iendcol)
	116
	117	! The scaling used for the optical depth in the cloudy regions
	118	real(jprb) :: od_scaling(2:nregions,nlev,istartcol:iendcol)
	119
	120	! Directional overlap matrices defined at all layer interfaces
	121	! including top-of-atmosphere and the surface
	122	real(jprb), dimension(nregions,nregions,nlev+1, &
	123	& istartcol:iendcol) :: u_matrix, v_matrix
	124
	125	! Two-stream variables
	126	real(jprb), dimension(config%n_g_lw) :: gamma1, gamma2
	127
	128	! Diffuse reflection and transmission matrices of each layer
	129	real(jprb), dimension(config%n_g_lw, nregions, nlev) :: reflectance, transmittance
	130
	131	! Emission by a layer into the upwelling or downwelling diffuse
	132	! streams
	133	real(jprb), dimension(config%n_g_lw, nregions, nlev) &
[4489]	134	& :: source_up, source_dn
[3908]	135
[4489]	136	! Clear-sky reflectance and transmittance
	137	real(jprb), dimension(config%n_g_lw, nlev) &
	138	& :: ref_clear, trans_clear
	139
[3908]	140	! ...clear-sky equivalent
	141	real(jprb), dimension(config%n_g_lw, nlev) &
[4489]	142	& :: source_up_clear, source_dn_clear
[3908]	143
	144	! Total albedo of the atmosphere/surface just above a layer
	145	! interface with respect to downwelling diffuse radiation at that
	146	! interface, where level index = 1 corresponds to the
	147	! top-of-atmosphere
	148	real(jprb), dimension(config%n_g_lw, nregions, nlev+1) :: total_albedo
	149
	150	! Upwelling radiation just above a layer interface due to emission
	151	! below that interface, where level index = 1 corresponds to the
	152	! top-of-atmosphere
	153	real(jprb), dimension(config%n_g_lw, nregions, nlev+1) :: total_source
	154
	155	! Total albedo and source of the atmosphere just below a layer interface
	156	real(jprb), dimension(config%n_g_lw, nregions) &
	157	& :: total_albedo_below, total_source_below
	158
	159	! Downwelling flux below and above an interface between
	160	! layers into a plane perpendicular to the direction of the sun
	161	real(jprb), dimension(config%n_g_lw, nregions) &
	162	& :: flux_dn, flux_dn_below, flux_up
	163
	164	! ...clear-sky equivalent (no distinction between "above/below")
[4489]	165	real(jprb), dimension(config%n_g_lw, nlev+1) &
[3908]	166	& :: flux_dn_clear, flux_up_clear
	167
	168	! Clear-sky equivalent, but actually its reciprocal to replace
	169	! some divisions by multiplications
	170	real(jprb), dimension(config%n_g_lw, nregions) :: inv_denom
	171
	172	! Identify clear-sky layers, with pseudo layers for outer space
	173	! and below the ground, both treated as single-region clear skies
	174	logical :: is_clear_sky_layer(0:nlev+1)
	175
[4489]	176	! Index of the highest cloudy layer
	177	integer :: i_cloud_top
	178
[3908]	179	integer :: jcol, jlev, jg, jreg, jreg2, ng
	180
	181	real(jprb) :: hook_handle
	182
	183	if (lhook) call dr_hook('radiation_tripleclouds_lw:solver_tripleclouds_lw',0,hook_handle)
	184
	185	! --------------------------------------------------------
	186	! Section 1: Prepare general variables and arrays
	187	! --------------------------------------------------------
	188	! Copy array dimensions to local variables for convenience
	189	ng = config%n_g_lw
	190
	191	! Compute the wavelength-independent region fractions and
	192	! optical-depth scalings
	193	call calc_region_properties(nlev,nregions,istartcol,iendcol, &
	194	& config%i_cloud_pdf_shape == IPdfShapeGamma, &
	195	& cloud%fraction, cloud%fractional_std, region_fracs, &
	196	& od_scaling, config%cloud_fraction_threshold)
	197
	198	! Compute wavelength-independent overlap matrices u_matrix and v_matrix
	199	call calc_overlap_matrices(nlev,nregions,istartcol,iendcol, &
	200	& region_fracs, cloud%overlap_param, &
	201	& u_matrix, v_matrix, &
	202	& decorrelation_scaling=config%cloud_inhom_decorr_scaling, &
	203	& cloud_fraction_threshold=config%cloud_fraction_threshold, &
	204	& use_beta_overlap=config%use_beta_overlap, &
	205	& cloud_cover=flux%cloud_cover_lw)
	206
	207	! Main loop over columns
[5158]	208	DO jcol = istartcol, iendcol
[3908]	209	! --------------------------------------------------------
	210	! Section 2: Prepare column-specific variables and arrays
	211	! --------------------------------------------------------
	212
	213	! Define which layers contain cloud; assume that
	214	! cloud%crop_cloud_fraction has already been called
	215	is_clear_sky_layer = .true.
[4489]	216	i_cloud_top = nlev+1
[5158]	217	DO jlev = 1,nlev
[3908]	218	if (cloud%fraction(jcol,jlev) > 0.0_jprb) then
	219	is_clear_sky_layer(jlev) = .false.
[4489]	220	! Get index to the first cloudy layer from the top
	221	if (i_cloud_top > jlev) then
	222	i_cloud_top = jlev
	223	end if
[3908]	224	end if
	225	end do
[4489]	226	if (config%do_lw_aerosol_scattering) then
	227	! This is actually the first layer in which we need to
	228	! consider scattering
	229	i_cloud_top = 1
	230	end if
[3908]	231
	232	! --------------------------------------------------------
[4489]	233	! Section 3: Clear-sky calculation
[3908]	234	! --------------------------------------------------------
[4489]	235
	236	if (.not. config%do_lw_aerosol_scattering) then
	237	! No scattering in clear-sky flux calculation
[5158]	238	DO jlev = 1,nlev
[4489]	239	! Array-wise assignments
	240	gamma1 = 0.0_jprb
	241	gamma2 = 0.0_jprb
	242	call calc_no_scattering_transmittance_lw(ng, od(:,jlev,jcol), &
	243	& planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
	244	& trans_clear(:,jlev), source_up_clear(:,jlev), source_dn_clear(:,jlev))
	245	ref_clear(:,jlev) = 0.0_jprb
	246	end do
	247	! Simple down-then-up method to compute fluxes
	248	call calc_fluxes_no_scattering_lw(ng, nlev, &
	249	& trans_clear, source_up_clear, source_dn_clear, &
	250	& emission(:,jcol), albedo(:,jcol), &
	251	& flux_up_clear, flux_dn_clear)
	252	else
	253	! Scattering in clear-sky flux calculation
[5158]	254	DO jlev = 1,nlev
[4489]	255	! Array-wise assignments
	256	gamma1 = 0.0_jprb
	257	gamma2 = 0.0_jprb
	258	call calc_two_stream_gammas_lw(ng, &
	259	& ssa(:,jlev,jcol), g(:,jlev,jcol), gamma1, gamma2)
	260	call calc_reflectance_transmittance_lw(ng, &
	261	& od(:,jlev,jcol), gamma1, gamma2, &
	262	& planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
	263	& ref_clear(:,jlev), trans_clear(:,jlev), &
	264	& source_up_clear(:,jlev), source_dn_clear(:,jlev))
	265	end do
	266	! Use adding method to compute fluxes
	267	call adding_ica_lw(ng, nlev, &
	268	& ref_clear, trans_clear, source_up_clear, source_dn_clear, &
	269	& emission(:,jcol), albedo(:,jcol), &
	270	& flux_up_clear, flux_dn_clear)
	271	end if
	272
	273	if (config%do_clear) then
	274	! Sum over g-points to compute broadband fluxes
	275	flux%lw_up_clear(jcol,:) = sum(flux_up_clear,1)
	276	flux%lw_dn_clear(jcol,:) = sum(flux_dn_clear,1)
	277	! Store surface spectral downwelling fluxes
	278	flux%lw_dn_surf_clear_g(:,jcol) = flux_dn_clear(:,nlev+1)
	279	! Save the spectral fluxes if required
	280	if (config%do_save_spectral_flux) then
[5158]	281	DO jlev = 1,nlev+1
[4489]	282	call indexed_sum(flux_up_clear(:,jlev), &
	283	& config%i_spec_from_reordered_g_lw, &
	284	& flux%lw_up_clear_band(:,jcol,jlev))
	285	call indexed_sum(flux_dn_clear(:,jlev), &
	286	& config%i_spec_from_reordered_g_lw, &
	287	& flux%lw_dn_clear_band(:,jcol,jlev))
	288	end do
	289	end if
	290	end if
	291
	292	! --------------------------------------------------------
	293	! Section 4: Loop over cloudy layers to compute reflectance and transmittance
	294	! --------------------------------------------------------
[3908]	295	! In this section the reflectance, transmittance and sources
	296	! are computed for each layer
[4489]	297
	298	! Firstly, ensure clear-sky transmittance is valid for whole
	299	! depth of the atmosphere, because even above cloud it is used
	300	! by the LW derivatives
	301	transmittance(:,1,:) = trans_clear(:,:)
	302	! Dummy values in cloudy regions above cloud top
	303	if (i_cloud_top > 0) then
	304	transmittance(:,2:,1:min(i_cloud_top,nlev)) = 1.0_jprb
	305	end if
[3908]	306
[5158]	307	DO jlev = i_cloud_top,nlev ! Start at cloud top and work down
[4489]	308
[3908]	309	! Array-wise assignments
	310	gamma1 = 0.0_jprb
	311	gamma2 = 0.0_jprb
	312
[4489]	313	! Copy over clear-sky properties
	314	reflectance(:,1,jlev) = ref_clear(:,jlev)
	315	source_up(:,1,jlev) = source_up_clear(:,jlev) ! Scaled later by region size
	316	source_dn(:,1,jlev) = source_dn_clear(:,jlev) ! Scaled later by region size
[3908]	317	nreg = nregions
	318	if (is_clear_sky_layer(jlev)) then
	319	nreg = 1
	320	reflectance(:,2:,jlev) = 0.0_jprb
[4489]	321	transmittance(:,2:,jlev) = 1.0_jprb
	322	source_up(:,2:,jlev) = 0.0_jprb
	323	source_dn(:,2:,jlev) = 0.0_jprb
	324	else
[5158]	325	DO jreg = 2,nreg
[3908]	326	! Cloudy sky
	327	! Add scaled cloud optical depth to clear-sky value
	328	od_cloud_new = od_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	329	& * od_scaling(jreg,jlev,jcol)
	330	od_total = od(:,jlev,jcol) + od_cloud_new
	331
	332	if (config%do_lw_cloud_scattering) then
	333	ssa_total = 0.0_jprb
	334	g_total = 0.0_jprb
	335	if (config%do_lw_aerosol_scattering) then
	336	where (od_total > 0.0_jprb)
	337	ssa_total = (ssa(:,jlev,jcol)*od(:,jlev,jcol) &
	338	& + ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	339	& * od_cloud_new) &
	340	& / od_total
	341	end where
[5185]	342	where (ssa_total > 0.0_jprb .AND. od_total > 0.0_jprb)
[3908]	343	g_total = (g(:,jlev,jcol)ssa(:,jlev,jcol)od(:,jlev,jcol) &
	344	& + g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	345	& * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	346	& * od_cloud_new) &
	347	& / (ssa_total*od_total)
	348	end where
	349	else
	350	where (od_total > 0.0_jprb)
	351	ssa_total = ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	352	& * od_cloud_new / od_total
	353	end where
[5185]	354	where (ssa_total > 0.0_jprb .AND. od_total > 0.0_jprb)
[3908]	355	g_total = g_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	356	& * ssa_cloud(config%i_band_from_reordered_g_lw,jlev,jcol) &
	357	& * od_cloud_new / (ssa_total*od_total)
	358	end where
	359	end if
	360	call calc_two_stream_gammas_lw(ng, &
	361	& ssa_total, g_total, gamma1, gamma2)
	362	call calc_reflectance_transmittance_lw(ng, &
	363	& od_total, gamma1, gamma2, &
	364	& planck_hl(:,jlev,jcol), planck_hl(:,jlev+1,jcol), &
	365	& reflectance(:,jreg,jlev), transmittance(:,jreg,jlev), &
[4489]	366	& source_up(:,jreg,jlev), source_dn(:,jreg,jlev))
[3908]	367	else
	368	! No-scattering case: use simpler functions for
	369	! transmission and emission
	370	call calc_no_scattering_transmittance_lw(ng, od_total, &
	371	& planck_hl(:,jlev,jcol), planck_hl(:,jlev+1, jcol), &
[4489]	372	& transmittance(:,jreg,jlev), source_up(:,jreg,jlev), source_dn(:,jreg,jlev))
[3908]	373	reflectance(:,jreg,jlev) = 0.0_jprb
	374	end if
[4489]	375	end do
[3908]	376	! Emission is scaled by the size of each region
[5158]	377	DO jreg = 1,nregions
[4489]	378	source_up(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_up(:,jreg,jlev)
	379	source_dn(:,jreg,jlev) = region_fracs(jreg,jlev,jcol) * source_dn(:,jreg,jlev)
[3908]	380	end do
	381	end if
	382
	383	end do ! Loop over levels
	384
	385	! --------------------------------------------------------
[4489]	386	! Section 5: Compute total sources and albedos at each half level
[3908]	387	! --------------------------------------------------------
	388
	389	total_albedo(:,:,:) = 0.0_jprb
	390	total_source(:,:,:) = 0.0_jprb
	391
	392	! Calculate the upwelling radiation emitted by the surface, and
	393	! copy the surface albedo into total_albedo
[5158]	394	DO jreg = 1,nregions
	395	DO jg = 1,ng
[3908]	396	! region_fracs(jreg,nlev,jcol) is the fraction of each region in the
	397	! lowest model level
	398	total_source(jg,jreg,nlev+1) = region_fracs(jreg,nlev,jcol)*emission(jg,jcol)
	399	total_albedo(jg,jreg,nlev+1) = albedo(jg,jcol)
	400	end do
	401	end do
	402
	403	! Work up from the surface computing the total albedo of the
	404	! atmosphere and the total upwelling due to emission below each
	405	! level below using the adding method
[5158]	406	DO jlev = nlev,i_cloud_top,-1
[3908]	407
	408	total_albedo_below = 0.0_jprb
	409
	410	if (is_clear_sky_layer(jlev)) then
	411	total_albedo_below = 0.0_jprb
	412	total_source_below = 0.0_jprb
[5158]	413	DO jg = 1,ng
[4489]	414	inv_denom(jg,1) = 1.0_jprb &
	415	& / (1.0_jprb - total_albedo(jg,1,jlev+1)*reflectance(jg,1,jlev))
	416	total_albedo_below(jg,1) = reflectance(jg,1,jlev) &
	417	& + transmittance(jg,1,jlev)transmittance(jg,1,jlev)total_albedo(jg,1,jlev+1) &
	418	& * inv_denom(jg,1)
	419	total_source_below(jg,1) = source_up(jg,1,jlev) &
	420	& + transmittance(jg,1,jlev)*(total_source(jg,1,jlev+1) &
	421	& + total_albedo(jg,1,jlev+1)*source_dn(jg,1,jlev)) &
	422	& * inv_denom(jg,1)
	423	end do
[3908]	424	else
	425	inv_denom = 1.0_jprb / (1.0_jprb - total_albedo(:,:,jlev+1)*reflectance(:,:,jlev))
	426	total_albedo_below = reflectance(:,:,jlev) &
	427	& + transmittance(:,:,jlev)transmittance(:,:,jlev)total_albedo(:,:,jlev+1) &
	428	& * inv_denom
[4489]	429	total_source_below = source_up(:,:,jlev) &
[3908]	430	& + transmittance(:,:,jlev)*(total_source(:,:,jlev+1) &
[4489]	431	& + total_albedo(:,:,jlev+1)*source_dn(:,:,jlev)) &
[3908]	432	& * inv_denom
	433	end if
	434
	435	! Account for cloud overlap when converting albedo below a
	436	! layer interface to the equivalent values just above
[5185]	437	if (is_clear_sky_layer(jlev) .AND. is_clear_sky_layer(jlev-1)) then
[3908]	438	total_albedo(:,:,jlev) = total_albedo_below(:,:)
	439	total_source(:,:,jlev) = total_source_below(:,:)
	440	else
	441	total_source(:,:,jlev) = singlemat_x_vec(ng,ng,nregions,&
	442	& u_matrix(:,:,jlev,jcol), total_source_below)
	443	! Use overlap matrix and exclude "anomalous" horizontal
	444	! transport described by Shonk & Hogan (2008). Therefore,
	445	! the operation we perform is essentially diag(total_albedo)
	446	! = matmul(transpose(v_matrix), diag(total_albedo_below)).
[5158]	447	DO jreg = 1,nregions
	448	DO jreg2 = 1,nregions
[3908]	449	total_albedo(:,jreg,jlev) &
	450	& = total_albedo(:,jreg,jlev) &
	451	& + total_albedo_below(:,jreg2) &
	452	& * v_matrix(jreg2,jreg,jlev,jcol)
	453
	454	end do
	455	end do
	456
	457	end if
	458
	459	end do ! Reverse loop over levels
	460
	461	! --------------------------------------------------------
[4489]	462	! Section 6: Copy over downwelling fluxes above cloud top
[3908]	463	! --------------------------------------------------------
[5158]	464	DO jlev = 1,i_cloud_top
[4489]	465	if (config%do_clear) then
	466	! Clear-sky fluxes have already been averaged: use these
	467	flux%lw_dn(jcol,jlev) = flux%lw_dn_clear(jcol,jlev)
	468	if (config%do_save_spectral_flux) then
	469	flux%lw_dn_band(:,jcol,jlev) = flux%lw_dn_clear_band(:,jcol,jlev)
	470	end if
	471	else
	472	flux%lw_dn(jcol,:) = sum(flux_dn_clear(:,jlev))
	473	if (config%do_save_spectral_flux) then
	474	call indexed_sum(flux_dn_clear(:,jlev), &
	475	& config%i_spec_from_reordered_g_lw, &
	476	& flux%lw_dn_band(:,jcol,jlev))
	477	end if
	478	end if
	479	end do
[3908]	480
[4489]	481	! --------------------------------------------------------
	482	! Section 7: Compute fluxes up to top-of-atmosphere
	483	! --------------------------------------------------------
[3908]	484
[4489]	485	! Compute the fluxes just above the highest cloud
	486	flux_up(:,1) = total_source(:,1,i_cloud_top) &
	487	& + total_albedo(:,1,i_cloud_top)*flux_dn_clear(:,i_cloud_top)
	488	flux_up(:,2:) = 0.0_jprb
	489	flux%lw_up(jcol,i_cloud_top) = sum(flux_up(:,1))
[3908]	490	if (config%do_save_spectral_flux) then
[4489]	491	call indexed_sum(flux_up(:,1), &
[3908]	492	& config%i_spec_from_reordered_g_lw, &
[4489]	493	& flux%lw_up_band(:,jcol,i_cloud_top))
	494	end if
[5158]	495	DO jlev = i_cloud_top-1,1,-1
[4489]	496	flux_up(:,1) = trans_clear(:,jlev)*flux_up(:,1) + source_up_clear(:,jlev)
	497	flux%lw_up(jcol,jlev) = sum(flux_up(:,1))
	498	if (config%do_save_spectral_flux) then
	499	call indexed_sum(flux_up(:,1), &
[3908]	500	& config%i_spec_from_reordered_g_lw, &
[4489]	501	& flux%lw_up_band(:,jcol,jlev))
[3908]	502	end if
[4489]	503	end do
[3908]	504
[4489]	505	! --------------------------------------------------------
	506	! Section 8: Compute fluxes down to surface
	507	! --------------------------------------------------------
	508
	509	! Copy over downwelling spectral fluxes at top of first
	510	! scattering layer, using overlap matrix to translate to the
	511	! regions of the first layer of cloud
[5158]	512	DO jreg = 1,nregions
[4489]	513	flux_dn(:,jreg) = v_matrix(jreg,1,i_cloud_top,jcol)*flux_dn_clear(:,i_cloud_top)
	514	end do
	515
[3908]	516	! Final loop back down through the atmosphere to compute fluxes
[5158]	517	DO jlev = i_cloud_top,nlev
[3908]	518
	519	if (is_clear_sky_layer(jlev)) then
[5158]	520	DO jg = 1,ng
[4489]	521	flux_dn(jg,1) = (transmittance(jg,1,jlev)*flux_dn(jg,1) &
	522	& + reflectance(jg,1,jlev)*total_source(jg,1,jlev+1) + source_dn(jg,1,jlev) ) &
	523	& / (1.0_jprb - reflectance(jg,1,jlev)*total_albedo(jg,1,jlev+1))
	524	flux_up(jg,1) = total_source(jg,1,jlev+1) + flux_dn(jg,1)*total_albedo(jg,1,jlev+1)
	525	end do
[3908]	526	flux_dn(:,2:) = 0.0_jprb
	527	flux_up(:,2:) = 0.0_jprb
	528	else
	529	flux_dn = (transmittance(:,:,jlev)*flux_dn &
[4489]	530	& + reflectance(:,:,jlev)*total_source(:,:,jlev+1) + source_dn(:,:,jlev) ) &
[3908]	531	& / (1.0_jprb - reflectance(:,:,jlev)*total_albedo(:,:,jlev+1))
	532	flux_up = total_source(:,:,jlev+1) + flux_dn*total_albedo(:,:,jlev+1)
	533	end if
	534
	535	if (.not. (is_clear_sky_layer(jlev) &
[5185]	536	& .AND. is_clear_sky_layer(jlev+1))) then
[3908]	537	! Account for overlap rules in translating fluxes just above
	538	! a layer interface to the values just below
	539	flux_dn_below = singlemat_x_vec(ng,ng,nregions, &
	540	& v_matrix(:,:,jlev+1,jcol), flux_dn)
	541	flux_dn = flux_dn_below
	542	end if ! Otherwise the fluxes in each region are the same so
	543	! nothing to do
	544
	545	! Store the broadband fluxes
	546	flux%lw_up(jcol,jlev+1) = sum(sum(flux_up,1))
	547	flux%lw_dn(jcol,jlev+1) = sum(sum(flux_dn,1))
	548
	549	! Save the spectral fluxes if required
	550	if (config%do_save_spectral_flux) then
	551	call indexed_sum(sum(flux_up,2), &
	552	& config%i_spec_from_reordered_g_lw, &
	553	& flux%lw_up_band(:,jcol,jlev+1))
	554	call indexed_sum(sum(flux_dn,2), &
	555	& config%i_spec_from_reordered_g_lw, &
	556	& flux%lw_dn_band(:,jcol,jlev+1))
[4489]	557	end if
[3908]	558
	559	end do ! Final loop over levels
	560
	561	! Store surface spectral downwelling fluxes, which at this point
	562	! are at the surface
	563	flux%lw_dn_surf_g(:,jcol) = sum(flux_dn,2)
	564
	565	! Compute the longwave derivatives needed by Hogan and Bozzo
	566	! (2015) approximate radiation update scheme
	567	if (config%do_lw_derivatives) then
	568	! Note that at this point flux_up contains the spectral
	569	! fluxes into the regions of the lowest layer; we sum over
	570	! regions first to provide a simple spectral flux upwelling
	571	! from the surface
	572	call calc_lw_derivatives_region(ng, nlev, nregions, jcol, transmittance, &
	573	& u_matrix(:,:,:,jcol), sum(flux_up,2), flux%lw_derivatives)
	574	end if
	575
	576	end do ! Loop over columns
	577
	578	if (lhook) call dr_hook('radiation_tripleclouds_lw:solver_tripleclouds_lw',1,hook_handle)
	579
	580	end subroutine solver_tripleclouds_lw
	581
	582	end module radiation_tripleclouds_lw

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

Download in other formats:

Original Format