Context Navigation

radiation_ifs_rrtm.F90 @ 4629

Last change on this file since 4629 was 4489, checked in by lguez, 19 months ago
Merge LMDZ_ECRad branch back into trunk!
File size: 35.8 KB

Rev	Line
[3908]	1	! radiation_ifs_rrtm.F90 - Interface to IFS implementation of RRTM-G
	2	!
	3	! (C) Copyright 2015- ECMWF.
	4	!
	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.
	7	!
	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.
	11	!
	12	! Author: Robin Hogan
	13	! Email: r.j.hogan@ecmwf.int
	14	!
	15	! Modifications
	16	! 2017-04-11 R. Hogan Receive "surface" dummy argument
	17	! 2017-09-08 R. Hogan Reverted some changes
	18	! 2017-10-18 R. Hogan Added planck_function public function
	19	! 2018-01-11 R. Hogan Added optional spectral scaling of incoming solar radiation
	20	! 2018-02-22 R. Hogan Optimized reverse indexing of heights
	21	! 2018-05-05 R. Hogan gas_optics can be called for reduced number of levels
	22	! 2019-01-02 R. Hogan Initialize shortwave props to zero in case sun below horizon
	23
	24	module radiation_ifs_rrtm
	25
	26	implicit none
	27
	28	public :: setup_gas_optics, gas_optics, planck_function, set_gas_units
	29
	30	contains
	31
	32	!---------------------------------------------------------------------
	33	! Setup the IFS implementation of RRTM-G gas absorption model
	34	subroutine setup_gas_optics(config, directory)
	35
	36	use yoerrtm, only : jpglw
	37	use yoesrtm, only : jpgsw
	38	use yoerrtftr, only : ngb_lw => ngb
	39	use yoesrtm, only : ngb_sw => ngbsw
	40	use yomhook, only : lhook, dr_hook
	41
	42	use radiation_config
	43
	44	type(config_type), intent(inout), target :: config
	45	character(len=*), intent(in) :: directory
	46
	47	integer :: irep ! For implied do
	48
	49	integer, parameter :: RRTM_GPOINT_REORDERING_LW(140) = (/ &
	50	& 89, 90, 139, 77, 137, 69, 131, 97, 91, 70, 78, 71, 53, 72, 123, 54, 79, 98, &
	51	& 92, 55, 80, 132, 124, 81, 73, 56, 99, 82, 57, 23, 125, 100, 24, 74, 93, 58, 25, &
	52	& 83, 126, 75, 26, 11, 101, 133, 59, 27, 76, 140, 12, 84, 102, 94, 28, 127, 85, &
	53	& 13, 39, 60, 86, 103, 87, 109, 14, 29, 115, 40, 95, 15, 61, 88, 41, 110, 104, 1, &
	54	& 116, 42, 30, 134, 128, 138, 96, 62, 16, 43, 117, 63, 111, 44, 2, 64, 31, 65, &
	55	& 105, 17, 45, 66, 118, 32, 3, 33, 67, 18, 129, 135, 46, 112, 34, 106, 68, 35, 4, &
	56	& 119, 36, 47, 107, 19, 37, 38, 113, 48, 130, 5, 120, 49, 108, 20, 50, 51, 114, &
	57	& 21, 121, 52, 136, 122, 6, 22, 7, 8, 9, 10 &
	58	& /)
	59	integer, parameter :: RRTM_GPOINT_REORDERING_SW(112) = (/ &
	60	& 35, 45, 19, 27, 36, 57, 20, 46, 58, 21, 28, 67, 55, 68, 37, 1, 69, 22, 29, 59, &
	61	& 78, 101, 79, 77, 70, 76, 47, 75, 30, 81, 60, 102, 80, 82, 23, 2, 83, 84, 85, &
	62	& 86, 103, 61, 31, 87, 56, 38, 71, 48, 88, 3, 62, 89, 24, 7, 49, 32, 104, 72, 90, &
	63	& 63, 39, 4, 8, 50, 91, 64, 40, 33, 25, 51, 95, 96, 73, 65, 9, 41, 97, 92, 105, &
	64	& 52, 5, 98, 10, 42, 99, 100, 66, 11, 74, 34, 53, 26, 6, 106, 12, 43, 13, 54, 93, &
	65	& 44, 107, 94, 14, 108, 15, 16, 109, 17, 18, 110, 111, 112 &
	66	& /)
	67
	68	real(jprb) :: hook_handle
	69
[4489]	70	!#include "surdi.intfb.h"
[3908]	71	#include "surrtab.intfb.h"
	72	#include "surrtpk.intfb.h"
	73	#include "surrtrf.intfb.h"
	74	#include "rrtm_init_140gp.intfb.h"
	75	#include "srtm_init.intfb.h"
	76
	77	if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',0,hook_handle)
	78
	79	! The IFS implementation of RRTMG uses many global variables. In
	80	! the IFS these will have been set up already; otherwise set them
	81	! up now.
	82	if (config%do_setup_ifsrrtm) then
[4489]	83	!call SURDI
[3908]	84	call SURRTAB
	85	call SURRTPK
	86	call SURRTRF
	87	call RRTM_INIT_140GP(directory)
	88	call SRTM_INIT(directory)
	89	end if
	90
