Context Navigation

radiation_ecckd.F90 @ 5435

Last change on this file since 5435 was 4677, checked in by idelkadi, 16 months ago

Implementation in the LMDZ code of the double call of the ECRAD radiative transfer code to estimate the 3D radiative effect of clouds.

This double call of Ecrad is controlled by the ok_3Deffect logic key.
If this key is enabled, 2 files of parameter configuration "namelists" for ECRAD are required at runtime: namelist_ecrad and namelist_ecrad_s2.
If this key is deactivated, the configuration and initialization part (reading namelist and netcdf files) is performed only once during simulation (1st call to ECRAD). Otherwise, configuration and initialization are performed each time Ecrad is called.

File size: 20.4 KB

Line
1	! radiation_ecckd.F90 - ecCKD generalized gas optics model
2	!
3	! (C) Copyright 2020- 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	! License: see the COPYING file for details
15	!
16
17	module radiation_ecckd
18
19	use parkind1, only : jprb
20	use radiation_gas_constants
21	use radiation_ecckd_gas
22	use radiation_spectral_definition, only : spectral_definition_type
23
24	implicit none
25
26	public
27
28	!---------------------------------------------------------------------
29	! This derived type contains all the data needed to describe a
30	! correlated k-distribution gas optics model created using the ecCKD
31	! tool
32	type ckd_model_type
33
34	! Gas information
35
36	! Number of gases
37	integer :: ngas = 0
38	! Array of individual-gas data objects
39	type(ckd_gas_type), allocatable :: single_gas(:)
40	! Mapping from the "gas codes" in the radiation_gas_constants
41	! module to an index to the single_gas array, where zero means gas
42	! not present (or part of a composite gas)
43	integer :: i_gas_mapping(0:NMaxGases)
44
45	! Coordinates of main look-up table for absorption coeffts
46
47	! Number of pressure and temperature points
48	integer :: npress = 0
49	integer :: ntemp = 0
50	! Natural logarithm of first (lowest) pressure (Pa) and increment
51	real(jprb) :: log_pressure1, d_log_pressure
52	! First temperature profile (K), dimensioned (npress)
53	real(jprb), allocatable :: temperature1(:)
54	! Temperature increment (K)
55	real(jprb) :: d_temperature
56
57	! Look-up table for Planck function
58
59	! Number of entries
60	integer :: nplanck = 0
61	! Temperature of first element of look-up table and increment (K)
62	real(jprb), allocatable :: temperature1_planck
63	real(jprb), allocatable :: d_temperature_planck
64	! Planck function (black body flux into a horizontal plane) in W
65	! m-2, dimensioned (ng,nplanck)
66	real(jprb), allocatable :: planck_function(:,:)
67
68	! Normalized solar irradiance in each g point dimensioned (ng)
69	real(jprb), allocatable :: norm_solar_irradiance(:)
70
71	! Rayleigh molar scattering coefficient in m2 mol-1 in each g
72	! point
73	real(jprb), allocatable :: rayleigh_molar_scat(:)
74
75	! ! Spectral mapping of g points
76
77	! ! Number of wavenumber intervals
78	! integer :: nwav = 0
79
80	! Number of k terms / g points
81	integer :: ng = 0
82
83	! Spectral definition describing bands and g points
84	type(spectral_definition_type) :: spectral_def
85
86	! Shortwave: true, longwave: false
87	logical :: is_sw
88
89	contains
90
91	procedure :: read => read_ckd_model
92	procedure :: calc_optical_depth => calc_optical_depth_ckd_model
93	procedure :: print => print_ckd_model
94	procedure :: calc_planck_function
95	procedure :: calc_incoming_sw
96	! procedure :: deallocate => deallocate_ckd_model
97
98	end type ckd_model_type
99
100
101	contains
102
103	!---------------------------------------------------------------------
104	! Read a complete ecCKD gas optics model from a NetCDF file
105	! "filename"
106	subroutine read_ckd_model(this, filename, iverbose)
107
108	use easy_netcdf, only : netcdf_file
109	!use radiation_io, only : nulerr, radiation_abort
110	use yomhook, only : lhook, dr_hook
111
112	class(ckd_model_type), intent(inout) :: this
113	character(len=*), intent(in) :: filename
114	integer, optional, intent(in) :: iverbose
115
116	type(netcdf_file) :: file
117
118	real(jprb), allocatable :: pressure_lut(:)
119	real(jprb), allocatable :: temperature_full(:,:)
120	real(jprb), allocatable :: temperature_planck(:)
121
122	character(len=512) :: constituent_id
123
124	integer :: iverbose_local
125
126	! Loop counters
127	integer :: jgas, jjgas
128
129	integer :: istart, inext, nchar, i_gas_code
130
131	real(jprb) :: hook_handle
132
133	if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',0,hook_handle)
134
135	if (present(iverbose)) then
136	iverbose_local = iverbose
137	else
138	iverbose_local = 3
139	end if
140
141	call file%open(trim(filename), iverbose=iverbose_local)
142
143	! Read temperature and pressure coordinate variables
144	call file%get('pressure', pressure_lut)
145	this%log_pressure1 = log(pressure_lut(1))
146	this%npress = size(pressure_lut)
147	this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1
148	call file%get('temperature', temperature_full)
149	if (allocated(this%temperature1)) deallocate(this%temperature1)
150	allocate(this%temperature1(this%npress));
151	this%temperature1 = temperature_full(:,1)
152	this%d_temperature = temperature_full(1,2)-temperature_full(1,1)
153	this%ntemp = size(temperature_full,2)
154	deallocate(temperature_full)
155
156	! Read Planck function, or solar irradiance and Rayleigh
157	! scattering coefficient
158	if (file%exists('solar_irradiance')) then
159	this%is_sw = .true.
160	call file%get('solar_irradiance', this%norm_solar_irradiance)
161	this%norm_solar_irradiance = this%norm_solar_irradiance &
162	& / sum(this%norm_solar_irradiance)
163	call file%get('rayleigh_molar_scattering_coeff', &
164	& this%rayleigh_molar_scat)
165	else
166	this%is_sw = .false.
167	call file%get('temperature_planck', temperature_planck)
168	this%nplanck = size(temperature_planck)
169	this%temperature1_planck = temperature_planck(1)
170	this%d_temperature_planck = temperature_planck(2) - temperature_planck(1)
171	deallocate(temperature_planck)
172	call file%get('planck_function', this%planck_function)
173	end if
174
175	! Read the spectral definition information into a separate
176	! derived type
177	call this%spectral_def%read(file);
178	this%ng = this%spectral_def%ng
179
180	! Read gases
181	call file%get('n_gases', this%ngas)
182	if (allocated(this%single_gas)) deallocate(this%single_gas)
183	allocate(this%single_gas(this%ngas))
184	call file%get_global_attribute('constituent_id', constituent_id)
185	nchar = len(trim(constituent_id))
186	istart = 1
187	this%i_gas_mapping = 0
188	do jgas = 1, this%ngas
189	if (jgas < this%ngas) then
190	inext = istart + scan(constituent_id(istart:nchar), ' ')
191	else
192	inext = nchar+2
193	end if
194	! Find gas code
195	i_gas_code = 0
196	do jjgas = 1, NMaxGases
197	if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then
198	i_gas_code = jjgas
199	exit
200	end if
201	end do
202	! if (i_gas_code == 0) then
203	! write(nulerr,'(a,a,a)') '*** Error: Gas "', constituent_id(istart:inext-2), &
204	! & '" not understood'
205	! call radiation_abort('Radiation configuration error')
206	! end if
207	this%i_gas_mapping(i_gas_code) = jgas;
208	call this%single_gas(jgas)%read(file, constituent_id(istart:inext-2), i_gas_code)
209	istart = inext
210	end do
211
212	if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',1,hook_handle)
213
214	end subroutine read_ckd_model
215
216	!---------------------------------------------------------------------
217	! Print a description of the correlated k-distribution model to the
218	! "nulout" unit
219	subroutine print_ckd_model(this)
220
221	use radiation_io, only : nulout
222	use radiation_gas_constants
223
224	class(ckd_model_type), intent(in) :: this
225
226	integer :: jgas
227
228	if (this%is_sw) then
229	write(nulout,'(a)',advance='no') 'ecCKD shortwave gas optics model: '
230	else
231	write(nulout,'(a)',advance='no') 'ecCKD longwave gas optics model: '
232	end if
233
234	write(nulout,'(i0,a,i0,a,i0,a,i0,a)') &
235	& nint(this%spectral_def%wavenumber1(1)), '-', &
236	& nint(this%spectral_def%wavenumber2(size(this%spectral_def%wavenumber2))), &
237	& ' cm-1, ', this%ng, ' g-points in ', this%spectral_def%nband, ' bands'
238	write(nulout,'(a,i0,a,i0,a,i0,a)') ' Look-up table sizes: ', this%npress, ' pressures, ', &
239	& this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries'
240	write(nulout, '(a)') ' Gases:'
241	do jgas = 1,this%ngas
242	if (this%single_gas(jgas)%i_gas_code > 0) then
243	write(nulout, '(a,a,a)', advance='no') ' ', &
244	& trim(GasName(this%single_gas(jgas)%i_gas_code)), ': '
245	else
246	write(nulout, '(a)', advance='no') ' Composite of well-mixed background gases: '
247	end if
248	select case (this%single_gas(jgas)%i_conc_dependence)
249	case (IConcDependenceNone)
250	write(nulout, '(a)') 'no concentration dependence'
251	case (IConcDependenceLinear)
252	write(nulout, '(a)') 'linear concentration dependence'
253	case (IConcDependenceRelativeLinear)
254	write(nulout, '(a,e14.6)') 'linear concentration dependence relative to a mole fraction of ', &
255	& this%single_gas(jgas)%reference_mole_frac
256	case (IConcDependenceLUT)
257	write(nulout, '(a,i0,a,e14.6,a,e13.6)') 'look-up table with ', this%single_gas(jgas)%n_mole_frac, &
258	& ' log-spaced mole fractions in range ', exp(this%single_gas(jgas)%log_mole_frac1), &
259	& ' to ', exp(this%single_gas(jgas)%log_mole_frac1 &
260	& + this%single_gas(jgas)%n_mole_frac*this%single_gas(jgas)%d_log_mole_frac)
261	end select
262	end do
263
264	end subroutine print_ckd_model
265
266
267	!---------------------------------------------------------------------
268	! Compute layerwise optical depth for each g point for ncol columns
269	! at nlev layers
270	subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
271	& pressure_hl, temperature_fl, mole_fraction_fl, &
272	& optical_depth_fl, rayleigh_od_fl)
273
274	use yomhook, only : lhook, dr_hook
275	use radiation_constants, only : AccelDueToGravity
276
277	! Input variables
278
279	class(ckd_model_type), intent(in), target :: this
280	! Number of columns, levels and input gases
281	integer, intent(in) :: ncol, nlev, nmaxgas, istartcol, iendcol
282	! Pressure at half levels (Pa), dimensioned (ncol,nlev+1)
283	real(jprb), intent(in) :: pressure_hl(ncol,nlev+1)
284	! Temperature at full levels (K), dimensioned (ncol,nlev)
285	real(jprb), intent(in) :: temperature_fl(istartcol:iendcol,nlev)
286	! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
287	real(jprb), intent(in) :: mole_fraction_fl(ncol,nlev,nmaxgas)
288
289	! Output variables
290
291	! Layer absorption optical depth for each g point
292	real(jprb), intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
293	! In the shortwave only, the Rayleigh scattering optical depth
294	real(jprb), optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
295
296	! Local variables
297
298	real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:)
299
300	! Natural logarithm of pressure at full levels
301	real(jprb) :: log_pressure_fl(nlev)
302
303	! Optical depth of single gas at one point in space versus
304	! spectral interval
305	!real(jprb) :: od_single_gas(this%ng)
306
307	real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1
308
309	! Indices and weights in temperature, pressure and concentration interpolation
310	real(jprb) :: pindex1, tindex1, cindex1
311	real(jprb) :: pw1(nlev), pw2(nlev), tw1(nlev), tw2(nlev), cw1(nlev), cw2(nlev)
312	integer :: ip1(nlev), it1(nlev), ic1(nlev)
313
314	! Natural logarithm of mole fraction at one point
315	real(jprb) :: log_conc
316
317	! Minimum mole fraction in look-up-table
318	real(jprb) :: mole_frac1
319
320	integer :: jcol, jlev, jgas, igascode
321
322	real(jprb) :: hook_handle
323
324	if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',0,hook_handle)
325
326	global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
327
328	do jcol = istartcol,iendcol
329
330	log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)))
331
332	do jlev = 1,nlev
333	! Find interpolation points in pressure
334	pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) &
335	& / this%d_log_pressure
336	pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb))
337	ip1(jlev) = int(pindex1)
338	pw2(jlev) = pindex1 - ip1(jlev)
339	pw1(jlev) = 1.0_jprb - pw2(jlev)
340
341	! Find interpolation points in temperature
342	temperature1 = pw1(jlev)*this%temperature1(ip1(jlev)) &
343	& + pw2(jlev)*this%temperature1(ip1(jlev)+1)
344	tindex1 = (temperature_fl(jcol,jlev) - temperature1) &
345	& / this%d_temperature
346	tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb))
347	it1(jlev) = int(tindex1)
348	tw2(jlev) = tindex1 - it1(jlev)
349	tw1(jlev) = 1.0_jprb - tw2(jlev)
350
351	! Concentration multiplier
352	simple_multiplier(jlev) = global_multiplier &
353	& * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev))
354	end do
355
356	optical_depth_fl(:,:,jcol) = 0.0_jprb
357
358	do jgas = 1,this%ngas
359
360	associate (single_gas => this%single_gas(jgas))
361	igascode = this%single_gas(jgas)%i_gas_code
362
363	select case (single_gas%i_conc_dependence)
364
365	case (IConcDependenceLinear)
366	molar_abs => this%single_gas(jgas)%molar_abs
367	multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode)
368
369	do jlev = 1,nlev
370	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
371	& + (multiplier(jlev)tw1(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
372	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
373	& + (multiplier(jlev)tw2(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
374	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
375	end do
376
377	case (IConcDependenceRelativeLinear)
378	molar_abs => this%single_gas(jgas)%molar_abs
379	multiplier = simple_multiplier * (mole_fraction_fl(jcol,:,igascode) &
380	& - single_gas%reference_mole_frac)
381	do jlev = 1,nlev
382	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
383	& + (multiplier(jlev)tw1(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
384	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
385	& + (multiplier(jlev)tw2(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
386	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
387	end do
388
389	case (IConcDependenceNone)
390	! Composite gases
391	molar_abs => this%single_gas(jgas)%molar_abs
392	do jlev = 1,nlev
393	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
394	& + (simple_multiplier(jlev)tw1(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
395	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
396	& + (simple_multiplier(jlev)tw2(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
397	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
398	end do
399
400	case (IConcDependenceLUT)
401	! Logarithmic interpolation in concentration space
402	molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
403	mole_frac1 = exp(single_gas%log_mole_frac1)
404	do jlev = 1,nlev
405	! Take care of mole_fraction == 0
406	log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
407	cindex1 = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
408	cindex1 = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
409	ic1(jlev) = int(cindex1)
410	cw2(jlev) = cindex1 - ic1(jlev)
411	cw1(jlev) = 1.0_jprb - cw2(jlev)
412	end do
413	! od_single_gas = cw1 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1) &
414	! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1)) &
415	! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1) &
416	! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1))) &
417	! & +cw2 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1+1) &
418	! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1+1)) &
419	! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) &
420	! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1)))
421	do jlev = 1,nlev
422	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
423	& + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( &
424	& (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) &
425	& +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) &
426	& +(cw1(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)) &
427	& +(cw1(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)) &
428	& +(cw2(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)+1) &
429	& +(cw2(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)+1) &
430	& +(cw2(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)+1) &
431	& +(cw2(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)+1))
432	end do
433	end select
434
435	end associate
436
437	end do
438
439	! Ensure the optical depth is not negative
440	optical_depth_fl(:,:,jcol) = max(0.0_jprb, optical_depth_fl(:,:,jcol))
441
442	! Rayleigh scattering
443	if (this%is_sw .and. present(rayleigh_od_fl)) then
444	do jlev = 1,nlev
445	rayleigh_od_fl(:,jlev,jcol) = global_multiplier &
446	& * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat
447	end do
448	end if
449
450	end do
451
452	if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',1,hook_handle)
453
454	end subroutine calc_optical_depth_ckd_model
455
456	!---------------------------------------------------------------------
457	! Calculate the Planck function integrated across each of the g
458	! points of this correlated k-distribution model, for a given
459	! temperature, where Planck function is defined as the flux emitted
460	! by a black body (rather than radiance)
461	subroutine calc_planck_function(this, nt, temperature, planck)
462
463	class(ckd_model_type), intent(in) :: this
464	integer, intent(in) :: nt
465	real(jprb), intent(in) :: temperature(:) ! K
466	real(jprb), intent(out) :: planck(this%ng,nt) ! W m-2
467
468	real(jprb) :: tindex1, tw1, tw2
469	integer :: it1, jt
470
471	do jt = 1,nt
472	tindex1 = (temperature(jt) - this%temperature1_planck) &
473	& * (1.0_jprb / this%d_temperature_planck)
474	if (tindex1 >= 0) then
475	! Normal interpolation, and extrapolation for high temperatures
476	tindex1 = 1.0_jprb + tindex1
477	it1 = min(int(tindex1), this%nplanck-1)
478	tw2 = tindex1 - it1
479	tw1 = 1.0_jprb - tw2
480	planck(:,jt) = tw1 * this%planck_function(:,it1) &
481	& + tw2 * this%planck_function(:,it1+1)
482	else
483	! Interpolate linearly to zero
484	planck(:,jt) = this%planck_function(:,1) &
485	& * (temperature(jt)/this%temperature1_planck)
486	end if
487	end do
488
489	end subroutine calc_planck_function
490
491
492	!---------------------------------------------------------------------
493	! Return the spectral solar irradiance integrated over each g point
494	! of a solar correlated k-distribution model, given the
495	! total_solar_irradiance
496	subroutine calc_incoming_sw(this, total_solar_irradiance, &
497	& spectral_solar_irradiance)
498
499	class(ckd_model_type), intent(in) :: this
500	real(jprb), intent(in) :: total_solar_irradiance ! W m-2
501	real(jprb), intent(inout) :: spectral_solar_irradiance(:,:) ! W m-2
502
503	spectral_solar_irradiance &
504	& = spread(total_solar_irradiance * this%norm_solar_irradiance, &
505	& 2, size(spectral_solar_irradiance,2))
506
507	end subroutine calc_incoming_sw
508
509	end module radiation_ecckd

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

Download in other formats:

Original Format