Context Navigation

radiation_ecckd.F90 @ 5159

Last change on this file since 5159 was 5159, checked in by abarral, 3 months ago
Put dimensions.h and paramet.h into modules
File size: 40.1 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	#include "ecrad_config.h"
18
19	module radiation_ecckd
20
21	use parkind1, only : jprb
22	use radiation_gas_constants
23	use radiation_ecckd_gas
24	use radiation_spectral_definition, only : spectral_definition_type
25
26	implicit none
27
28	public
29
30	!---------------------------------------------------------------------
31	! This derived type contains all the data needed to describe a
32	! correlated k-distribution gas optics model created using the ecCKD
33	! tool
34	type ckd_model_type
35
36	! Gas information
37
38	! Number of gases
39	integer :: ngas = 0
40	! Array of individual-gas data objects
41	type(ckd_gas_type), allocatable :: single_gas(:)
42	! Mapping from the "gas codes" in the radiation_gas_constants
43	! module to an index to the single_gas array, where zero means gas
44	! not present (or part of a composite gas)
45	integer :: i_gas_mapping(0:NMaxGases)
46
47	! Coordinates of main look-up table for absorption coeffts
48
49	! Number of pressure and temperature points
50	integer :: npress = 0
51	integer :: ntemp = 0
52	! Natural logarithm of first (lowest) pressure (Pa) and increment
53	real(jprb) :: log_pressure1, d_log_pressure
54	! First temperature profile (K), dimensioned (npress)
55	real(jprb), allocatable :: temperature1(:)
56	! Temperature increment (K)
57	real(jprb) :: d_temperature
58
59	! Look-up table for Planck function
60
61	! Number of entries
62	integer :: nplanck = 0
63	! Temperature of first element of look-up table and increment (K)
64	real(jprb), allocatable :: temperature1_planck
65	real(jprb), allocatable :: d_temperature_planck
66	! Planck function (black body flux into a horizontal plane) in W
67	! m-2, dimensioned (ng,nplanck)
68	real(jprb), allocatable :: planck_function(:,:)
69
70	! Normalized solar irradiance in each g point, dimensioned (ng)
71	real(jprb), allocatable :: norm_solar_irradiance(:)
72
73	! Normalized amplitude of variations in the solar irradiance
74	! through the solar cycle in each g point, dimensioned (ng).
75	! Since the user always provides the solar irradiance SI
76	! integrated across the spectrum, this variable must sum to zero:
77	! this ensures that the solar irradiance in each g-point is
78	! SSI=SI*(norm_solar_irradiance +
79	! A*norm_amplitude_solar_irradiance) for any A, where A denotes
80	! the amplitude of deviations from the mean solar spectrum,
81	! typically between -1.0 and 1.0 and provided by
82	! single_level%solar_spectral_multiplier.
83	real(jprb), allocatable :: norm_amplitude_solar_irradiance(:)
84
85	! Rayleigh molar scattering coefficient in m2 mol-1 in each g
86	! point
87	real(jprb), allocatable :: rayleigh_molar_scat(:)
88
89	! ! Spectral mapping of g points
90
91	! ! Number of wavenumber intervals
92	! integer :: nwav = 0
93
94	! Number of k terms / g points
95	integer :: ng = 0
96
97	! Spectral definition describing bands and g points
98	type(spectral_definition_type) :: spectral_def
99
100	! Shortwave: true, longwave: false
101	logical :: is_sw
102
103	contains
104
105	procedure :: read => read_ckd_model
106	procedure :: read_spectral_solar_cycle
107	! Vectorized version of the optical depth look-up performs better on
108	! NEC, but slower on x86
109	#ifdef DWD_VECTOR_OPTIMIZATIONS
110	procedure :: calc_optical_depth => calc_optical_depth_ckd_model_vec
111	#else
112	procedure :: calc_optical_depth => calc_optical_depth_ckd_model
113	#endif
114	procedure :: print => print_ckd_model
115	procedure :: calc_planck_function
116	procedure :: calc_incoming_sw
117	! procedure :: deallocate => deallocate_ckd_model
118
119	end type ckd_model_type
120
121
122	contains
123
124	!---------------------------------------------------------------------
125	! Read a complete ecCKD gas optics model from a NetCDF file
126	! "filename"
127	subroutine read_ckd_model(this, filename, iverbose)
128
129	#ifdef EASY_NETCDF_READ_MPI
130	use easy_netcdf_read_mpi, only : netcdf_file
131	#else
132	use easy_netcdf, only : netcdf_file
133	#endif
134	!use radiation_io, only : nulerr, radiation_abort
135	use yomhook, only : lhook, dr_hook, jphook
136
137	class(ckd_model_type), intent(inout) :: this
138	character(len=*), intent(in) :: filename
139	integer, optional, intent(in) :: iverbose
140
141	type(netcdf_file) :: file
142
143	real(jprb), allocatable :: pressure_lut(:)
144	real(jprb), allocatable :: temperature_full(:,:)
145	real(jprb), allocatable :: temperature_planck(:)
146
147	character(len=512) :: constituent_id
148
149	integer :: iverbose_local
150
151	! Loop counters
152	integer :: jgas, jjgas
153
154	integer :: istart, inext, nchar, i_gas_code
155
156	real(jphook) :: hook_handle
157
158	if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',0,hook_handle)
159
160	if (present(iverbose)) then
161	iverbose_local = iverbose
162	else
163	iverbose_local = 3
164	end if
165
166	call file%open(trim(filename), iverbose=iverbose_local)
167
168	! Read temperature and pressure coordinate variables
169	call file%get('pressure', pressure_lut)
170	this%log_pressure1 = log(pressure_lut(1))
171	this%npress = size(pressure_lut)
172	this%d_log_pressure = log(pressure_lut(2)) - this%log_pressure1
173	call file%get('temperature', temperature_full)
174	allocate(this%temperature1(this%npress));
175	this%temperature1 = temperature_full(:,1)
176	this%d_temperature = temperature_full(1,2)-temperature_full(1,1)
177	this%ntemp = size(temperature_full,2)
178	deallocate(temperature_full)
179
180	! Read Planck function, or solar irradiance and Rayleigh
181	! scattering coefficient
182	if (file%exists('solar_irradiance')) then
183	this%is_sw = .true.
184	call file%get('solar_irradiance', this%norm_solar_irradiance)
185	this%norm_solar_irradiance = this%norm_solar_irradiance &
186	& / sum(this%norm_solar_irradiance)
187	call file%get('rayleigh_molar_scattering_coeff', &
188	& this%rayleigh_molar_scat)
189	else
190	this%is_sw = .false.
191	call file%get('temperature_planck', temperature_planck)
192	this%nplanck = size(temperature_planck)
193	this%temperature1_planck = temperature_planck(1)
194	this%d_temperature_planck = temperature_planck(2) - temperature_planck(1)
195	deallocate(temperature_planck)
196	call file%get('planck_function', this%planck_function)
197	end if
198
199	! Read the spectral definition information into a separate
200	! derived type
201	call this%spectral_def%read(file);
202	this%ng = this%spectral_def%ng
203
204	! Read gases
205	call file%get('n_gases', this%ngas)
206	allocate(this%single_gas(this%ngas))
207	call file%get_global_attribute('constituent_id', constituent_id)
208	nchar = len(trim(constituent_id))
209	istart = 1
210	this%i_gas_mapping = 0
211	DO jgas = 1, this%ngas
212	if (jgas < this%ngas) then
213	inext = istart + scan(constituent_id(istart:nchar), ' ')
214	else
215	inext = nchar+2
216	end if
217	! Find gas code
218	i_gas_code = 0
219	DO jjgas = 1, NMaxGases
220	if (constituent_id(istart:inext-2) == trim(GasLowerCaseName(jjgas))) then
221	i_gas_code = jjgas
222	exit
223	end if
224	end do
225	! if (i_gas_code == 0) then
226	! write(nulerr,'(a,a,a)') '*** Error: Gas "', constituent_id(istart:inext-2), &
227	! & '" not understood'
228	! call radiation_abort('Radiation configuration error')
229	! end if
230	this%i_gas_mapping(i_gas_code) = jgas;
231	call this%single_gas(jgas)%read(file, constituent_id(istart:inext-2), i_gas_code)
232	istart = inext
233	end do
234
235	call file%close()
236
237	if (lhook) call dr_hook('radiation_ecckd:read_ckd_model',1,hook_handle)
238
239	end subroutine read_ckd_model
240
241	!---------------------------------------------------------------------
242	! Print a description of the correlated k-distribution model to the
243	! "nulout" unit
244	subroutine print_ckd_model(this)
245
246	use radiation_io, only : nulout
247	use radiation_gas_constants
248
249	class(ckd_model_type), intent(in) :: this
250
251	integer :: jgas
252
253	if (this%is_sw) then
254	write(nulout,'(a)',advance='no') 'ecCKD shortwave gas optics model: '
255	else
256	write(nulout,'(a)',advance='no') 'ecCKD longwave gas optics model: '
257	end if
258
259	write(nulout,'(i0,a,i0,a,i0,a,i0,a)') &
260	& nint(this%spectral_def%wavenumber1(1)), '-', &
261	& nint(this%spectral_def%wavenumber2(size(this%spectral_def%wavenumber2))), &
262	& ' cm-1, ', this%ng, ' g-points in ', this%spectral_def%nband, ' bands'
263	write(nulout,'(a,i0,a,i0,a,i0,a)') ' Look-up table sizes: ', this%npress, ' pressures, ', &
264	& this%ntemp, ' temperatures, ', this%nplanck, ' planck-function entries'
265	write(nulout, '(a)') ' Gases:'
266	DO jgas = 1,this%ngas
267	if (this%single_gas(jgas)%i_gas_code > 0) then
268	write(nulout, '(a,a,a)', advance='no') ' ', &
269	& trim(GasName(this%single_gas(jgas)%i_gas_code)), ': '
270	else
271	write(nulout, '(a)', advance='no') ' Composite of well-mixed background gases: '
272	end if
273	select case (this%single_gas(jgas)%i_conc_dependence)
274	case (IConcDependenceNone)
275	write(nulout, '(a)') 'no concentration dependence'
276	case (IConcDependenceLinear)
277	write(nulout, '(a)') 'linear concentration dependence'
278	case (IConcDependenceRelativeLinear)
279	write(nulout, '(a,e14.6)') 'linear concentration dependence relative to a mole fraction of ', &
280	& this%single_gas(jgas)%reference_mole_frac
281	case (IConcDependenceLUT)
282	write(nulout, '(a,i0,a,e14.6,a,e13.6)') 'look-up table with ', this%single_gas(jgas)%n_mole_frac, &
283	& ' log-spaced mole fractions in range ', exp(this%single_gas(jgas)%log_mole_frac1), &
284	& ' to ', exp(this%single_gas(jgas)%log_mole_frac1 &
285	& + this%single_gas(jgas)%n_mole_frac*this%single_gas(jgas)%d_log_mole_frac)
286	end select
287	end do
288
289	end subroutine print_ckd_model
290
291
292	!---------------------------------------------------------------------
293	! Read the amplitude of the spectral variations associated with the
294	! solar cycle and map to g-points
295	subroutine read_spectral_solar_cycle(this, filename, iverbose, use_updated_solar_spectrum)
296
297	#ifdef EASY_NETCDF_READ_MPI
298	use easy_netcdf_read_mpi, only : netcdf_file
299	#else
300	use easy_netcdf, only : netcdf_file
301	#endif
302	use radiation_io, only : nulout, nulerr, radiation_abort
303	use yomhook, only : lhook, dr_hook, jphook
304
305	! Reference total solar irradiance (W m-2)
306	real(jprb), parameter :: ReferenceTSI = 1361.0_jprb
307
308	class(ckd_model_type), intent(inout) :: this
309	character(len=*), intent(in) :: filename
310	integer, optional, intent(in) :: iverbose
311	! Do we update the mean solar spectral irradiance for each g-point
312	! based on the contents of the file?
313	logical, optional, intent(in) :: use_updated_solar_spectrum
314
315	type(netcdf_file) :: file
316
317	! Solar spectral irradiance, its amplitude and wavenumber
318	! coordinate variable, read from NetCDF file
319	real(jprb), allocatable :: wavenumber(:) ! cm-1
320	real(jprb), allocatable :: ssi(:) ! W m-2 cm
321	real(jprb), allocatable :: ssi_amplitude(:) ! W m-2 cm
322
323	! As above but on the wavenumber grid delimited by
324	! this%wavenumber1 and this%wavenumber2
325	real(jprb), allocatable :: ssi_grid(:)
326	real(jprb), allocatable :: ssi_amplitude_grid(:)
327	real(jprb), allocatable :: wavenumber_grid(:)
328
329	! Old normalized solar irradiance in case it gets changed and we
330	! need to report the amplitude of the change
331	real(jprb), allocatable :: old_norm_solar_irradiance(:)
332
333	real(jprb) :: dwav_grid
334
335	! Number of input wavenumbers, and number on ecCKD model's grid
336	integer :: nwav, nwav_grid
337	! Corresponding loop indices
338	integer :: jwav, jwav_grid, jg
339
340	integer :: iband
341
342	integer :: iverbose_local
343
344	real(jphook) :: hook_handle
345
346	if (lhook) call dr_hook('radiation_ecckd:read_spectral_solar_cycle',0,hook_handle)
347
348	if (present(iverbose)) then
349	iverbose_local = iverbose
350	else
351	iverbose_local = 3
352	end if
353
354	call file%open(trim(filename), iverbose=iverbose_local)
355
356	call file%get('wavenumber', wavenumber)
357	call file%get('mean_solar_spectral_irradiance', ssi)
358	call file%get('ssi_solar_cycle_amplitude', ssi_amplitude)
359
360	call file%close()
361
362	nwav = size(wavenumber)
363
364	nwav_grid = size(this%spectral_def%wavenumber1)
365	allocate(ssi_grid(nwav_grid))
366	allocate(ssi_amplitude_grid(nwav_grid))
367	allocate(wavenumber_grid(nwav_grid))
368	wavenumber_grid = 0.5_jprb * (this%spectral_def%wavenumber1+this%spectral_def%wavenumber2)
369	dwav_grid = this%spectral_def%wavenumber2(1)-this%spectral_def%wavenumber1(1)
370
371	ssi_grid = 0.0_jprb
372	ssi_amplitude_grid = 0.0_jprb
373
374	! Interpolate input SSI to regular wavenumber grid
375	DO jwav_grid = 1,nwav_grid
376	DO jwav = 1,nwav-1
377	if (wavenumber(jwav) < wavenumber_grid(jwav_grid) &
378	& .and. wavenumber(jwav+1) >= wavenumber_grid(jwav_grid)) then
379	! Linear interpolation - this is not perfect
380	ssi_grid(jwav_grid) = (ssi(jwav)*(wavenumber(jwav+1)-wavenumber_grid(jwav_grid)) &
381	& +ssi(jwav+1)*(wavenumber_grid(jwav_grid)-wavenumber(jwav))) &
382	& * dwav_grid / (wavenumber(jwav+1)-wavenumber(jwav))
383	ssi_amplitude_grid(jwav_grid) = (ssi_amplitude(jwav)*(wavenumber(jwav+1)-wavenumber_grid(jwav_grid)) &
384	& +ssi_amplitude(jwav+1)*(wavenumber_grid(jwav_grid)-wavenumber(jwav))) &
385	& * dwav_grid / (wavenumber(jwav+1)-wavenumber(jwav))
386	exit
387	end if
388	end do
389	end do
390
391	! Optionally update the solar irradiances in each g-point, and the
392	! spectral solar irradiance on the wavenumber grid corresponding
393	! to gpoint_fraction
394	allocate(old_norm_solar_irradiance(nwav_grid))
395	old_norm_solar_irradiance = this%norm_solar_irradiance
396	if (present(use_updated_solar_spectrum)) then
397	if (use_updated_solar_spectrum) then
398	if (.not. allocated(this%spectral_def%solar_spectral_irradiance)) then
399	write(nulerr,'(a)') 'Cannot use_updated_solar_spectrum unless gas optics model is from ecCKD >= 1.4'
400	call radiation_abort()
401	end if
402	this%norm_solar_irradiance = old_norm_solar_irradiance &
403	& * matmul(ssi_grid,this%spectral_def%gpoint_fraction) &
404	& / matmul(this%spectral_def%solar_spectral_irradiance,this%spectral_def%gpoint_fraction)
405	this%norm_solar_irradiance = this%norm_solar_irradiance / sum(this%norm_solar_irradiance)
406	this%spectral_def%solar_spectral_irradiance = ssi_grid
407	end if
408	end if
409
410	! Map on to g-points
411	this%norm_amplitude_solar_irradiance &
412	& = this%norm_solar_irradiance &
413	& * matmul(ssi_amplitude_grid, this%spectral_def%gpoint_fraction) &
414	& / matmul(ssi_grid,this%spectral_def%gpoint_fraction)
415
416	! Remove the mean from the solar-cycle fluctuations, since the
417	! user will scale with total solar irradiance
418	this%norm_amplitude_solar_irradiance &
419	& = (this%norm_solar_irradiance+this%norm_amplitude_solar_irradiance) &
420	& / sum(this%norm_solar_irradiance+this%norm_amplitude_solar_irradiance) &
421	& - this%norm_solar_irradiance
422
423	! Print the spectral solar irradiance per g point, and solar cycle amplitude
424	if (iverbose_local >= 2) then
425	write(nulout,'(a,f6.1,a)') 'G-point, solar irradiance for nominal TSI = ', &
426	& ReferenceTSI, ' W m-2, solar cycle amplitude (at solar maximum), update to original solar irradiance'
427	iband = 0
428	DO jg = 1,this%ng
429	if (this%spectral_def%i_band_number(jg) > iband) then
430	iband = this%spectral_def%i_band_number(jg)
431	write(nulout, '(i2,f10.4,f7.3,a,f8.4,a,i2,a,f7.1,a,f7.1,a)') &
432	& jg, ReferenceTSI*this%norm_solar_irradiance(jg), &
433	& 100.0_jprb * this%norm_amplitude_solar_irradiance(jg) &
434	& / this%norm_solar_irradiance(jg), '% ', &
435	& 100.0_jprb * (this%norm_solar_irradiance(jg) &
436	& / old_norm_solar_irradiance(jg) - 1.0_jprb), '% Band ', iband, ': ', &
437	& this%spectral_def%wavenumber1_band(iband), '-', &
438	& this%spectral_def%wavenumber2_band(iband), ' cm-1'
439	else
440	write(nulout, '(i2,f10.4,f7.3,a,f8.4,a)') jg, ReferenceTSI*this%norm_solar_irradiance(jg), &
441	& 100.0_jprb * this%norm_amplitude_solar_irradiance(jg) &
442	& / this%norm_solar_irradiance(jg), '% ', &
443	& 100.0_jprb * (this%norm_solar_irradiance(jg) &
444	& / old_norm_solar_irradiance(jg) - 1.0_jprb), '%'
445	end if
446	end do
447	end if
448
449	if (lhook) call dr_hook('radiation_ecckd:read_spectral_solar_cycle',1,hook_handle)
450
451	end subroutine read_spectral_solar_cycle
452
453
454	!---------------------------------------------------------------------
455	! Compute layerwise optical depth for each g point for ncol columns
456	! at nlev layers
457	subroutine calc_optical_depth_ckd_model(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
458	& pressure_hl, temperature_fl, mole_fraction_fl, &
459	& optical_depth_fl, rayleigh_od_fl, concentration_scaling)
460
461	use yomhook, only : lhook, dr_hook, jphook
462	use radiation_constants, only : AccelDueToGravity
463
464	! Input variables
465
466	class(ckd_model_type), intent(in), target :: this
467	! Number of columns, levels and input gases
468	integer, intent(in) :: ncol, nlev, nmaxgas, istartcol, iendcol
469	! Pressure at half levels (Pa), dimensioned (ncol,nlev+1)
470	real(jprb), intent(in) :: pressure_hl(ncol,nlev+1)
471	! Temperature at full levels (K), dimensioned (ncol,nlev)
472	real(jprb), intent(in) :: temperature_fl(istartcol:iendcol,nlev)
473	! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
474	real(jprb), intent(in) :: mole_fraction_fl(ncol,nlev,nmaxgas)
475	! Optional concentration scaling of each gas
476	real(jprb), optional, intent(in) :: concentration_scaling(nmaxgas)
477
478	! Output variables
479
480	! Layer absorption optical depth for each g point
481	real(jprb), intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
482	! In the shortwave only, the Rayleigh scattering optical depth
483	real(jprb), optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
484
485	! Local variables
486
487	real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:)
488
489	! Natural logarithm of pressure at full levels
490	real(jprb) :: log_pressure_fl(nlev)
491
492	! Optical depth of single gas at one point in space versus
493	! spectral interval
494	!real(jprb) :: od_single_gas(this%ng)
495
496	real(jprb) :: multiplier(nlev), simple_multiplier(nlev), global_multiplier, temperature1
497	real(jprb) :: scaling
498
499	! Indices and weights in temperature, pressure and concentration interpolation
500	real(jprb) :: pindex1, tindex1, cindex1
501	real(jprb) :: pw1(nlev), pw2(nlev), tw1(nlev), tw2(nlev), cw1(nlev), cw2(nlev)
502	integer :: ip1(nlev), it1(nlev), ic1(nlev)
503
504	! Natural logarithm of mole fraction at one point
505	real(jprb) :: log_conc
506
507	! Minimum mole fraction in look-up-table
508	real(jprb) :: mole_frac1
509
510	integer :: jcol, jlev, jgas, igascode
511
512	real(jphook) :: hook_handle
513
514	if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',0,hook_handle)
515
516	global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
517
518	DO jcol = istartcol,iendcol
519
520	log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,1:nlev)+pressure_hl(jcol,2:nlev+1)))
521
522	DO jlev = 1,nlev
523	! Find interpolation points in pressure
524	pindex1 = (log_pressure_fl(jlev)-this%log_pressure1) &
525	& / this%d_log_pressure
526	pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb))
527	ip1(jlev) = int(pindex1)
528	pw2(jlev) = pindex1 - ip1(jlev)
529	pw1(jlev) = 1.0_jprb - pw2(jlev)
530
531	! Find interpolation points in temperature
532	temperature1 = pw1(jlev)*this%temperature1(ip1(jlev)) &
533	& + pw2(jlev)*this%temperature1(ip1(jlev)+1)
534	tindex1 = (temperature_fl(jcol,jlev) - temperature1) &
535	& / this%d_temperature
536	tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb))
537	it1(jlev) = int(tindex1)
538	tw2(jlev) = tindex1 - it1(jlev)
539	tw1(jlev) = 1.0_jprb - tw2(jlev)
540
541	! Concentration multiplier
542	simple_multiplier(jlev) = global_multiplier &
543	& * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev))
544	end do
545
546	optical_depth_fl(:,:,jcol) = 0.0_jprb
547
548	DO jgas = 1,this%ngas
549
550	associate (single_gas => this%single_gas(jgas))
551	igascode = this%single_gas(jgas)%i_gas_code
552
553	select case (single_gas%i_conc_dependence)
554
555	case (IConcDependenceLinear)
556	molar_abs => this%single_gas(jgas)%molar_abs
557	multiplier = simple_multiplier * mole_fraction_fl(jcol,:,igascode)
558
559	if (present(concentration_scaling)) then
560	multiplier = multiplier * concentration_scaling(igascode)
561	end if
562
563	DO jlev = 1,nlev
564	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
565	& + (multiplier(jlev)tw1(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
566	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
567	& + (multiplier(jlev)tw2(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
568	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
569	end do
570
571	case (IConcDependenceRelativeLinear)
572	molar_abs => this%single_gas(jgas)%molar_abs
573
574	if (present(concentration_scaling)) then
575	multiplier = simple_multiplier &
576	& * (mole_fraction_fl(jcol,:,igascode)*concentration_scaling(igascode) &
577	& - single_gas%reference_mole_frac)
578	else
579	multiplier = simple_multiplier * (mole_fraction_fl(jcol,:,igascode) &
580	& - single_gas%reference_mole_frac)
581	end if
582
583	DO jlev = 1,nlev
584	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
585	& + (multiplier(jlev)tw1(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
586	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
587	& + (multiplier(jlev)tw2(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
588	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
589	end do
590
591	case (IConcDependenceNone)
592	! Composite gases
593	molar_abs => this%single_gas(jgas)%molar_abs
594	DO jlev = 1,nlev
595	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
596	& + (simple_multiplier(jlev)tw1(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)) &
597	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev))) &
598	& + (simple_multiplier(jlev)tw2(jlev)) (pw1(jlev) * molar_abs(:,ip1(jlev),it1(jlev)+1) &
599	& +pw2(jlev) * molar_abs(:,ip1(jlev)+1,it1(jlev)+1))
600	end do
601
602	case (IConcDependenceLUT)
603
604	if (present(concentration_scaling)) then
605	scaling = concentration_scaling(igascode)
606	else
607	scaling = 1.0_jprb
608	end if
609
610	! Logarithmic interpolation in concentration space
611	molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
612	mole_frac1 = exp(single_gas%log_mole_frac1)
613	DO jlev = 1,nlev
614	! Take care of mole_fraction == 0
615	log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode)*scaling, mole_frac1))
616	cindex1 = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
617	cindex1 = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
618	ic1(jlev) = int(cindex1)
619	cw2(jlev) = cindex1 - ic1(jlev)
620	cw1(jlev) = 1.0_jprb - cw2(jlev)
621	end do
622	! od_single_gas = cw1 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1) &
623	! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1)) &
624	! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1) &
625	! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1))) &
626	! & +cw2 * (tw1 * (pw1 * molar_abs_conc(:,ip1,it1,ic1+1) &
627	! & +pw2 * molar_abs_conc(:,ip1+1,it1,ic1+1)) &
628	! & +tw2 * (pw1 * molar_abs_conc(:,ip1,it1+1,ic1+1) &
629	! & +pw2 * molar_abs_conc(:,ip1+1,it1+1,ic1+1)))
630	DO jlev = 1,nlev
631	optical_depth_fl(:,jlev,jcol) = optical_depth_fl(:,jlev,jcol) &
632	& + (simple_multiplier(jlev) * mole_fraction_fl(jcol,jlev,igascode) * scaling) * ( &
633	& (cw1(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)) &
634	& +(cw1(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)) &
635	& +(cw1(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)) &
636	& +(cw1(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)) &
637	& +(cw2(jlev) * tw1(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev),ic1(jlev)+1) &
638	& +(cw2(jlev) * tw1(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev),ic1(jlev)+1) &
639	& +(cw2(jlev) * tw2(jlev) * pw1(jlev)) * molar_abs_conc(:,ip1(jlev),it1(jlev)+1,ic1(jlev)+1) &
640	& +(cw2(jlev) * tw2(jlev) * pw2(jlev)) * molar_abs_conc(:,ip1(jlev)+1,it1(jlev)+1,ic1(jlev)+1))
641	end do
642	end select
643
644	end associate
645
646	end do
647
648	! Ensure the optical depth is not negative
649	optical_depth_fl(:,:,jcol) = max(0.0_jprb, optical_depth_fl(:,:,jcol))
650
651	! Rayleigh scattering
652	if (this%is_sw .and. present(rayleigh_od_fl)) then
653	DO jlev = 1,nlev
654	rayleigh_od_fl(:,jlev,jcol) = global_multiplier &
655	& * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat
656	end do
657	end if
658
659	end do
660
661	if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth',1,hook_handle)
662
663	end subroutine calc_optical_depth_ckd_model
664
665
666	!---------------------------------------------------------------------
667	! Vectorized variant of above routine
668	subroutine calc_optical_depth_ckd_model_vec(this, ncol, nlev, istartcol, iendcol, nmaxgas, &
669	& pressure_hl, temperature_fl, mole_fraction_fl, &
670	& optical_depth_fl, rayleigh_od_fl)
671
672	use yomhook, only : lhook, dr_hook, jphook
673	use radiation_constants, only : AccelDueToGravity
674
675	! Input variables
676
677	class(ckd_model_type), intent(in), target :: this
678	! Number of columns, levels and input gases
679	integer, intent(in) :: ncol, nlev, nmaxgas, istartcol, iendcol
680	! Pressure at half levels (Pa), dimensioned (ncol,nlev+1)
681	real(jprb), intent(in) :: pressure_hl(ncol,nlev+1)
682	! Temperature at full levels (K), dimensioned (ncol,nlev)
683	real(jprb), intent(in) :: temperature_fl(istartcol:iendcol,nlev)
684	! Gas mole fractions at full levels (mol mol-1), dimensioned (ncol,nlev,nmaxgas)
685	real(jprb), intent(in) :: mole_fraction_fl(ncol,nlev,nmaxgas)
686
687	! Output variables
688
689	! Layer absorption optical depth for each g point
690	real(jprb), intent(out) :: optical_depth_fl(this%ng,nlev,istartcol:iendcol)
691	! In the shortwave only, the Rayleigh scattering optical depth
692	real(jprb), optional, intent(out) :: rayleigh_od_fl(this%ng,nlev,istartcol:iendcol)
693
694	! Local variables
695
696	real(jprb), pointer :: molar_abs(:,:,:), molar_abs_conc(:,:,:,:)
697
698	! Natural logarithm of pressure at full levels
699	real(jprb) :: log_pressure_fl
700
701	! Optical depth of single gas at one point in space versus
702	! spectral interval
703	!real(jprb) :: od_single_gas(this%ng)
704
705	real(jprb) :: multiplier, simple_multiplier(ncol,nlev), global_multiplier, temperature1
706
707	! Indices and weights in temperature, pressure and concentration interpolation
708	real(jprb) :: pindex1, tindex1, cindex1
709	real(jprb) :: pw1(ncol,nlev), pw2(ncol,nlev), tw1(ncol,nlev), tw2(ncol,nlev), cw1(ncol,nlev), cw2(ncol,nlev)
710	integer :: ip1(ncol,nlev), it1(ncol,nlev), ic1(ncol,nlev)
711
712	! Natural logarithm of mole fraction at one point
713	real(jprb) :: log_conc
714
715	! Minimum mole fraction in look-up-table
716	real(jprb) :: mole_frac1
717
718	! Layer absorption optical depth for each g point (memory layout adjusted to vectorization)
719	real(jprb) :: od_fl(ncol,this%ng,nlev)
720
721	integer :: jcol, jlev, jgas, igascode, jg
722
723	real(jphook) :: hook_handle
724
725	if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth_vec',0,hook_handle)
726
727	global_multiplier = 1.0_jprb / (AccelDueToGravity * 0.001_jprb * AirMolarMass)
728
729	od_fl(:,:,:) = 0.0_jprb
730
731	DO jlev = 1,nlev
732	DO jcol = istartcol,iendcol
733
734	log_pressure_fl = log(0.5_jprb * (pressure_hl(jcol,jlev)+pressure_hl(jcol,jlev+1)))
735
736	! Find interpolation points in pressure
737	pindex1 = (log_pressure_fl-this%log_pressure1) &
738	& / this%d_log_pressure
739	pindex1 = 1.0_jprb + max(0.0_jprb, min(pindex1, this%npress-1.0001_jprb))
740	ip1(jcol,jlev) = int(pindex1)
741	pw2(jcol,jlev) = pindex1 - ip1(jcol,jlev)
742	pw1(jcol,jlev) = 1.0_jprb - pw2(jcol,jlev)
743
744	! Find interpolation points in temperature
745	temperature1 = pw1(jcol,jlev)*this%temperature1(ip1(jcol,jlev)) &
746	& + pw2(jcol,jlev)*this%temperature1(ip1(jcol,jlev)+1)
747	tindex1 = (temperature_fl(jcol,jlev) - temperature1) &
748	& / this%d_temperature
749	tindex1 = 1.0_jprb + max(0.0_jprb, min(tindex1, this%ntemp-1.0001_jprb))
750	it1(jcol,jlev) = int(tindex1)
751	tw2(jcol,jlev) = tindex1 - it1(jcol,jlev)
752	tw1(jcol,jlev) = 1.0_jprb - tw2(jcol,jlev)
753
754	! Concentration multiplier
755	simple_multiplier(jcol,jlev) = global_multiplier &
756	& * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev))
757	end do
758	end do
759
760	DO jgas = 1,this%ngas
761
762	associate (single_gas => this%single_gas(jgas))
763	igascode = this%single_gas(jgas)%i_gas_code
764
765	select case (single_gas%i_conc_dependence)
766
767	case (IConcDependenceLinear)
768	molar_abs => this%single_gas(jgas)%molar_abs
769
770	DO jlev = 1,nlev
771	DO jg = 1, this%ng
772	DO jcol = istartcol,iendcol
773	multiplier = simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)
774
775	od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
776	& + (multipliertw1(jcol,jlev)) (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) &
777	& +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) &
778	& + (multipliertw2(jcol,jlev)) (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) &
779	& +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1))
780	end do
781	end do
782	end do
783
784	case (IConcDependenceRelativeLinear)
785	molar_abs => this%single_gas(jgas)%molar_abs
786
787	DO jlev = 1,nlev
788	DO jg = 1, this%ng
789	DO jcol = istartcol,iendcol
790	multiplier = simple_multiplier(jcol,jlev) * (mole_fraction_fl(jcol,jlev,igascode) &
791	& - single_gas%reference_mole_frac)
792
793	od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
794	& + (multipliertw1(jcol,jlev)) (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) &
795	& +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) &
796	& + (multipliertw2(jcol,jlev)) (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) &
797	& +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1))
798	end do
799	end do
800	end do
801
802	case (IConcDependenceNone)
803	! Composite gases
804	molar_abs => this%single_gas(jgas)%molar_abs
805
806	DO jlev = 1,nlev
807	DO jg = 1, this%ng
808	DO jcol = istartcol,iendcol
809	od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
810	& + (simple_multiplier(jcol,jlev)tw1(jcol,jlev)) &
811	& (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)) &
812	& +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev))) &
813	& + (simple_multiplier(jcol,jlev)tw2(jcol,jlev)) &
814	& (pw1(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev),it1(jcol,jlev)+1) &
815	& +pw2(jcol,jlev) * molar_abs(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1))
816	end do
817	end do
818	end do
819
820	case (IConcDependenceLUT)
821	! Logarithmic interpolation in concentration space
822	molar_abs_conc => this%single_gas(jgas)%molar_abs_conc
823	mole_frac1 = exp(single_gas%log_mole_frac1)
824
825	DO jlev = 1,nlev
826	DO jcol = istartcol,iendcol
827	! Take care of mole_fraction == 0
828	log_conc = log(max(mole_fraction_fl(jcol,jlev,igascode), mole_frac1))
829	cindex1 = (log_conc - single_gas%log_mole_frac1) / single_gas%d_log_mole_frac
830	cindex1 = 1.0_jprb + max(0.0_jprb, min(cindex1, single_gas%n_mole_frac-1.0001_jprb))
831	ic1(jcol,jlev) = int(cindex1)
832	cw2(jcol,jlev) = cindex1 - ic1(jcol,jlev)
833	cw1(jcol,jlev) = 1.0_jprb - cw2(jcol,jlev)
834	end do
835	end do
836
837	DO jlev = 1,nlev
838	DO jg = 1, this%ng
839	!NEC$ select_vector
840	DO jcol = istartcol,iendcol
841
842	od_fl(jcol,jg,jlev) = od_fl(jcol,jg,jlev) &
843	& + (simple_multiplier(jcol,jlev) * mole_fraction_fl(jcol,jlev,igascode)) * ( &
844	& (cw1(jcol,jlev) * tw1(jcol,jlev) * pw1(jcol,jlev)) * &
845	& molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev),ic1(jcol,jlev)) &
846	& +(cw1(jcol,jlev) * tw1(jcol,jlev) * pw2(jcol,jlev)) * &
847	& molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev),ic1(jcol,jlev)) &
848	& +(cw1(jcol,jlev) * tw2(jcol,jlev) * pw1(jcol,jlev)) * &
849	& molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev)+1,ic1(jcol,jlev)) &
850	& +(cw1(jcol,jlev) * tw2(jcol,jlev) * pw2(jcol,jlev)) * &
851	& molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1,ic1(jcol,jlev)) &
852	& +(cw2(jcol,jlev) * tw1(jcol,jlev) * pw1(jcol,jlev)) * &
853	& molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev),ic1(jcol,jlev)+1) &
854	& +(cw2(jcol,jlev) * tw1(jcol,jlev) * pw2(jcol,jlev)) * &
855	& molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev),ic1(jcol,jlev)+1) &
856	& +(cw2(jcol,jlev) * tw2(jcol,jlev) * pw1(jcol,jlev)) * &
857	& molar_abs_conc(jg,ip1(jcol,jlev),it1(jcol,jlev)+1,ic1(jcol,jlev)+1) &
858	& +(cw2(jcol,jlev) * tw2(jcol,jlev) * pw2(jcol,jlev)) * &
859	& molar_abs_conc(jg,ip1(jcol,jlev)+1,it1(jcol,jlev)+1,ic1(jcol,jlev)+1))
860	end do
861	end do
862	end do
863	end select
864
865	end associate
866
867	! Ensure the optical depth is not negative
868	DO jcol = istartcol,iendcol
869	DO jlev = 1,nlev
870	DO jg = 1, this%ng
871	optical_depth_fl(jg,jlev,jcol) = max(0.0_jprb, od_fl(jcol,jg,jlev))
872	end do
873	end do
874	end do
875
876	! Rayleigh scattering
877	if (this%is_sw .and. present(rayleigh_od_fl)) then
878	DO jcol = istartcol,iendcol
879	DO jlev = 1,nlev
880	DO jg = 1, this%ng
881	rayleigh_od_fl(jg,jlev,jcol) = global_multiplier &
882	& * (pressure_hl(jcol,jlev+1) - pressure_hl(jcol,jlev)) * this%rayleigh_molar_scat(jg)
883	end do
884	end do
885	end do
886	end if
887
888	end do
889
890	if (lhook) call dr_hook('radiation_ecckd:calc_optical_depth_vec',1,hook_handle)
891
892	end subroutine calc_optical_depth_ckd_model_vec
893
894
895	!---------------------------------------------------------------------
896	! Calculate the Planck function integrated across each of the g
897	! points of this correlated k-distribution model, for a given
898	! temperature, where Planck function is defined as the flux emitted
899	! by a black body (rather than radiance)
900	subroutine calc_planck_function(this, nt, temperature, planck)
901
902	class(ckd_model_type), intent(in) :: this
903	integer, intent(in) :: nt
904	real(jprb), intent(in) :: temperature(:) ! K
905	real(jprb), intent(out) :: planck(this%ng,nt) ! W m-2
906
907	real(jprb) :: tindex1, tw1, tw2
908	integer :: it1, jt
909
910	DO jt = 1,nt
911	tindex1 = (temperature(jt) - this%temperature1_planck) &
912	& * (1.0_jprb / this%d_temperature_planck)
913	if (tindex1 >= 0) then
914	! Normal interpolation, and extrapolation for high temperatures
915	tindex1 = 1.0_jprb + tindex1
916	it1 = min(int(tindex1), this%nplanck-1)
917	tw2 = tindex1 - it1
918	tw1 = 1.0_jprb - tw2
919	planck(:,jt) = tw1 * this%planck_function(:,it1) &
920	& + tw2 * this%planck_function(:,it1+1)
921	else
922	! Interpolate linearly to zero
923	planck(:,jt) = this%planck_function(:,1) &
924	& * (temperature(jt)/this%temperature1_planck)
925	end if
926	end do
927
928	end subroutine calc_planck_function
929
930
931	!---------------------------------------------------------------------
932	! Return the spectral solar irradiance integrated over each g point
933	! of a solar correlated k-distribution model, given the
934	! total_solar_irradiance
935	subroutine calc_incoming_sw(this, total_solar_irradiance, &
936	& spectral_solar_irradiance, &
937	& solar_spectral_multiplier)
938
939	use radiation_io, only : nulerr, radiation_abort
940
941	class(ckd_model_type), intent(in) :: this
942	real(jprb), intent(in) :: total_solar_irradiance ! W m-2
943	real(jprb), intent(inout) :: spectral_solar_irradiance(:,:) ! W m-2
944	real(jprb), optional, intent(in) :: solar_spectral_multiplier
945
946	if (.not. present(solar_spectral_multiplier)) then
947	spectral_solar_irradiance &
948	& = spread(total_solar_irradiance * this%norm_solar_irradiance, &
949	& 2, size(spectral_solar_irradiance,2))
950	else if (allocated(this%norm_amplitude_solar_irradiance)) then
951	spectral_solar_irradiance &
952	& = spread(total_solar_irradiance * (this%norm_solar_irradiance &
953	& + solar_spectral_multiplier*this%norm_amplitude_solar_irradiance), &
954	& 2, size(spectral_solar_irradiance,2))
955	else if (solar_spectral_multiplier == 0.0_jprb) then
956	spectral_solar_irradiance &
957	& = spread(total_solar_irradiance * this%norm_solar_irradiance, &
958	& 2, size(spectral_solar_irradiance,2))
959	else
960	write(nulerr, '(a)') '*** Error in calc_incoming_sw: no information present on solar cycle'
961	call radiation_abort()
962	end if
963
964	end subroutine calc_incoming_sw
965
966	end module radiation_ecckd

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

Download in other formats:

Original Format