Context Navigation

radiation_ecckd.F90 @ 4773

Last change on this file since 4773 was 4773, checked in by idelkadi, 5 months ago

Update of Ecrad in LMDZ The same organization of the Ecrad offline version is retained in order to facilitate the updating of Ecrad in LMDZ and the comparison between online and offline results. version 1.6.1 of Ecrad (https://github.com/lguez/ecrad.git)

Implementation of the double call of Ecrad in LMDZ

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

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

Download in other formats:

Original Format