[4489]	91	! Cloud and aerosol properties can only be defined per band
	92	config%do_cloud_aerosol_per_sw_g_point = .false.
	93	config%do_cloud_aerosol_per_lw_g_point = .false.
	94
[3908]	95	config%n_g_sw = jpgsw
	96	config%n_g_lw = jpglw
	97	config%n_bands_sw = 14
	98	config%n_bands_lw = 16
	99
	100	! Wavenumber ranges of each band may be needed so that the user
	101	! can compute UV and photosynthetically active radiation for a
	102	! particular wavelength range
[4489]	103	call config%gas_optics_sw%spectral_def%allocate_bands_only( &
	104	& [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
	105	& 8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
	106	& [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
	107	& 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
	108	call config%gas_optics_lw%spectral_def%allocate_bands_only( &
	109	& [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
	110	& 1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
	111	& [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
	112	& 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
	113
[3908]	114	allocate(config%i_band_from_g_sw (config%n_g_sw))
	115	allocate(config%i_band_from_g_lw (config%n_g_lw))
	116	allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
	117	allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
	118	allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
	119	allocate(config%i_g_from_reordered_g_lw(config%n_g_lw))
	120
	121	! Shortwave starts at 16: need to start at 1
	122	config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
	123	config%i_band_from_g_lw = ngb_lw
	124
	125	if (config%i_solver_sw == ISolverSpartacus) then
	126	! SPARTACUS requires g points ordered in approximately
	127	! increasing order of optical depth
	128	config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
	129	else
	130	! Implied-do for no reordering
	131	! config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
	132	config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
	133	end if
	134
	135	if (config%i_solver_lw == ISolverSpartacus) then
	136	! SPARTACUS requires g points ordered in approximately
	137	! increasing order of optical depth
	138	config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
	139	else
	140	! Implied-do for no reordering
	141	config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
	142	end if
	143
	144	config%i_band_from_reordered_g_sw &
	145	= config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
	146
	147	config%i_band_from_reordered_g_lw &
	148	= config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
	149
	150	! The i_spec_* variables are used solely for storing spectral
	151	! data, and this can either be by band or by g-point
	152	if (config%do_save_spectral_flux) then
	153	if (config%do_save_gpoint_flux) then
	154	config%n_spec_sw = config%n_g_sw
	155	config%n_spec_lw = config%n_g_lw
	156	config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
	157	config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
	158	else
	159	config%n_spec_sw = config%n_bands_sw
	160	config%n_spec_lw = config%n_bands_lw
	161	config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
	162	config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
	163	end if
	164	else
	165	config%n_spec_sw = 0
	166	config%n_spec_lw = 0
	167	nullify(config%i_spec_from_reordered_g_sw)
	168	nullify(config%i_spec_from_reordered_g_lw)
	169	end if
	170
	171	if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',1,hook_handle)
	172
	173	end subroutine setup_gas_optics
	174
	175
	176	!---------------------------------------------------------------------
	177	! Scale gas mixing ratios according to required units
	178	subroutine set_gas_units(gas)
	179
	180	use radiation_gas, only : gas_type, IMassMixingRatio
	181	type(gas_type), intent(inout) :: gas
	182
	183	call gas%set_units(IMassMixingRatio)
	184
	185	end subroutine set_gas_units
	186
	187
	188	!---------------------------------------------------------------------
	189	! Compute gas optical depths, shortwave scattering, Planck function
	190	! and incoming shortwave radiation at top-of-atmosphere
	191	subroutine gas_optics(ncol,nlev,istartcol,iendcol, &
	192	& config, single_level, thermodynamics, gas, &
	193	& od_lw, od_sw, ssa_sw, lw_albedo, planck_hl, lw_emission, &
	194	& incoming_sw)
	195
	196	use parkind1, only : jprb, jpim
	197
	198	USE PARRRTM , ONLY : JPBAND, JPXSEC, JPINPX
	199	USE YOERRTM , ONLY : JPGPT_LW => JPGPT
	200	USE YOESRTM , ONLY : JPGPT_SW => JPGPT
	201	!USE YOMDIMV , ONLY : YRDIMV
	202	use yomhook , only : lhook, dr_hook
	203
	204	use radiation_config, only : config_type, ISolverSpartacus
	205	use radiation_thermodynamics, only : thermodynamics_type
	206	use radiation_single_level, only : single_level_type
	207	use radiation_gas
	208
	209	integer, intent(in) :: ncol ! number of columns
	210	integer, intent(in) :: nlev ! number of model levels
	211	integer, intent(in) :: istartcol, iendcol ! range of columns to process
	212	type(config_type), intent(in) :: config
	213	type(single_level_type), intent(in) :: single_level
	214	type(thermodynamics_type),intent(in) :: thermodynamics
	215	type(gas_type), intent(in) :: gas
	216
	217	! Longwave albedo of the surface
	218	real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
	219	& intent(in), optional :: lw_albedo
	220
	221	! Gaseous layer optical depth in longwave and shortwave, and
	222	! shortwave single scattering albedo (i.e. fraction of extinction
	223	! due to Rayleigh scattering) at each g-point
	224	real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: &
	225	& od_lw
	226	real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
	227	& od_sw, ssa_sw
	228
	229	! The Planck function (emitted flux from a black body) at half
	230	! levels at each longwave g-point
	231	real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), &
	232	& intent(out), optional :: planck_hl
	233	! Planck function for the surface (W m-2)
	234	real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
	235	& intent(out), optional :: lw_emission
	236
	237	! The incoming shortwave flux into a plane perpendicular to the
	238	! incoming radiation at top-of-atmosphere in each of the shortwave
	239	! g-points
	240	real(jprb), dimension(config%n_g_sw,istartcol:iendcol), &
	241	& intent(out), optional :: incoming_sw
	242
	243	real(jprb) :: incoming_sw_scale(istartcol:iendcol)
	244
	245	! The variables in capitals are used in the same way as the
	246	! equivalent routine in the IFS
	247
	248	real(jprb) :: ZOD_LW(JPGPT_LW,nlev,istartcol:iendcol) ! Note ordering of dimensions
	249	real(jprb) :: ZOD_SW(istartcol:iendcol,nlev,JPGPT_SW)
	250	real(jprb) :: ZSSA_SW(istartcol:iendcol,nlev,JPGPT_SW)
	251	real(jprb) :: ZINCSOL(istartcol:iendcol,JPGPT_SW)
	252
	253	real(jprb) :: ZCOLMOL(istartcol:iendcol,nlev)
	254	real(jprb) :: ZCOLDRY(istartcol:iendcol,nlev)
	255	real(jprb) :: ZWBRODL(istartcol:iendcol,nlev) !BROADENING GASES,column density (mol/cm2)
	256	real(jprb) :: ZCOLBRD(istartcol:iendcol,nlev) !BROADENING GASES, column amount
	257	real(jprb) :: ZWKL(istartcol:iendcol,JPINPX,nlev)
	258
	259	real(jprb) :: ZWX(istartcol:iendcol,JPXSEC,nlev) ! Amount of trace gases
	260
	261	real(jprb) :: ZFLUXFAC, ZPI
	262
	263	! - from AER
	264	real(jprb) :: ZTAUAERL(istartcol:iendcol,nlev,JPBAND)
	265
	266	!- from INTFAC
	267	real(jprb) :: ZFAC00(istartcol:iendcol,nlev)
	268	real(jprb) :: ZFAC01(istartcol:iendcol,nlev)
	269	real(jprb) :: ZFAC10(istartcol:iendcol,nlev)
	270	real(jprb) :: ZFAC11(istartcol:iendcol,nlev)
	271
	272	!- from FOR
	273	real(jprb) :: ZFORFAC(istartcol:iendcol,nlev)
	274	real(jprb) :: ZFORFRAC(istartcol:iendcol,nlev)
	275	integer :: INDFOR(istartcol:iendcol,nlev)
	276
	277	!- from MINOR
	278	integer :: INDMINOR(istartcol:iendcol,nlev)
	279	real(jprb) :: ZSCALEMINOR(istartcol:iendcol,nlev)
	280	real(jprb) :: ZSCALEMINORN2(istartcol:iendcol,nlev)
	281	real(jprb) :: ZMINORFRAC(istartcol:iendcol,nlev)
	282
	283	real(jprb) :: &
	284	& ZRAT_H2OCO2(istartcol:iendcol,nlev),ZRAT_H2OCO2_1(istartcol:iendcol,nlev), &
	285	& ZRAT_H2OO3(istartcol:iendcol,nlev) ,ZRAT_H2OO3_1(istartcol:iendcol,nlev), &
	286	& ZRAT_H2ON2O(istartcol:iendcol,nlev),ZRAT_H2ON2O_1(istartcol:iendcol,nlev), &
	287	& ZRAT_H2OCH4(istartcol:iendcol,nlev),ZRAT_H2OCH4_1(istartcol:iendcol,nlev), &
	288	& ZRAT_N2OCO2(istartcol:iendcol,nlev),ZRAT_N2OCO2_1(istartcol:iendcol,nlev), &
	289	& ZRAT_O3CO2(istartcol:iendcol,nlev) ,ZRAT_O3CO2_1(istartcol:iendcol,nlev)
	290
	291	!- from INTIND
	292	integer :: JP(istartcol:iendcol,nlev)
	293	integer :: JT(istartcol:iendcol,nlev)
	294	integer :: JT1(istartcol:iendcol,nlev)
	295
	296	!- from PRECISE
	297	real(jprb) :: ZONEMINUS, ZONEMINUS_ARRAY(istartcol:iendcol)
	298
	299	!- from PROFDATA
	300	real(jprb) :: ZCOLH2O(istartcol:iendcol,nlev)
	301	real(jprb) :: ZCOLCO2(istartcol:iendcol,nlev)
	302	real(jprb) :: ZCOLO3(istartcol:iendcol,nlev)
	303	real(jprb) :: ZCOLN2O(istartcol:iendcol,nlev)
	304	real(jprb) :: ZCOLCH4(istartcol:iendcol,nlev)
	305	real(jprb) :: ZCOLO2(istartcol:iendcol,nlev)
	306	real(jprb) :: ZCO2MULT(istartcol:iendcol,nlev)
	307	integer :: ILAYTROP(istartcol:iendcol)
	308	integer :: ILAYSWTCH(istartcol:iendcol)
	309	integer :: ILAYLOW(istartcol:iendcol)
	310
	311	!- from PROFILE
	312	real(jprb) :: ZPAVEL(istartcol:iendcol,nlev)
	313	real(jprb) :: ZTAVEL(istartcol:iendcol,nlev)
	314	real(jprb) :: ZPZ(istartcol:iendcol,0:nlev)
	315	real(jprb) :: ZTZ(istartcol:iendcol,0:nlev)
	316
	317	!- from SELF
	318	real(jprb) :: ZSELFFAC(istartcol:iendcol,nlev)
	319	real(jprb) :: ZSELFFRAC(istartcol:iendcol,nlev)
	320	integer :: INDSELF(istartcol:iendcol,nlev)
	321
	322	!- from SP
	323	real(jprb) :: ZPFRAC(istartcol:iendcol,JPGPT_LW,nlev)
	324
	325	!- from SURFACE
	326	integer :: IREFLECT(istartcol:iendcol)
	327
	328	real(jprb) :: pressure_fl(ncol, nlev), temperature_fl(ncol, nlev)
	329
	330	! If nlev is less than the number of heights at which gas mixing
	331	! ratios are stored, then we assume that the lower part of the
	332	! atmosphere is required. This enables nlev=1 to be passed in to
	333	! the routine, in which case the gas properties of the lowest
	334	! layer are provided, useful for canopy radiative transfer.
	335	integer :: istartlev, iendlev
	336
	337	integer :: jlev, jgreorder, jg, ig, iband, jcol
	338
	339	real(jprb) :: hook_handle
	340
	341	#include "rrtm_prepare_gases.intfb.h"
	342	#include "rrtm_setcoef_140gp.intfb.h"
	343	#include "rrtm_gas_optical_depth.intfb.h"
	344	#include "srtm_setcoef.intfb.h"
	345	#include "srtm_gas_optical_depth.intfb.h"
	346
	347	if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',0,hook_handle)
	348
	349	! Compute start and end levels for indexing the gas mixing ratio
	350	! and thermodynamics arrays
	351	iendlev = ubound(gas%mixing_ratio,2)
	352	istartlev = iendlev - nlev + 1
	353
	354	ZPI = 2.0_jprb*ASIN(1.0_jprb)
	355	ZFLUXFAC = ZPI * 1.E+4
	356	ZONEMINUS = 1.0_jprb - 1.0e-6_jprb
	357	ZONEMINUS_ARRAY = ZONEMINUS
	358
	359	! if (.not. associated(YRDIMV)) then
	360	! allocate(YRDIMV)
	361	! YRDIMV%NFLEVG = nlev
	362	! end if
	363
[4489]	364	do jlev=1,nlev
	365	do jcol= istartcol,iendcol
	366	pressure_fl(jcol,jlev) &
	367	& = 0.5_jprb * (thermodynamics%pressure_hl(jcol,jlev+istartlev-1) &
	368	& +thermodynamics%pressure_hl(jcol,jlev+istartlev))
	369	temperature_fl(jcol,jlev) &
	370	& = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev+istartlev-1) &
	371	& +thermodynamics%temperature_hl(jcol,jlev+istartlev))
	372	end do
	373	end do
[3908]	374
	375	! Check we have gas mixing ratios in the right units
	376	call gas%assert_units(IMassMixingRatio)
	377
	378	! Warning: O2 is hard-coded within the following function so the
	379	! user-provided concentrations of this gas are ignored for both
	380	! the longwave and shortwave
	381	CALL RRTM_PREPARE_GASES &
	382	& ( istartcol, iendcol, ncol, nlev, &
	383	& thermodynamics%pressure_hl(:,istartlev:iendlev+1), &
	384	& pressure_fl, &
	385	& thermodynamics%temperature_hl(:,istartlev:iendlev+1), &
	386	& temperature_fl, &
	387	& gas%mixing_ratio(:,istartlev:iendlev,IH2O), &
	388	& gas%mixing_ratio(:,istartlev:iendlev,ICO2), &
	389	& gas%mixing_ratio(:,istartlev:iendlev,ICH4), &
	390	& gas%mixing_ratio(:,istartlev:iendlev,IN2O), &
	391	& gas%mixing_ratio(:,istartlev:iendlev,INO2), &
	392	& gas%mixing_ratio(:,istartlev:iendlev,ICFC11), &
	393	& gas%mixing_ratio(:,istartlev:iendlev,ICFC12), &
	394	& gas%mixing_ratio(:,istartlev:iendlev,IHCFC22), &
	395	& gas%mixing_ratio(:,istartlev:iendlev,ICCl4), &
	396	& gas%mixing_ratio(:,istartlev:iendlev,IO3), &
	397	& ZCOLDRY, ZWBRODL,ZWKL, ZWX, &
	398	& ZPAVEL , ZTAVEL , ZPZ , ZTZ, IREFLECT)
	399
	400	CALL RRTM_SETCOEF_140GP &
	401	&( istartcol, iendcol, nlev , ZCOLDRY , ZWBRODL , ZWKL , &
	402	& ZFAC00 , ZFAC01 , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
	403	& ZCOLH2O, ZCOLCO2 , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
	404	& ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
	405	& INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
	406	& ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
	407	& ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
	408	& ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)
	409
[4489]	410	ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
[3908]	411
	412	CALL RRTM_GAS_OPTICAL_DEPTH &
	413	&( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
	414	& ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
	415	& JP, JT, JT1, ZONEMINUS ,&
	416	& ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
	417	& ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
	418	& INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
	419	& ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
	420	& ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
	421	& ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)
	422
	423	if (present(lw_albedo)) then
	424
	425	call planck_function_atmos(nlev, istartcol, iendcol, config, &
	426	& thermodynamics, ZPFRAC, planck_hl)
	427
	428	if (single_level%is_simple_surface) then
	429	call planck_function_surf(istartcol, iendcol, config, &
	430	& single_level%skin_temperature, ZPFRAC(:,:,1), &
	431	& lw_emission)
	432
	433	! The following can be used to extract the parameters defined at
	434	! the top of the planck_function routine below:
	435	!write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
	436	!& sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
	437
	438	! lw_emission at this point is actually the planck function of
	439	! the surface
	440	lw_emission = lw_emission * (1.0_jprb - lw_albedo)
	441	else
[4489]	442	! Longwave emission has already been computed
[3908]	443	if (config%use_canopy_full_spectrum_lw) then
	444	lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
	445	else
	446	lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
	447	& config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
	448	end if
	449	end if
	450
	451	end if
	452
	453	if (config%i_solver_lw == ISolverSpartacus) then
	454	! if (.true.) then
	455	! We need to rearrange the gas optics info in memory: reordering
	456	! the g points in order of approximately increasing optical
	457	! depth (for efficient 3D processing on only the regions of the
	458	! spectrum that are optically thin for gases) and reorder in
	459	! pressure since the the functions above treat pressure
	460	! decreasing with increasing index. Note that the output gas
	461	! arrays have dimensions in a different order to the inputs,
	462	! so there is some inefficiency here.
	463	do jgreorder = 1,config%n_g_lw
	464	iband = config%i_band_from_reordered_g_lw(jgreorder)
	465	ig = config%i_g_from_reordered_g_lw(jgreorder)
	466
	467	! Top-of-atmosphere half level
	468	do jlev = 1,nlev
	469	do jcol = istartcol,iendcol
	470	! Some g points can return negative optical depths;
	471	! specifically original g points 54-56 which causes
	472	! unphysical single-scattering albedo when combined with
	473	! aerosol
	474	od_lw(jgreorder,jlev,jcol) &
	475	& = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
	476	end do
	477	end do
	478	end do
	479	else
	480	! G points have not been reordered
	481	do jcol = istartcol,iendcol
	482	do jlev = 1,nlev
	483	! Check for negative optical depth
	484	od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
	485	end do
	486	end do
	487	end if
	488
	489	CALL SRTM_SETCOEF &
	490	& ( istartcol, iendcol, nlev,&
	491	& ZPAVEL , ZTAVEL,&
	492	& ZCOLDRY , ZWKL,&
	493	& ILAYTROP,&
	494	& ZCOLCH4 , ZCOLCO2 , ZCOLH2O , ZCOLMOL , ZCOLO2 , ZCOLO3,&
	495	& ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF, &
	496	& ZFAC00 , ZFAC01 , ZFAC10 , ZFAC11,&
	497	& JP , JT , JT1 , single_level%cos_sza(istartcol:iendcol) &
	498	& )
	499
	500	! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
	501	! is below the horizon, so we do it here
	502	ZOD_SW(istartcol:iendcol,:,:) = 0.0_jprb
	503	ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
	504	ZINCSOL(istartcol:iendcol,:) = 0.0_jprb
	505
	506	CALL SRTM_GAS_OPTICAL_DEPTH &
	507	&( istartcol, iendcol , nlev , ZONEMINUS_ARRAY,&
	508	& single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
	509	& ZCOLCH4 , ZCOLCO2 , ZCOLH2O, ZCOLMOL , ZCOLO2 , ZCOLO3,&
	510	& ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
	511	& ZFAC00 , ZFAC01 , ZFAC10 , ZFAC11 ,&
	512	& JP , JT , JT1 ,&
	513	& ZOD_SW , ZSSA_SW , ZINCSOL )
	514
	515	! Scale the incoming solar per band, if requested
	516	if (config%use_spectral_solar_scaling) then
[4489]	517	do jg = 1,JPGPT_SW
	518	do jcol = istartcol,iendcol
	519	ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
	520	& single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
	521	end do
	522	end do
[3908]	523	end if
	524
	525	! Scaling factor to ensure that the total solar irradiance is as
	526	! requested. Note that if the sun is below the horizon then
	527	! ZINCSOL will be zero.
	528	if (present(incoming_sw)) then
	529	do jcol = istartcol,iendcol
	530	if (single_level%cos_sza(jcol) > 0.0_jprb) then
[4489]	531	! Added for DWD (2020)
	532	!NEC$ nounroll
[3908]	533	incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
[4489]	534	else
	535	incoming_sw_scale(jcol) = 1.0_jprb
[3908]	536	end if
	537	end do
	538	end if
	539
	540	if (config%i_solver_sw == ISolverSpartacus) then
	541	! if (.true.) then
	542	! Account for reordered g points
	543	do jgreorder = 1,config%n_g_sw
	544	ig = config%i_g_from_reordered_g_sw(jgreorder)
	545	do jlev = 1,nlev
	546	do jcol = istartcol,iendcol
	547	! Check for negative optical depth
	548	od_sw (jgreorder,nlev+1-jlev,jcol) &
	549	& = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
	550	ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
	551	end do
	552	end do
	553	if (present(incoming_sw)) then
	554	incoming_sw(jgreorder,:) &
	555	& = incoming_sw_scale(:) * ZINCSOL(:,ig)
	556	end if
	557	end do
	558	else
	559	! G points have not been reordered
[4489]	560	do jcol = istartcol,iendcol
[3908]	561	do jlev = 1,nlev
[4489]	562	do jg = 1,config%n_g_sw
[3908]	563	! Check for negative optical depth
	564	od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
	565	ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
	566	end do
	567	end do
	568	if (present(incoming_sw)) then
[4489]	569	do jg = 1,config%n_g_sw
	570	incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
	571	end do
[3908]	572	end if
	573	end do
	574	end if
	575
	576	if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',1,hook_handle)
	577
	578	end subroutine gas_optics
	579
	580
	581	!---------------------------------------------------------------------
	582	! Compute Planck function of the atmosphere
	583	subroutine planck_function_atmos(nlev,istartcol,iendcol, &
	584	config, thermodynamics, PFRAC, &
	585	planck_hl)
	586
	587	use parkind1, only : jprb, jpim
	588
	589	USE YOERRTM , ONLY : JPGPT_LW => JPGPT
	590	use yoerrtwn, only : totplnk, delwave
	591
	592	use yomhook, only : lhook, dr_hook
	593
	594	use radiation_config, only : config_type, ISolverSpartacus
	595	use radiation_thermodynamics, only : thermodynamics_type
	596	!use radiation_gas
	597
	598	integer, intent(in) :: nlev ! number of model levels
	599	integer, intent(in) :: istartcol, iendcol ! range of columns to process
	600	type(config_type), intent(in) :: config
	601	type(thermodynamics_type),intent(in) :: thermodynamics
	602	real(jprb), intent(in) :: PFRAC(istartcol:iendcol,JPGPT_LW,nlev)
	603
	604	! The Planck function (emitted flux from a black body) at half
	605	! levels at each longwave g-point
	606	real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), intent(out) :: &
	607	& planck_hl
	608
	609	! Planck function values per band
	610	real(jprb), dimension(istartcol:iendcol, config%n_bands_lw) :: planck_store
	611
	612	! Look-up table variables for Planck function
	613	real(jprb), dimension(istartcol:iendcol) :: frac
	614	integer, dimension(istartcol:iendcol) :: ind
	615
	616	! Temperature (K) of a half-level
	617	real(jprb) :: temperature
	618
[4489]	619	real(jprb) :: factor, planck_tmp(istartcol:iendcol,config%n_g_lw)
[3908]	620	real(jprb) :: ZFLUXFAC
	621
	622	integer :: jlev, jgreorder, jg, ig, iband, jband, jcol, ilevoffset
	623
	624	real(jprb) :: hook_handle
	625
	626	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_atmos',0,hook_handle)
	627
	628	ZFLUXFAC = 2.0_jprbASIN(1.0_jprb) 1.0e4_jprb
	629
	630	! nlev may be less than the number of original levels, in which
	631	! case we assume that the user wants the lower part of the
	632	! atmosphere
	633	ilevoffset = ubound(thermodynamics%temperature_hl,2)-nlev-1
	634
	635	! Work out interpolations: for each half level, the index of the
	636	! lowest interpolation bound, and the fraction into interpolation
	637	! interval
	638	do jlev = 1,nlev+1
	639	do jcol = istartcol,iendcol
	640	temperature = thermodynamics%temperature_hl(jcol,jlev+ilevoffset)
	641	if (temperature < 339.0_jprb .and. temperature >= 160.0_jprb) then
	642	! Linear interpolation between -113 and 66 degC
	643	ind(jcol) = int(temperature - 159.0_jprb)
	644	frac(jcol) = temperature - int(temperature)
	645	else if(temperature >= 339.0_jprb) then
	646	! Extrapolation above 66 degC
	647	ind(jcol) = 180
	648	frac(jcol) = temperature - 339.0_jprb
	649	else
	650	! Cap below -113 degC (to avoid possible negative Planck
	651	! function values)
	652	ind(jcol) = 1
	653	frac(jcol) = 0.0_jprb
	654	end if
	655	end do
	656
	657	! Calculate Planck functions per band
	658	do jband = 1,config%n_bands_lw
	659	factor = zfluxfac * delwave(jband)
	660	do jcol = istartcol,iendcol
	661	planck_store(jcol,jband) = factor &
	662	& * (totplnk(ind(jcol),jband) &
	663	& + frac(jcol)*(totplnk(ind(jcol)+1,jband)-totplnk(ind(jcol),jband)))
	664	end do
	665	end do
	666
	667	if (config%i_solver_lw == ISolverSpartacus) then
	668	! We need to rearrange the gas optics info in memory:
	669	! reordering the g points in order of approximately increasing
	670	! optical depth (for efficient 3D processing on only the
	671	! regions of the spectrum that are optically thin for gases)
	672	! and reorder in pressure since the the functions above treat
	673	! pressure decreasing with increasing index.
	674	if (jlev == 1) then
	675	! Top-of-atmosphere half level - note that PFRAC is on model
	676	! levels not half levels
	677	do jgreorder = 1,config%n_g_lw
	678	iband = config%i_band_from_reordered_g_lw(jgreorder)
	679	ig = config%i_g_from_reordered_g_lw(jgreorder)
	680	planck_hl(jgreorder,1,:) = planck_store(:,iband) &
	681	& * PFRAC(:,ig,nlev)
	682	end do
	683	else
	684	do jgreorder = 1,config%n_g_lw
	685	iband = config%i_band_from_reordered_g_lw(jgreorder)
	686	ig = config%i_g_from_reordered_g_lw(jgreorder)
	687	planck_hl(jgreorder,jlev,:) &
	688	& = planck_store(:,iband) &
	689	& * PFRAC(:,ig,nlev+2-jlev)
	690	end do
	691	end if
	692	else
	693	! G points have not been reordered
	694	if (jlev == 1) then
	695	! Top-of-atmosphere half level - note that PFRAC is on model
	696	! levels not half levels
	697	do jg = 1,config%n_g_lw
	698	iband = config%i_band_from_g_lw(jg)
	699	planck_hl(jg,1,:) = planck_store(:,iband) * PFRAC(:,jg,nlev)
	700	end do
	701	else
	702	do jg = 1,config%n_g_lw
	703	iband = config%i_band_from_g_lw(jg)
[4489]	704	planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
[3908]	705	end do
[4489]	706	do jcol = istartcol,iendcol
	707	planck_hl(:,jlev,jcol) = planck_tmp(jcol,:)
	708	end do
[3908]	709	end if
	710	end if
	711	end do
	712
	713	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_atmos',1,hook_handle)
	714
	715	end subroutine planck_function_atmos
	716
	717
	718	!---------------------------------------------------------------------
	719	! Compute Planck function of the surface
	720	subroutine planck_function_surf(istartcol, iendcol, config, temperature, PFRAC, &
	721	& planck_surf)
	722
	723	use parkind1, only : jprb, jpim
	724
	725	USE YOERRTM , ONLY : JPGPT_LW => JPGPT
	726	use yoerrtwn, only : totplnk, delwave
	727
	728	use yomhook, only : lhook, dr_hook
	729
	730	use radiation_config, only : config_type, ISolverSpartacus
	731	! use radiation_gas
	732
	733	integer, intent(in) :: istartcol, iendcol ! range of columns to process
	734	type(config_type), intent(in) :: config
	735	real(jprb), intent(in) :: temperature(:)
	736
	737	real(jprb), intent(in) :: PFRAC(istartcol:iendcol,JPGPT_LW)
	738
	739	! Planck function of the surface (W m-2)
	740	real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
	741	& intent(out) :: planck_surf
	742
	743	! Planck function values per band
	744	real(jprb), dimension(istartcol:iendcol, config%n_bands_lw) :: planck_store
	745
	746	! Look-up table variables for Planck function
	747	real(jprb), dimension(istartcol:iendcol) :: frac
	748	integer, dimension(istartcol:iendcol) :: ind
	749
	750	! Temperature (K)
	751	real(jprb) :: Tsurf
	752
	753	real(jprb) :: factor
	754	real(jprb) :: ZFLUXFAC
	755
	756	integer :: jgreorder, jg, ig, iband, jband, jcol
	757
	758	real(jprb) :: hook_handle
	759
	760	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_surf',0,hook_handle)
	761
	762	ZFLUXFAC = 2.0_jprbASIN(1.0_jprb) 1.0e4_jprb
	763
	764	! Work out surface interpolations
	765	do jcol = istartcol,iendcol
	766	Tsurf = temperature(jcol)
	767	if (Tsurf < 339.0_jprb .and. Tsurf >= 160.0_jprb) then
	768	! Linear interpolation between -113 and 66 degC
	769	ind(jcol) = int(Tsurf - 159.0_jprb)
	770	frac(jcol) = Tsurf - int(Tsurf)
	771	else if(Tsurf >= 339.0_jprb) then
	772	! Extrapolation above 66 degC
	773	ind(jcol) = 180
	774	frac(jcol) = Tsurf - 339.0_jprb
	775	else
	776	! Cap below -113 degC (to avoid possible negative Planck
	777	! function values)
	778	ind(jcol) = 1
	779	frac(jcol) = 0.0_jprb
	780	end if
	781	end do
	782
	783	! Calculate Planck functions per band
	784	do jband = 1,config%n_bands_lw
	785	factor = zfluxfac * delwave(jband)
	786	do jcol = istartcol,iendcol
	787	planck_store(jcol,jband) = factor &
	788	& * (totplnk(ind(jcol),jband) &
	789	& + frac(jcol)*(totplnk(ind(jcol)+1,jband)-totplnk(ind(jcol),jband)))
	790	end do
	791	end do
	792
	793	if (config%i_solver_lw == ISolverSpartacus) then
	794	! We need to rearrange the gas optics info in memory: reordering
	795	! the g points in order of approximately increasing optical
	796	! depth (for efficient 3D processing on only the regions of
	797	! the spectrum that are optically thin for gases) and reorder
	798	! in pressure since the the functions above treat pressure
	799	! decreasing with increasing index.
	800	do jgreorder = 1,config%n_g_lw
	801	iband = config%i_band_from_reordered_g_lw(jgreorder)
	802	ig = config%i_g_from_reordered_g_lw(jgreorder)
	803	planck_surf(jgreorder,:) = planck_store(:,iband) * PFRAC(:,ig)
	804	end do
	805	else
	806	! G points have not been reordered
	807	do jg = 1,config%n_g_lw
	808	iband = config%i_band_from_g_lw(jg)
	809	planck_surf(jg,:) = planck_store(:,iband) * PFRAC(:,jg)
	810	end do
	811	end if
	812
	813	if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_surf',1,hook_handle)
	814
	815	end subroutine planck_function_surf
	816
	817
	818	!---------------------------------------------------------------------
	819	! Externally facing function for computing the Planck function
	820	! without reference to any gas profile; typically this would be used
	821	! for computing the emission by facets of a complex surface. Note
	822	! that this uses fixed "PFRAC" values, obtained by averaging over
	823	! those derived from RRTM-G for near-surface conditions over a line
	824	! of meridian from the ECMWF model.
	825	subroutine planck_function(config, temperature, planck_surf)
	826
	827	use parkind1, only : jprb, jpim
	828
	829	use radiation_config, only : config_type
	830
	831	type(config_type), intent(in) :: config
	832	real(jprb), intent(in) :: temperature
	833
	834	! Planck function of the surface (W m-2)
	835	real(jprb), dimension(config%n_g_lw), &
	836	& intent(out) :: planck_surf
	837
	838	! Fraction of each band contributed by each g-point within
	839	! it. Since there are 16 bands, this array sums to 16
	840	real(jprb), parameter, dimension(1,140) :: frac &
	841	= reshape( (/ 0.21227E+00, 0.18897E+00, 0.25491E+00, 0.17864E+00, 0.11735E+00, 0.38298E-01, 0.57871E-02, &
	842	& 0.31753E-02, 0.53169E-03, 0.76476E-04, 0.16388E+00, 0.15241E+00, 0.14290E+00, 0.12864E+00, &
	843	& 0.11615E+00, 0.10047E+00, 0.80013E-01, 0.60445E-01, 0.44918E-01, 0.63395E-02, 0.32942E-02, &
	844	& 0.54541E-03, 0.15380E+00, 0.15194E+00, 0.14339E+00, 0.13138E+00, 0.11701E+00, 0.10081E+00, &
	845	& 0.82296E-01, 0.61735E-01, 0.41918E-01, 0.45918E-02, 0.37743E-02, 0.30121E-02, 0.22500E-02, &
	846	& 0.14490E-02, 0.55410E-03, 0.78364E-04, 0.15938E+00, 0.15146E+00, 0.14213E+00, 0.13079E+00, &
	847	& 0.11672E+00, 0.10053E+00, 0.81566E-01, 0.61126E-01, 0.41150E-01, 0.44488E-02, 0.36950E-02, &
	848	& 0.29101E-02, 0.21357E-02, 0.19609E-02, 0.14134E+00, 0.14390E+00, 0.13913E+00, 0.13246E+00, &
	849	& 0.12185E+00, 0.10596E+00, 0.87518E-01, 0.66164E-01, 0.44862E-01, 0.49402E-02, 0.40857E-02, &
	850	& 0.32288E-02, 0.23613E-02, 0.15406E-02, 0.58258E-03, 0.82171E-04, 0.29127E+00, 0.28252E+00, &
	851	& 0.22590E+00, 0.14314E+00, 0.45494E-01, 0.71792E-02, 0.38483E-02, 0.65712E-03, 0.29810E+00, &
	852	& 0.27559E+00, 0.11997E+00, 0.10351E+00, 0.84515E-01, 0.62253E-01, 0.41050E-01, 0.44217E-02, &
	853	& 0.36946E-02, 0.29113E-02, 0.34290E-02, 0.55993E-03, 0.31441E+00, 0.27586E+00, 0.21297E+00, &
	854	& 0.14064E+00, 0.45588E-01, 0.65665E-02, 0.34232E-02, 0.53199E-03, 0.19811E+00, 0.16833E+00, &
	855	& 0.13536E+00, 0.11549E+00, 0.10649E+00, 0.93264E-01, 0.75720E-01, 0.56405E-01, 0.41865E-01, &
	856	& 0.59331E-02, 0.26510E-02, 0.40040E-03, 0.32328E+00, 0.26636E+00, 0.21397E+00, 0.14038E+00, &
	857	& 0.52142E-01, 0.38852E-02, 0.14601E+00, 0.13824E+00, 0.27703E+00, 0.22388E+00, 0.15446E+00, &
	858	& 0.48687E-01, 0.98054E-02, 0.18870E-02, 0.11961E+00, 0.12106E+00, 0.13215E+00, 0.13516E+00, &
	859	& 0.25249E+00, 0.16542E+00, 0.68157E-01, 0.59725E-02, 0.49258E+00, 0.33651E+00, 0.16182E+00, &
	860	& 0.90984E-02, 0.95202E+00, 0.47978E-01, 0.91716E+00, 0.82857E-01, 0.77464E+00, 0.22536E+00 /), (/ 1,140 /) )
	861
	862	call planck_function_surf(1, 1, config, spread(temperature,1,1), &
	863	& frac, planck_surf)
	864
	865	end subroutine planck_function
	866
	867	end module radiation_ifs_rrtm

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

Download in other formats:

Original Format