source: LMDZ6/branches/contrails/libf/phylmd/ecrad/radiation/radiation_interface.F90 @ 5456

Last change on this file since 5456 was 4853, checked in by idelkadi, 9 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

  • version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
  • interface routines between lmdz and ecrad grouped in a new "lmdz" directory
  • double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
  • clean-up in the read.F routine: delete unitules arguments
  • modification of default gas model in namelist (default: ECCKD)
File size: 28.9 KB
Line 
1! radiation_interface.F90 - Public interface to radiation scheme
2!
3! (C) Copyright 2014- ECMWF.
4!
5! This software is licensed under the terms of the Apache Licence Version 2.0
6! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
7!
8! In applying this licence, ECMWF does not waive the privileges and immunities
9! granted to it by virtue of its status as an intergovernmental organisation
10! nor does it submit to any jurisdiction.
11!
12! Author:  Robin Hogan
13! Email:   r.j.hogan@ecmwf.int
14!
15! Modifications
16!   2017-04-11  R. Hogan  Changes to enable generalized surface description
17!   2017-09-08  R. Hogan  Reverted some changes
18!
19! To use the radiation scheme, create a configuration_type object,
20! call "setup_radiation" on it once to load the look-up-tables and
21! data describing how gas and hydrometeor absorption/scattering are to
22! be represented, and call "radiation" multiple times on different
23! input profiles.
24
25module radiation_interface
26
27  implicit none
28
29  public  :: setup_radiation, set_gas_units, radiation
30  private :: radiation_reverse
31
32contains
33
34  !---------------------------------------------------------------------
35  ! Load the look-up-tables and data describing how gas and
36  ! hydrometeor absorption/scattering are to be represented
37  subroutine setup_radiation(config)
38
39    use parkind1,         only : jprb
40    use yomhook,          only : lhook, dr_hook, jphook
41    use radiation_io,     only : nulerr, radiation_abort
42    use radiation_config, only : config_type, ISolverMcICA, &
43         &   IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD
44    use radiation_spectral_definition, only &
45         &  : SolarReferenceTemperature, TerrestrialReferenceTemperature
46    ! Currently there are two gas absorption models: RRTMG (default)
47    ! and monochromatic
48    use radiation_monochromatic,  only : &
49         &   setup_gas_optics_mono     => setup_gas_optics, &
50         &   setup_cloud_optics_mono   => setup_cloud_optics, &
51         &   setup_aerosol_optics_mono => setup_aerosol_optics
52    use radiation_ifs_rrtm,       only :  setup_gas_optics_rrtmg => setup_gas_optics
53    use radiation_ecckd_interface,only :  setup_gas_optics_ecckd => setup_gas_optics
54    use radiation_cloud_optics,   only :  setup_cloud_optics
55    use radiation_general_cloud_optics, only :  setup_general_cloud_optics
56    use radiation_aerosol_optics, only :  setup_aerosol_optics
57
58   
59    type(config_type), intent(inout) :: config
60
61    real(jphook) :: hook_handle
62
63    if (lhook) call dr_hook('radiation_interface:setup_radiation',0,hook_handle)
64
65    ! Consolidate configuration data, including setting data file
66    ! names
67    call config%consolidate()
68
69    ! Load the look-up tables from files in the specified directory
70    if (config%i_gas_model_sw == IGasModelMonochromatic) then
71      call setup_gas_optics_mono(config, trim(config%directory_name))
72    else
73      ! Note that we can run RRTMG and ECCKD for different parts of
74      ! the spectrum: the setup routines only configure the relevant
75      ! part.
76      if (config%i_gas_model_sw == IGasModelIFSRRTMG .or. config%i_gas_model_lw == IGasModelIFSRRTMG) then
77        call setup_gas_optics_rrtmg(config, trim(config%directory_name))
78      end if
79      if (config%i_gas_model_sw == IGasModelECCKD .or. config%i_gas_model_lw == IGasModelECCKD) then
80        call setup_gas_optics_ecckd(config)
81      end if
82    end if
83
84    if (config%do_lw_aerosol_scattering &
85         & .and. .not. config%do_lw_cloud_scattering) then
86      write(nulerr, '(a)') '*** Error: longwave aerosol scattering requires longwave cloud scattering'
87      call radiation_abort('Radiation configuration error')
88    end if
89
90   
91    ! Whether or not the "radiation" subroutine needs ssa_lw and g_lw
92    ! arrays depends on whether longwave scattering by aerosols is to
93    ! be included.  If not, one of the array dimensions will be set to
94    ! zero.
95    if (config%do_lw_aerosol_scattering) then
96      config%n_g_lw_if_scattering = config%n_g_lw
97    else
98      config%n_g_lw_if_scattering = 0
99    end if
100
101    ! Whether or not the "radiation" subroutine needs ssa_lw_cloud and
102    ! g_lw_cloud arrays depends on whether longwave scattering by
103    ! hydrometeors is to be included.  If not, one of the array
104    ! dimensions will be set to zero.
105    if (config%do_lw_cloud_scattering) then
106      config%n_bands_lw_if_scattering = config%n_bands_lw
107    else
108      config%n_bands_lw_if_scattering = 0
109    end if
110
111    ! If we have longwave scattering and McICA then even if there is
112    ! no aerosol, it is convenient if single scattering albedo and
113    ! g factor arrays are allocated before the call to
114    ! solver_lw as they will be needed.
115    if (config%do_lw_cloud_scattering &
116         &  .and. config%i_solver_lw == ISolverMcICA) then
117      config%n_g_lw_if_scattering = config%n_g_lw
118    end if
119
120    ! Consolidate the albedo/emissivity intervals with the shortwave
121    ! and longwave spectral bands
122    if (config%do_sw) then
123      call config%consolidate_sw_albedo_intervals
124    end if
125    if (config%do_lw) then
126      call config%consolidate_lw_emiss_intervals
127    end if
128
129    if (config%do_clouds) then
130      if (config%i_gas_model_sw == IGasModelMonochromatic) then
131        !      call setup_cloud_optics_mono(config)
132      elseif (config%use_general_cloud_optics) then
133        call setup_general_cloud_optics(config)
134      else
135        call setup_cloud_optics(config)
136      end if
137    end if
138
139    if (config%use_aerosols) then
140      if (config%i_gas_model_sw == IGasModelMonochromatic) then
141!        call setup_aerosol_optics_mono(config)
142      else
143        call setup_aerosol_optics(config)
144      end if
145    end if
146
147    ! Load cloud water PDF look-up table for McICA
148    if (         config%i_solver_sw == ISolverMcICA &
149         &  .or. config%i_solver_lw == ISolverMcICA) then
150      call config%pdf_sampler%setup(config%cloud_pdf_file_name, &
151           &                        iverbose=config%iverbosesetup)
152    end if
153
154    if (lhook) call dr_hook('radiation_interface:setup_radiation',1,hook_handle)
155
156  end subroutine setup_radiation
157
158
159  !---------------------------------------------------------------------
160  ! Scale the gas mixing ratios so that they have the units (and
161  ! possibly scale factors) required by the specific gas absorption
162  ! model.  This subroutine simply passes the gas object on to the
163  ! module of the currently active gas model.
164  subroutine set_gas_units(config, gas)
165   
166    use radiation_config
167    use radiation_gas,             only : gas_type
168    use radiation_monochromatic,   only : set_gas_units_mono  => set_gas_units
169    use radiation_ifs_rrtm,        only : set_gas_units_ifs   => set_gas_units
170    use radiation_ecckd_interface, only : set_gas_units_ecckd => set_gas_units
171
172    type(config_type), intent(in)    :: config
173    type(gas_type),    intent(inout) :: gas
174
175    if (config%i_gas_model_sw == IGasModelMonochromatic) then
176      call set_gas_units_mono(gas)
177    elseif (config%i_gas_model_sw == IGasModelIFSRRTMG &
178         &  .or. config%i_gas_model_lw == IGasModelIFSRRTMG) then
179      ! Convert to mass-mixing ratio for RRTMG; note that ecCKD can
180      ! work with this but performs an internal scaling
181      call set_gas_units_ifs(gas)
182    else
183      ! Use volume mixing ratio preferred by ecCKD
184      call set_gas_units_ecckd(gas)
185    end if
186
187  end subroutine set_gas_units
188
189
190  !---------------------------------------------------------------------
191  ! Run the radiation scheme according to the configuration in the
192  ! config object. There are ncol profiles of which only istartcol to
193  ! iendcol are to be processed, and there are nlev model levels.  The
194  ! output fluxes are written to the flux object, and all other
195  ! objects contain the input variables.  The variables may be defined
196  ! either in order of increasing or decreasing pressure, but if in
197  ! order of decreasing pressure then radiation_reverse will be called
198  ! to reverse the order for the computation and then reverse the
199  ! order of the output fluxes to match the inputs.
200  subroutine radiation(ncol, nlev, istartcol, iendcol, config, &
201       &  single_level, thermodynamics, gas, cloud, aerosol, flux)
202
203    use parkind1,                 only : jprb
204    use yomhook,                  only : lhook, dr_hook, jphook
205
206    use radiation_io,             only : nulout
207    use radiation_config,         only : config_type, &
208         &   IGasModelMonochromatic, IGasModelIFSRRTMG, IGasModelECCKD, &
209         &   ISolverMcICA, ISolverSpartacus, ISolverHomogeneous, &
210         &   ISolverTripleclouds
211    use radiation_single_level,   only : single_level_type
212    use radiation_thermodynamics, only : thermodynamics_type
213    use radiation_gas,            only : gas_type
214    use radiation_cloud,          only : cloud_type
215    use radiation_aerosol,        only : aerosol_type
216    use radiation_flux,           only : flux_type
217    use radiation_spartacus_sw,   only : solver_spartacus_sw
218    use radiation_spartacus_lw,   only : solver_spartacus_lw
219    use radiation_tripleclouds_sw,only : solver_tripleclouds_sw
220    use radiation_tripleclouds_lw,only : solver_tripleclouds_lw
221    use radiation_mcica_sw,       only : solver_mcica_sw
222    use radiation_mcica_lw,       only : solver_mcica_lw
223    use radiation_cloudless_sw,   only : solver_cloudless_sw
224    use radiation_cloudless_lw,   only : solver_cloudless_lw
225    use radiation_homogeneous_sw, only : solver_homogeneous_sw
226    use radiation_homogeneous_lw, only : solver_homogeneous_lw
227    use radiation_save,           only : save_radiative_properties
228
229    ! Treatment of gas and hydrometeor optics
230    use radiation_monochromatic,  only : &
231         &   gas_optics_mono         => gas_optics, &
232         &   cloud_optics_mono       => cloud_optics, &
233         &   add_aerosol_optics_mono => add_aerosol_optics
234    use radiation_ifs_rrtm,       only : gas_optics_rrtmg => gas_optics
235    use radiation_ecckd_interface,only : gas_optics_ecckd => gas_optics
236    use radiation_cloud_optics,   only : cloud_optics
237    use radiation_general_cloud_optics, only : general_cloud_optics
238    use radiation_aerosol_optics, only : add_aerosol_optics
239
240    ! Inputs
241    integer, intent(in) :: ncol               ! number of columns
242    integer, intent(in) :: nlev               ! number of model levels
243    integer, intent(in) :: istartcol, iendcol ! range of columns to process
244    type(config_type),        intent(in)   :: config
245    type(single_level_type),  intent(in)   :: single_level
246    type(thermodynamics_type),intent(in)   :: thermodynamics
247    type(gas_type),           intent(in)   :: gas
248    type(cloud_type),         intent(inout):: cloud
249    type(aerosol_type),       intent(in)   :: aerosol
250    ! Output
251    type(flux_type),          intent(inout):: flux
252
253
254    ! Local variables
255
256    ! Layer optical depth, single scattering albedo and asymmetry factor of
257    ! gases and aerosols at each longwave g-point, where the latter
258    ! two variables are only defined if aerosol longwave scattering is
259    ! enabled (otherwise both are treated as zero).
260    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od_lw
261    real(jprb), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: &
262         &  ssa_lw, g_lw
263
264    ! Layer in-cloud optical depth, single scattering albedo and
265    ! asymmetry factor of hydrometeors in each longwave band, where
266    ! the latter two variables are only defined if hydrometeor
267    ! longwave scattering is enabled (otherwise both are treated as
268    ! zero).
269    real(jprb), dimension(config%n_bands_lw,nlev,istartcol:iendcol) :: od_lw_cloud
270    real(jprb), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol) :: &
271         &  ssa_lw_cloud, g_lw_cloud
272
273    ! Layer optical depth, single scattering albedo and asymmetry factor of
274    ! gases and aerosols at each shortwave g-point
275    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: od_sw, ssa_sw, g_sw
276
277    ! Layer in-cloud optical depth, single scattering albedo and
278    ! asymmetry factor of hydrometeors in each shortwave band
279    real(jprb), dimension(config%n_bands_sw,nlev,istartcol:iendcol)   :: &
280         &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud
281
282    ! The Planck function (emitted flux from a black body) at half
283    ! levels
284    real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: planck_hl
285
286    ! The longwave emission from and albedo of the surface in each
287    ! longwave g-point; note that these are weighted averages of the
288    ! values from individual tiles
289    real(jprb), dimension(config%n_g_lw, istartcol:iendcol) :: lw_emission
290    real(jprb), dimension(config%n_g_lw, istartcol:iendcol) :: lw_albedo
291
292    ! Direct and diffuse shortwave surface albedo in each shortwave
293    ! g-point; note that these are weighted averages of the values
294    ! from individual tiles
295    real(jprb), dimension(config%n_g_sw, istartcol:iendcol) :: sw_albedo_direct
296    real(jprb), dimension(config%n_g_sw, istartcol:iendcol) :: sw_albedo_diffuse
297
298    ! The incoming shortwave flux into a plane perpendicular to the
299    ! incoming radiation at top-of-atmosphere in each of the shortwave
300    ! g-points
301    real(jprb), dimension(config%n_g_sw,istartcol:iendcol) :: incoming_sw
302
303    character(len=100) :: rad_prop_file_name
304    character(*), parameter :: rad_prop_base_file_name = "radiative_properties"
305
306    real(jphook) :: hook_handle
307
308    if (lhook) call dr_hook('radiation_interface:radiation',0,hook_handle)
309
310    if (thermodynamics%pressure_hl(istartcol,2) &
311         &  < thermodynamics%pressure_hl(istartcol,1)) then
312      ! Input arrays are arranged in order of decreasing pressure /
313      ! increasing height: the following subroutine reverses them,
314      ! call the radiation scheme and then reverses the returned
315      ! fluxes
316      call radiation_reverse(ncol, nlev, istartcol, iendcol, config, &
317           &  single_level, thermodynamics, gas, cloud, aerosol, flux)
318    else
319
320      ! Input arrays arranged in order of increasing pressure /
321      ! decreasing height: progress normally
322
323      ! Extract surface albedos at each gridpoint
324      call single_level%get_albedos(istartcol, iendcol, config, &
325           &                        sw_albedo_direct, sw_albedo_diffuse, &
326           &                        lw_albedo)
327
328      ! Compute gas absorption optical depth in shortwave and
329      ! longwave, shortwave single scattering albedo (i.e. fraction of
330      ! extinction due to Rayleigh scattering), Planck functions and
331      ! incoming shortwave flux at each g-point, for the specified
332      ! range of atmospheric columns
333      if (config%i_gas_model_sw == IGasModelMonochromatic) then
334        call gas_optics_mono(ncol,nlev,istartcol,iendcol, config, &
335             &  single_level, thermodynamics, gas, lw_albedo, &
336             &  od_lw, od_sw, ssa_sw, &
337             &  planck_hl, lw_emission, incoming_sw)
338      else
339        if (config%i_gas_model_sw == IGasModelIFSRRTMG &
340             &   .or. config%i_gas_model_lw == IGasModelIFSRRTMG) then
341          call gas_optics_rrtmg(ncol,nlev,istartcol,iendcol, config, &
342               &  single_level, thermodynamics, gas, &
343               &  od_lw, od_sw, ssa_sw, lw_albedo=lw_albedo, &
344               &  planck_hl=planck_hl, lw_emission=lw_emission, &
345               &  incoming_sw=incoming_sw)
346        end if
347        if (config%i_gas_model_sw == IGasModelECCKD &
348             &   .or. config%i_gas_model_lw == IGasModelECCKD) then
349          call gas_optics_ecckd(ncol,nlev,istartcol,iendcol, config, &
350               &  single_level, thermodynamics, gas, &
351               &  od_lw, od_sw, ssa_sw, lw_albedo=lw_albedo, &
352               &  planck_hl=planck_hl, lw_emission=lw_emission, &
353               &  incoming_sw=incoming_sw)
354        end if
355      end if
356
357      if (config%do_clouds) then
358        ! Crop the cloud fraction to remove clouds that have too small
359        ! a fraction or water content; after this, we can safely
360        ! assume that a cloud is present if cloud%fraction > 0.0.
361        call cloud%crop_cloud_fraction(istartcol, iendcol, &
362             &            config%cloud_fraction_threshold, &
363             &            config%cloud_mixing_ratio_threshold)
364
365        ! Compute hydrometeor absorption/scattering properties in each
366        ! shortwave and longwave band
367        if (config%i_gas_model_sw == IGasModelMonochromatic) then
368          call cloud_optics_mono(nlev, istartcol, iendcol, &
369               &  config, thermodynamics, cloud, &
370               &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
371               &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
372        elseif (config%use_general_cloud_optics) then
373          call general_cloud_optics(nlev, istartcol, iendcol, &
374               &  config, thermodynamics, cloud, &
375               &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
376               &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
377        else
378          call cloud_optics(nlev, istartcol, iendcol, &
379               &  config, thermodynamics, cloud, &
380               &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
381               &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
382        end if
383      end if ! do_clouds
384
385      if (config%use_aerosols) then
386        if (config%i_gas_model_sw == IGasModelMonochromatic) then
387!          call add_aerosol_optics_mono(nlev,istartcol,iendcol, &
388!               &  config, thermodynamics, gas, aerosol, &
389!               &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw)
390        else
391          call add_aerosol_optics(nlev,istartcol,iendcol, &
392               &  config, thermodynamics, gas, aerosol, &
393               &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw)
394        end if
395      else
396        g_sw(:,:,istartcol:iendcol) = 0.0_jprb
397        if (config%do_lw_aerosol_scattering) then
398          ssa_lw(:,:,istartcol:iendcol) = 0.0_jprb
399          g_lw(:,:,istartcol:iendcol)   = 0.0_jprb
400        end if
401      end if
402
403      ! For diagnostic purposes, save these intermediate variables to
404      ! a NetCDF file
405      if (config%do_save_radiative_properties) then
406        if (istartcol == 1 .and. iendcol == ncol) then
407          rad_prop_file_name = rad_prop_base_file_name // ".nc"
408        else
409          write(rad_prop_file_name,'(a,a,i4.4,a,i4.4,a)') &
410               &  rad_prop_base_file_name, '_', istartcol, '-',iendcol,'.nc'
411        end if
412        call save_radiative_properties(trim(rad_prop_file_name), &
413             &  nlev, istartcol, iendcol, &
414             &  config, single_level, thermodynamics, cloud, &
415             &  planck_hl, lw_emission, lw_albedo, &
416             &  sw_albedo_direct, sw_albedo_diffuse, incoming_sw, &
417             &  od_lw, ssa_lw, g_lw, od_sw, ssa_sw, g_sw, &
418             &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
419             &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
420      end if
421
422      if (config%do_lw) then
423        if (config%iverbose >= 2) then
424          write(nulout,'(a)') 'Computing longwave fluxes'
425        end if
426
427        if (config%i_solver_lw == ISolverMcICA) then
428          ! Compute fluxes using the McICA longwave solver
429          call solver_mcica_lw(nlev,istartcol,iendcol, &
430               &  config, single_level, cloud, &
431               &  od_lw, ssa_lw, g_lw, od_lw_cloud, ssa_lw_cloud, &
432               &  g_lw_cloud, planck_hl, lw_emission, lw_albedo, flux)
433        else if (config%i_solver_lw == ISolverSPARTACUS) then
434          ! Compute fluxes using the SPARTACUS longwave solver
435          call solver_spartacus_lw(nlev,istartcol,iendcol, &
436               &  config, thermodynamics, cloud, &
437               &  od_lw, ssa_lw, g_lw, od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
438               &  planck_hl, lw_emission, lw_albedo, flux)
439        else if (config%i_solver_lw == ISolverTripleclouds) then
440          ! Compute fluxes using the Tripleclouds longwave solver
441          call solver_tripleclouds_lw(nlev,istartcol,iendcol, &
442               &  config, cloud, &
443               &  od_lw, ssa_lw, g_lw, od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
444               &  planck_hl, lw_emission, lw_albedo, flux)
445        elseif (config%i_solver_lw == ISolverHomogeneous) then
446          ! Compute fluxes using the homogeneous solver
447          call solver_homogeneous_lw(nlev,istartcol,iendcol, &
448               &  config, cloud, &
449               &  od_lw, ssa_lw, g_lw, od_lw_cloud, ssa_lw_cloud, &
450               &  g_lw_cloud, planck_hl, lw_emission, lw_albedo, flux)
451        else
452          ! Compute fluxes using the cloudless solver
453          call solver_cloudless_lw(nlev,istartcol,iendcol, &
454               &  config, od_lw, ssa_lw, g_lw, &
455               &  planck_hl, lw_emission, lw_albedo, flux)
456        end if
457      end if
458
459      if (config%do_sw) then
460        if (config%iverbose >= 2) then
461          write(nulout,'(a)') 'Computing shortwave fluxes'
462        end if
463
464        if (config%i_solver_sw == ISolverMcICA) then
465          ! Compute fluxes using the McICA shortwave solver
466          call solver_mcica_sw(nlev,istartcol,iendcol, &
467               &  config, single_level, cloud, &
468               &  od_sw, ssa_sw, g_sw, od_sw_cloud, ssa_sw_cloud, &
469               &  g_sw_cloud, sw_albedo_direct, sw_albedo_diffuse, &
470               &  incoming_sw, flux)
471        else if (config%i_solver_sw == ISolverSPARTACUS) then
472          ! Compute fluxes using the SPARTACUS shortwave solver
473          call solver_spartacus_sw(nlev,istartcol,iendcol, &
474               &  config, single_level, thermodynamics, cloud, &
475               &  od_sw, ssa_sw, g_sw, od_sw_cloud, ssa_sw_cloud, &
476               &  g_sw_cloud, sw_albedo_direct, sw_albedo_diffuse, &
477               &  incoming_sw, flux)
478        else if (config%i_solver_sw == ISolverTripleclouds) then
479          ! Compute fluxes using the Tripleclouds shortwave solver
480          call solver_tripleclouds_sw(nlev,istartcol,iendcol, &
481               &  config, single_level, cloud, &
482               &  od_sw, ssa_sw, g_sw, od_sw_cloud, ssa_sw_cloud, &
483               &  g_sw_cloud, sw_albedo_direct, sw_albedo_diffuse, &
484               &  incoming_sw, flux)
485        elseif (config%i_solver_sw == ISolverHomogeneous) then
486          ! Compute fluxes using the homogeneous solver
487          call solver_homogeneous_sw(nlev,istartcol,iendcol, &
488               &  config, single_level, cloud, &
489               &  od_sw, ssa_sw, g_sw, od_sw_cloud, ssa_sw_cloud, &
490               &  g_sw_cloud, sw_albedo_direct, sw_albedo_diffuse, &
491               &  incoming_sw, flux)
492        else
493          ! Compute fluxes using the cloudless solver
494          call solver_cloudless_sw(nlev,istartcol,iendcol, &
495               &  config, single_level, od_sw, ssa_sw, g_sw, &
496               &  sw_albedo_direct, sw_albedo_diffuse, &
497               &  incoming_sw, flux)
498        end if
499      end if
500
501      ! Store surface downwelling, and TOA, fluxes in bands from
502      ! fluxes in g points
503      call flux%calc_surface_spectral(config, istartcol, iendcol)
504      call flux%calc_toa_spectral    (config, istartcol, iendcol)
505
506    end if
507   
508    if (lhook) call dr_hook('radiation_interface:radiation',1,hook_handle)
509
510  end subroutine radiation
511
512
513  !---------------------------------------------------------------------
514  ! If the input arrays are arranged in order of decreasing pressure /
515  ! increasing height then this subroutine reverses them, calls the
516  ! radiation scheme and then reverses the returned fluxes. Since this
517  ! subroutine calls, and is called by "radiation", it must be in this
518  ! module to avoid circular dependencies.
519  subroutine radiation_reverse(ncol, nlev, istartcol, iendcol, config, &
520       &  single_level, thermodynamics, gas, cloud, aerosol, flux)
521 
522    use parkind1, only : jprb
523
524    use radiation_io,             only : nulout
525    use radiation_config,         only : config_type
526    use radiation_single_level,   only : single_level_type
527    use radiation_thermodynamics, only : thermodynamics_type
528    use radiation_gas,            only : gas_type
529    use radiation_cloud,          only : cloud_type
530    use radiation_aerosol,        only : aerosol_type
531    use radiation_flux,           only : flux_type
532
533    ! Inputs
534    integer, intent(in) :: ncol               ! number of columns
535    integer, intent(in) :: nlev               ! number of model levels
536    integer, intent(in) :: istartcol, iendcol ! range of columns to process
537    type(config_type),        intent(in) :: config
538    type(single_level_type),  intent(in) :: single_level
539    type(thermodynamics_type),intent(in) :: thermodynamics
540    type(gas_type),           intent(in) :: gas
541    type(cloud_type),         intent(in) :: cloud
542    type(aerosol_type),       intent(in) :: aerosol
543    ! Output
544    type(flux_type),          intent(inout):: flux
545
546    ! Reversed data structures
547    type(thermodynamics_type) :: thermodynamics_rev
548    type(gas_type)            :: gas_rev
549    type(cloud_type)          :: cloud_rev
550    type(aerosol_type)        :: aerosol_rev
551    type(flux_type)           :: flux_rev
552
553    ! Start and end levels for aerosol data
554    integer :: istartlev, iendlev
555
556    if (config%iverbose >= 2) then
557      write(nulout,'(a)') 'Reversing arrays to be in order of increasing pressure'
558    end if
559
560    ! Allocate reversed arrays
561    call thermodynamics_rev%allocate(ncol, nlev)
562    call cloud_rev%allocate(ncol, nlev)
563    call flux_rev%allocate(config, istartcol, iendcol, nlev)
564    if (allocated(aerosol%mixing_ratio)) then
565      istartlev = nlev + 1 - aerosol%iendlev
566      iendlev   = nlev + 1 - aerosol%istartlev
567      call aerosol_rev%allocate(ncol, istartlev, iendlev, &
568           &                    config%n_aerosol_types)
569    end if
570
571    ! Fill reversed thermodynamic arrays
572    thermodynamics_rev%pressure_hl(istartcol:iendcol,:) &
573         &  = thermodynamics%pressure_hl(istartcol:iendcol, nlev+1:1:-1)
574    thermodynamics_rev%temperature_hl(istartcol:iendcol,:) &
575         &  = thermodynamics%temperature_hl(istartcol:iendcol, nlev+1:1:-1)
576
577    ! Fill reversed gas arrays
578    call gas%reverse(istartcol, iendcol, gas_rev)
579
580    if (config%do_clouds) then
581      ! Fill reversed cloud arrays
582      cloud_rev%q_liq(istartcol:iendcol,:) &
583           &  = cloud%q_liq(istartcol:iendcol,nlev:1:-1)
584      cloud_rev%re_liq(istartcol:iendcol,:) &
585           &  = cloud%re_liq(istartcol:iendcol,nlev:1:-1)
586      cloud_rev%q_ice(istartcol:iendcol,:) &
587           &  = cloud%q_ice(istartcol:iendcol,nlev:1:-1)
588      cloud_rev%re_ice(istartcol:iendcol,:) &
589           &  = cloud%re_ice(istartcol:iendcol,nlev:1:-1)
590      cloud_rev%fraction(istartcol:iendcol,:) &
591           &  = cloud%fraction(istartcol:iendcol,nlev:1:-1)
592      cloud_rev%overlap_param(istartcol:iendcol,:) &
593           &  = cloud%overlap_param(istartcol:iendcol,nlev-1:1:-1)
594      if (allocated(cloud%fractional_std)) then
595        cloud_rev%fractional_std(istartcol:iendcol,:) &
596             &  = cloud%fractional_std(istartcol:iendcol,nlev:1:-1)
597      else
598        cloud_rev%fractional_std(istartcol:iendcol,:) = 0.0_jprb       
599      end if
600      if (allocated(cloud%inv_cloud_effective_size)) then
601        cloud_rev%inv_cloud_effective_size(istartcol:iendcol,:) &
602             &  = cloud%inv_cloud_effective_size(istartcol:iendcol,nlev:1:-1)
603      else
604        cloud_rev%inv_cloud_effective_size(istartcol:iendcol,:) = 0.0_jprb
605      end if
606    end if
607
608    if (allocated(aerosol%mixing_ratio)) then
609      aerosol_rev%mixing_ratio(:,istartlev:iendlev,:) &
610           &  = aerosol%mixing_ratio(:,aerosol%iendlev:aerosol%istartlev:-1,:)
611    end if
612
613    ! Run radiation scheme on reversed profiles
614    call radiation(ncol, nlev,istartcol,iendcol, &
615         &  config, single_level, thermodynamics_rev, gas_rev, &
616         &  cloud_rev, aerosol_rev, flux_rev)
617
618    ! Reorder fluxes
619    if (allocated(flux%lw_up)) then
620      flux%lw_up(istartcol:iendcol,:) &
621           &  = flux_rev%lw_up(istartcol:iendcol,nlev+1:1:-1)
622      flux%lw_dn(istartcol:iendcol,:) &
623           &  = flux_rev%lw_dn(istartcol:iendcol,nlev+1:1:-1)
624      if (allocated(flux%lw_up_clear)) then
625        flux%lw_up_clear(istartcol:iendcol,:) &
626             &  = flux_rev%lw_up_clear(istartcol:iendcol,nlev+1:1:-1)
627        flux%lw_dn_clear(istartcol:iendcol,:) &
628             &  = flux_rev%lw_dn_clear(istartcol:iendcol,nlev+1:1:-1)
629      end if
630    end if
631    if (allocated(flux%sw_up)) then
632      flux%sw_up(istartcol:iendcol,:) &
633           &  = flux_rev%sw_up(istartcol:iendcol,nlev+1:1:-1)
634      flux%sw_dn(istartcol:iendcol,:) &
635           &  = flux_rev%sw_dn(istartcol:iendcol,nlev+1:1:-1)
636      if (allocated(flux%sw_dn_direct)) then
637        flux%sw_dn_direct(istartcol:iendcol,:) &
638             &  = flux_rev%sw_dn_direct(istartcol:iendcol,nlev+1:1:-1)
639      end if
640      if (allocated(flux%sw_up_clear)) then
641        flux%sw_up_clear(istartcol:iendcol,:) &
642             &  = flux_rev%sw_up_clear(istartcol:iendcol,nlev+1:1:-1)
643        flux%sw_dn_clear(istartcol:iendcol,:) &
644             &  = flux_rev%sw_dn_clear(istartcol:iendcol,nlev+1:1:-1)
645        if (allocated(flux%sw_dn_direct_clear)) then
646          flux%sw_dn_direct_clear(istartcol:iendcol,:) &
647               &  = flux_rev%sw_dn_direct_clear(istartcol:iendcol,nlev+1:1:-1)
648        end if
649      end if
650    end if
651
652    ! Deallocate reversed arrays
653    call thermodynamics_rev%deallocate
654    call gas_rev%deallocate
655    call cloud_rev%deallocate
656    call flux_rev%deallocate
657    if (allocated(aerosol%mixing_ratio)) then
658      call aerosol_rev%deallocate
659    end if
660
661  end subroutine radiation_reverse
662
663end module radiation_interface
Note: See TracBrowser for help on using the repository browser.