source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_save.F90 @ 4484

Last change on this file since 4484 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

  • Inclusion of the ecrad directory containing the sources of the ECRAD code
    • interface routine : radiation_scheme.F90
  • Adaptation of compilation scripts :
    • compilation under CPP key CPP_ECRAD
    • compilation with option "-rad ecard" or "-ecard true"
    • The "-rad old/rtm/ecran" build option will need to replace the "-rrtm true" and "-ecrad true" options in the future.
  • Runing LMDZ simulations with ecrad, you need :
    • logical key iflag_rrtm = 2 in physiq.def
    • namelist_ecrad (DefLists?)
    • the directory "data" containing the configuration files is temporarily placed in ../libfphylmd/ecrad/
  • Compilation and execution are tested in the 1D case. The repository under svn would allow to continue the implementation work: tests, verification of the results, ...
File size: 43.0 KB
RevLine 
[3908]1! radiation_save.F90 - Save data to NetCDF files
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-22  R. Hogan  Adapt for new way of describing longwave properties
17!   2019-01-02  R. Hogan  Only save cloud properties if do_clouds==.true.
18
19module radiation_save
20
21  use parkind1, only : jprb
22
23  implicit none
24
25  ! Two available subroutines: save final fluxes and save intermediate
26  ! radiative properties
27  public :: save_fluxes, save_radiative_properties, save_inputs
28
29contains
30
31  !---------------------------------------------------------------------
32  ! Save fluxes in "flux" to NetCDF file_name, plus pressure from the
33  ! thermodynamics object
34  subroutine save_fluxes(file_name, config, thermodynamics, flux, &
35       &                 iverbose, is_hdf5_file, experiment_name)
36
37    use yomhook,                  only : lhook, dr_hook
38
39    use easy_netcdf
40
41    use radiation_io,             only : nulout
42    use radiation_config,         only : config_type, IGasModelMonochromatic
43    use radiation_thermodynamics, only : thermodynamics_type
44    use radiation_flux,           only : flux_type
45
46    character(len=*),           intent(in) :: file_name
47    type(config_type),          intent(in) :: config
48    type(thermodynamics_type),  intent(in) :: thermodynamics
49    type(flux_type),            intent(in) :: flux
50    integer,          optional, intent(in) :: iverbose
51    logical,          optional, intent(in) :: is_hdf5_file
52    character(len=*), optional, intent(in) :: experiment_name
53
54    type(netcdf_file)                      :: out_file
55    integer                                :: ncol, n_lev_plus1
56    character(5), parameter                :: default_lw_units_str = 'W m-2'
57    character(5)                           :: lw_units_str
58    integer                                :: i_local_verbose
59
60    real(jprb) :: hook_handle
61
62    if (lhook) call dr_hook('radiation_save:save_fluxes',0,hook_handle)
63   
64    if (present(iverbose)) then
65      i_local_verbose = iverbose
66    else
67      i_local_verbose = config%iverbose
68    end if
69
70    ! Work out array dimensions
71    if (config%do_sw) then
72      ncol = size(flux%sw_up,1)
73      n_lev_plus1 = size(flux%sw_up,2)
74    elseif (config%do_lw) then
75      ncol = size(flux%lw_up,1)
76      n_lev_plus1 = size(flux%lw_up,2)
77    else
78      if (i_local_verbose >= 1) then
79        write(nulout,'(a,a,a)') 'Warning: neither longwave nor shortwave computed so ', &
80             &                  file_name,' not written'
81      end if
82      return
83    end if
84
85    if (config%i_gas_model == IGasModelMonochromatic &
86         .and. config%mono_lw_wavelength > 0.0_jprb) then
87      lw_units_str = 'W m-3'
88    else
89      lw_units_str = default_lw_units_str
90    end if
91
92    ! Open the file
93    call out_file%create(trim(file_name), iverbose=i_local_verbose, is_hdf5_file=is_hdf5_file)
94
95    ! Variables stored internally with column varying fastest, but in
96    ! output file column varies most slowly so need to transpose
97    call out_file%transpose_matrices(.true.)
98
99    ! Spectral fluxes in memory are dimensioned (nband,ncol,nlev), but
100    ! are reoriented in the output file to be (nband,nlev,ncol), where
101    ! the convention here is first dimension varying fastest
102    call out_file%permute_3d_arrays( (/ 1, 3, 2 /) )
103
104    ! Define dimensions
105    call out_file%define_dimension("column", ncol)
106    call out_file%define_dimension("half_level", n_lev_plus1)
107
108    if (config%do_save_spectral_flux) then
109      if (config%do_lw) then
110        call out_file%define_dimension("band_lw", config%n_spec_lw)
111      end if
112      if (config%do_sw) then
113        call out_file%define_dimension("band_sw", config%n_spec_sw)
114      end if
115    else if (config%do_surface_sw_spectral_flux) then
116      if (config%do_sw) then
117        call out_file%define_dimension("band_sw", config%n_bands_sw)
118      end if
119    end if
120
121    if (config%do_lw .and. config%do_canopy_fluxes_lw) then
122      call out_file%define_dimension("canopy_band_lw", &
123           &  size(flux%lw_dn_surf_canopy, 1))
124    end if
125    if (config%do_sw .and. config%do_canopy_fluxes_sw) then
126      call out_file%define_dimension("canopy_band_sw", &
127           &  size(flux%sw_dn_diffuse_surf_canopy, 1))
128    end if
129
130    ! Put global attributes
131    call out_file%put_global_attributes( &
132         &   title_str="Radiative flux profiles from the ecRad offline radiation model", &
133         &   references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
134         &   //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
135         &   source_str="ecRad offline radiation model")
136
137    ! Save "experiment" global attribute if present and not empty
138    if (present(experiment_name)) then
139      if (experiment_name /= " ") then
140        call out_file%put_global_attribute("experiment", experiment_name)
141      end if
142    end if
143
144    ! Define variables
145    call out_file%define_variable("pressure_hl", &
146         &   dim2_name="column", dim1_name="half_level", &
147         &   units_str="Pa", long_name="Pressure", &
148         &   standard_name="air_pressure")
149
150    if (config%do_lw) then
151      call out_file%define_variable("flux_up_lw", &
152           &   dim2_name="column", dim1_name="half_level", &
153           &   units_str=lw_units_str, long_name="Upwelling longwave flux", &
154           &   standard_name="upwelling_longwave_flux_in_air")
155      call out_file%define_variable("flux_dn_lw", &
156           &   dim2_name="column", dim1_name="half_level", &
157           &   units_str=lw_units_str, long_name="Downwelling longwave flux", &
158           &   standard_name="downwelling_longwave_flux_in_air")
159      if (config%do_clear) then
160        call out_file%define_variable("flux_up_lw_clear", &
161             &   dim2_name="column", dim1_name="half_level", &
162             &   units_str=lw_units_str, &
163             &   long_name="Upwelling clear-sky longwave flux")
164        call out_file%define_variable("flux_dn_lw_clear", &
165             &   dim2_name="column", dim1_name="half_level", &
166             &   units_str=lw_units_str, &
167             &   long_name="Downwelling clear-sky longwave flux")
168      end if
169
170      if (config%do_lw_derivatives) then
171        call out_file%define_variable("lw_derivative", &
172             &  dim2_name="column", dim1_name="half_level", &
173             &  units_str="1", &
174             &  long_name="Derivative of upwelling LW flux w.r.t. surface value")
175      end if
176
177      if (config%do_save_spectral_flux) then
178        call out_file%define_variable("spectral_flux_up_lw", &
179             &   dim3_name="column", dim2_name="half_level", &
180             &   dim1_name="band_lw", units_str=lw_units_str, &
181             &   long_name="Spectral upwelling longwave flux")
182        call out_file%define_variable("spectral_flux_dn_lw", &
183             &   dim3_name="column", dim2_name="half_level", &
184             &   dim1_name="band_lw", units_str=lw_units_str, &
185             &   long_name="Spectral downwelling longwave flux")
186        if (config%do_clear) then
187          call out_file%define_variable("spectral_flux_up_lw_clear", &
188               &   dim3_name="column", dim2_name="half_level", &
189               &   dim1_name="band_lw", units_str=lw_units_str, &
190               &   long_name="Spectral upwelling clear-sky longwave flux")
191          call out_file%define_variable("spectral_flux_dn_lw_clear", &
192               &   dim3_name="column", dim2_name="half_level", &
193               &   dim1_name="band_lw", units_str=lw_units_str, &
194               &   long_name="Spectral downwelling clear-sky longwave flux")
195        end if
196      end if
197   
198      if (config%do_canopy_fluxes_lw) then
199        call out_file%define_variable("canopy_flux_dn_lw_surf", &
200             &   dim2_name="column", dim1_name="canopy_band_lw", units_str=lw_units_str, &
201             &   long_name="Surface downwelling longwave flux in canopy bands")
202      end if
203
204    end if
205
206    if (config%do_sw) then
207      call out_file%define_variable("flux_up_sw", &
208           &   dim2_name="column", dim1_name="half_level", &
209           &   units_str="W m-2", long_name="Upwelling shortwave flux", &
210           &   standard_name="upwelling_shortwave_flux_in_air")
211      call out_file%define_variable("flux_dn_sw", &
212           &   dim2_name="column", dim1_name="half_level", &
213           &   units_str="W m-2", long_name="Downwelling shortwave flux", &
214           &   standard_name="downwelling_shortwave_flux_in_air")
215      if (config%do_sw_direct) then
216        call out_file%define_variable("flux_dn_direct_sw", &
217             &   dim2_name="column", dim1_name="half_level", &
218             &   units_str="W m-2", &
219             &   long_name="Downwelling direct shortwave flux")
220      end if
221      if (config%do_clear) then
222        call out_file%define_variable("flux_up_sw_clear", &
223             &   dim2_name="column", dim1_name="half_level", &
224             &   units_str="W m-2", &
225             &   long_name="Upwelling clear-sky shortwave flux")
226        call out_file%define_variable("flux_dn_sw_clear", &
227             &   dim2_name="column", dim1_name="half_level", &
228             &   units_str="W m-2", &
229             &   long_name="Downwelling clear-sky shortwave flux")
230        if (config%do_sw_direct) then
231          call out_file%define_variable("flux_dn_direct_sw_clear", &
232               &   dim2_name="column", dim1_name="half_level", &
233               &   units_str="W m-2", &
234               &   long_name="Downwelling clear-sky direct shortwave flux")
235        end if
236      end if
237
238      if (config%do_save_spectral_flux) then
239        call out_file%define_variable("spectral_flux_up_sw", &
240             &   dim3_name="column", dim2_name="half_level", &
241             &   dim1_name="band_sw", units_str="W m-2", &
242             &   long_name="Spectral upwelling shortwave flux")
243        call out_file%define_variable("spectral_flux_dn_sw", &
244             &   dim3_name="column", dim2_name="half_level", &
245             &   dim1_name="band_sw", units_str="W m-2", &
246             &   long_name="Spectral downwelling shortwave flux")
247        if (config%do_sw_direct) then
248          call out_file%define_variable("spectral_flux_dn_direct_sw", &
249               &   dim3_name="column", dim2_name="half_level", &
250               &   dim1_name="band_sw", units_str="W m-2", &
251               &   long_name="Spectral downwelling direct shortwave flux")
252        end if
253        if (config%do_clear) then
254          call out_file%define_variable("spectral_flux_up_sw_clear", &
255               &   dim3_name="column", dim2_name="half_level", &
256               &   dim1_name="band_sw", units_str="W m-2", &
257               &   long_name="Spectral upwelling clear-sky shortwave flux")
258          call out_file%define_variable("spectral_flux_dn_sw_clear", &
259               &   dim3_name="column", dim2_name="half_level", &
260               &   dim1_name="band_sw", units_str="W m-2", &
261               &   long_name="Spectral downwelling clear-sky shortwave flux")
262          if (config%do_sw_direct) then
263            call out_file%define_variable("spectral_flux_dn_direct_sw_clear", &
264                 &   dim3_name="column", dim2_name="half_level", &
265                 &   dim1_name="band_sw", units_str="W m-2", &
266                 &   long_name="Spectral downwelling clear-sky direct shortwave flux")
267          end if
268        end if
269      else if (config%do_surface_sw_spectral_flux) then
270        call out_file%define_variable("spectral_flux_dn_sw_surf", &
271             &   dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
272             &   long_name="Spectral downwelling shortwave flux at surface")
273        call out_file%define_variable("spectral_flux_dn_direct_sw_surf", &
274             &   dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
275             &   long_name="Spectral downwelling direct shortwave flux at surface")
276        if (config%do_clear) then
277          call out_file%define_variable("spectral_flux_dn_sw_surf_clear", &
278               &   dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
279               &   long_name="Spectral downwelling clear-sky shortwave flux at surface")
280          call out_file%define_variable("spectral_flux_dn_direct_sw_surf_clear", &
281               &   dim2_name="column", dim1_name="band_sw", units_str="W m-2", &
282               &   long_name="Spectral downwelling clear-sky direct shortwave flux at surface")
283        end if
284      end if
285   
286      if (config%do_canopy_fluxes_sw) then
287        call out_file%define_variable("canopy_flux_dn_diffuse_sw_surf", &
288             &   dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
289             &   long_name="Surface downwelling diffuse shortwave flux in canopy bands")
290        call out_file%define_variable("canopy_flux_dn_direct_sw_surf", &
291             &   dim2_name="column", dim1_name="canopy_band_sw", units_str="W m-2", &
292             &   long_name="Surface downwelling direct shortwave flux in canopy bands")
293      end if
294
295    end if
296   
297    if (config%do_lw .and. config%do_clouds) then
298      call out_file%define_variable("cloud_cover_lw", &
299           &  dim1_name="column", units_str="1", &
300           &  long_name="Total cloud cover diagnosed by longwave solver", &
301           &  standard_name="cloud_area_fraction")
302    end if
303    if (config%do_sw .and. config%do_clouds) then
304      call out_file%define_variable("cloud_cover_sw", &
305           &  dim1_name="column", units_str="1", &
306           &  long_name="Total cloud cover diagnosed by shortwave solver", &
307           &  standard_name="cloud_area_fraction")
308    end if
309
310    ! Write variables
311
312    call out_file%put("pressure_hl", thermodynamics%pressure_hl)
313
314    if (config%do_lw) then
315      call out_file%put("flux_up_lw", flux%lw_up)
316      call out_file%put("flux_dn_lw", flux%lw_dn)
317      if (config%do_clear) then
318        call out_file%put("flux_up_lw_clear", flux%lw_up_clear)
319        call out_file%put("flux_dn_lw_clear", flux%lw_dn_clear)
320      end if
321
322      if (config%do_lw_derivatives) then
323        call out_file%put("lw_derivative", flux%lw_derivatives)
324      end if
325
326      if (config%do_save_spectral_flux) then
327        call out_file%put("spectral_flux_up_lw", flux%lw_up_band)
328        call out_file%put("spectral_flux_dn_lw", flux%lw_dn_band)
329        if (config%do_clear) then
330          call out_file%put("spectral_flux_up_lw_clear", flux%lw_up_clear_band)
331          call out_file%put("spectral_flux_dn_lw_clear", flux%lw_dn_clear_band)
332        end if
333      end if
334
335      if (config%do_canopy_fluxes_lw) then
336        call out_file%put("canopy_flux_dn_lw_surf", flux%lw_dn_surf_canopy, &
337             &            do_transp = .false.)
338      end if
339
340    end if
341
342    if (config%do_sw) then
343      call out_file%put("flux_up_sw", flux%sw_up)
344      call out_file%put("flux_dn_sw", flux%sw_dn)
345      if (config%do_sw_direct) then
346        call out_file%put("flux_dn_direct_sw", flux%sw_dn_direct)
347      end if
348      if (config%do_clear) then
349        call out_file%put("flux_up_sw_clear", flux%sw_up_clear)
350        call out_file%put("flux_dn_sw_clear", flux%sw_dn_clear)
351        if (config%do_sw_direct) then
352          call out_file%put("flux_dn_direct_sw_clear", flux%sw_dn_direct_clear)
353        end if
354      end if
355
356      if (config%do_save_spectral_flux) then
357        call out_file%put("spectral_flux_up_sw", flux%sw_up_band)
358        call out_file%put("spectral_flux_dn_sw", flux%sw_dn_band)
359        if (config%do_sw_direct) then
360          call out_file%put("spectral_flux_dn_direct_sw", &
361               &   flux%sw_dn_direct_band)
362        end if
363        if (config%do_clear) then
364          call out_file%put("spectral_flux_up_sw_clear", flux%sw_up_clear_band)
365          call out_file%put("spectral_flux_dn_sw_clear", flux%sw_dn_clear_band)
366          if (config%do_sw_direct) then
367            call out_file%put("spectral_flux_dn_direct_sw_clear", &
368                 &   flux%sw_dn_direct_clear_band)
369          end if
370        end if
371      else if (config%do_surface_sw_spectral_flux) then
372        call out_file%put("spectral_flux_dn_sw_surf", flux%sw_dn_surf_band, &
373               &   do_transp=.false.)
374        call out_file%put("spectral_flux_dn_direct_sw_surf", flux%sw_dn_direct_surf_band, &
375               &   do_transp=.false.)
376        if (config%do_clear) then
377          call out_file%put("spectral_flux_dn_sw_surf_clear", flux%sw_dn_surf_clear_band, &
378               &   do_transp=.false.)
379          call out_file%put("spectral_flux_dn_direct_sw_surf_clear", &
380               &            flux%sw_dn_direct_surf_clear_band, do_transp=.false.)
381        end if
382      end if
383
384      if (config%do_canopy_fluxes_sw) then
385        call out_file%put("canopy_flux_dn_diffuse_sw_surf", flux%sw_dn_diffuse_surf_canopy, &
386             &            do_transp = .false.)
387        call out_file%put("canopy_flux_dn_direct_sw_surf",  flux%sw_dn_direct_surf_canopy, &
388             &            do_transp = .false.)
389      end if
390
391    end if
392
393    if (config%do_lw .and. config%do_clouds) then
394      call out_file%put("cloud_cover_lw", flux%cloud_cover_lw)
395    end if
396    if (config%do_sw .and. config%do_clouds) then
397      call out_file%put("cloud_cover_sw", flux%cloud_cover_sw)
398    end if
399
400    ! Close file
401    call out_file%close()
402
403    if (lhook) call dr_hook('radiation_save:save_fluxes',1,hook_handle)
404
405  end subroutine save_fluxes
406 
407
408  !---------------------------------------------------------------------
409  ! Save intermediate radiative properties, specifically the
410  ! scattering and absorption properties at each g-point/band
411  subroutine save_radiative_properties(file_name, nlev, &
412       &  istartcol, iendcol, &
413       &  config, single_level, thermodynamics, cloud, &
414       &  planck_hl, lw_emission, lw_albedo, &
415       &  sw_albedo_direct, sw_albedo_diffuse, &
416       &  incoming_sw, &
417       &  od_lw, ssa_lw, g_lw, &
418       &  od_sw, ssa_sw, g_sw, &
419       &  od_lw_cloud, ssa_lw_cloud, g_lw_cloud, &
420       &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud)
421
422    use radiation_config,        only : config_type
423    use radiation_single_level,  only : single_level_type
424    use radiation_thermodynamics,only : thermodynamics_type
425    use radiation_cloud,         only : cloud_type
426    use easy_netcdf
427
428    character(len=*),         intent(in) :: file_name
429    type(config_type),        intent(in) :: config
430    type(single_level_type),  intent(in) :: single_level
431    type(thermodynamics_type),intent(in) :: thermodynamics
432    type(cloud_type),         intent(in) :: cloud
433
434    integer, intent(in) :: nlev, istartcol, iendcol
435
436    ! Input variables, as defined in radiation_interface.F90
437
438    ! Layer optical depth, single scattering albedo and asymmetry factor of
439    ! gases and aerosols at each shortwave g-point
440    real(jprb), intent(in), dimension(config%n_g_sw,nlev,istartcol:iendcol) :: od_sw, ssa_sw, g_sw
441
442   ! Layer optical depth, single scattering albedo and asymmetry factor of
443    ! hydrometeors in each shortwave band
444    real(jprb), intent(in), dimension(config%n_bands_sw,nlev,istartcol:iendcol)   :: &
445         &  od_sw_cloud, ssa_sw_cloud, g_sw_cloud
446
447    ! Direct and diffuse surface albedo, and the incoming shortwave
448    ! flux into a plane perpendicular to the incoming radiation at
449    ! top-of-atmosphere in each of the shortwave g-points
450    real(jprb), intent(in), dimension(config%n_g_sw,istartcol:iendcol) &
451         &  :: sw_albedo_direct, sw_albedo_diffuse, incoming_sw
452
453    ! Layer optical depth, single scattering albedo and asymmetry factor of
454    ! gases and aerosols at each longwave g-point, where the latter
455    ! two variables are only defined if aerosol longwave scattering is
456    ! enabled (otherwise both are treated as zero).
457    real(jprb), intent(in), dimension(config%n_g_lw,nlev,istartcol:iendcol) :: od_lw
458    real(jprb), intent(in), dimension(config%n_g_lw_if_scattering,nlev,istartcol:iendcol) :: &
459         &  ssa_lw, g_lw
460
461    ! Layer optical depth, single scattering albedo and asymmetry factor of
462    ! hydrometeors in each longwave band, where the latter two
463    ! variables are only defined if hydrometeor longwave scattering is
464    ! enabled (otherwise both are treated as zero).
465    real(jprb), intent(in), dimension(config%n_bands_lw,nlev,istartcol:iendcol) :: od_lw_cloud
466    real(jprb), intent(in), dimension(config%n_bands_lw_if_scattering,nlev,istartcol:iendcol) :: &
467         &  ssa_lw_cloud, g_lw_cloud
468
469    ! The Planck function (emitted flux from a black body) at half
470    ! levels and at the surface at each longwave g-point
471    real(jprb), intent(in), dimension(config%n_g_lw,nlev+1,istartcol:iendcol) :: planck_hl
472
473    ! Emission (Planck*emissivity) and albedo (1-emissivity) at the
474    ! surface at each longwave g-point
475    real(jprb), dimension(config%n_g_lw, istartcol:iendcol) :: lw_emission, lw_albedo
476
477    integer :: n_col_local ! Number of columns from istartcol to iendcol
478
479    ! Object for output NetCDF file
480    type(netcdf_file) :: out_file
481
482    n_col_local = iendcol + 1 - istartcol
483
484    ! Alas the NetCDF library is not thread-safe for writing, so we
485    ! must write radiative-property files serially
486
487    !$OMP CRITICAL
488
489    ! Open the file
490    call out_file%create(trim(file_name), iverbose=config%iverbose)
491
492    ! Configure matrix and 3D-array orientation
493    call out_file%transpose_matrices(.true.)
494
495    ! Sometimes the Planck function values are very large or small
496    ! even if the fluxes are within a manageable range
497    call out_file%double_precision(.true.)
498
499    ! Define dimensions
500    !    call out_file%define_dimension("column", n_col_local)
501    call out_file%define_dimension("column", 0) ! "Unlimited" dimension
502    call out_file%define_dimension("level", nlev)
503    call out_file%define_dimension("half_level", nlev+1)
504    if (config%do_clouds) then
505      call out_file%define_dimension("level_interface", nlev-1)
506    end if
507
508    if (config%do_lw) then
509      call out_file%define_dimension("gpoint_lw", config%n_g_lw)
510      if (config%do_clouds) then
511        call out_file%define_dimension("band_lw", config%n_bands_lw)
512      end if
513    end if
514    if (config%do_sw) then
515      call out_file%define_dimension("gpoint_sw", config%n_g_sw)
516      if (config%do_clouds) then
517        call out_file%define_dimension("band_sw", config%n_bands_sw)
518      end if
519    end if
520
521    ! Put global attributes
522    call out_file%put_global_attributes( &
523         &   title_str="Spectral radiative properties from the ecRad offline radiation model", &
524         &   references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
525         &   //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
526         &   source_str="ecRad offline radiation model")
527
528    ! Define variables
529    call out_file%define_variable("pressure_hl", &
530         &  dim2_name="column", dim1_name="half_level", &
531         &  units_str="Pa", long_name="Pressure on half-levels")
532
533    if (allocated(thermodynamics%h2o_sat_liq) .and. config%use_aerosols) then
534      call out_file%define_variable("q_sat_liquid", &
535           &  dim2_name="column", dim1_name="level", &
536           &  units_str="kg kg-1", long_name="Specific humidity at liquid saturation")
537    end if
538
539    if (config%do_sw) then
540      call out_file%define_variable("cos_solar_zenith_angle", &
541           &  dim1_name="column", units_str="1", &
542           &  long_name="Cosine of the solar zenith angle")
543    end if
544
545    if (config%do_clouds) then
546      call out_file%define_variable("cloud_fraction", &
547           &  dim2_name="column", dim1_name="level", &
548           &  units_str="1", long_name="Cloud fraction")
549      call out_file%define_variable("overlap_param", &
550           &  dim2_name="column", dim1_name="level_interface", &
551           &  units_str="1", long_name="Cloud overlap parameter")
552    end if
553
554    if (config%do_lw) then
555      call out_file%define_variable("planck_hl", &
556           &  dim3_name="column", dim2_name="half_level", dim1_name="gpoint_lw", &
557           &  units_str="W m-2", long_name="Planck function on half-levels")
558      call out_file%define_variable("lw_emission", &
559           &  dim2_name="column", dim1_name="gpoint_lw", &
560           &  units_str="W m-2", long_name="Longwave surface emission")
561      call out_file%define_variable("lw_emissivity", &
562           &  dim2_name="column", dim1_name="gpoint_lw", &
563           &  units_str="1", long_name="Surface longwave emissivity")
564
565      call out_file%define_variable("od_lw", &
566           &  dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
567           &  units_str="1", long_name="Clear-sky longwave optical depth")
568      if (config%do_lw_aerosol_scattering) then
569        call out_file%define_variable("ssa_lw", &
570           &  dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
571           &  units_str="1", long_name="Clear-sky longwave single scattering albedo")
572        call out_file%define_variable("asymmetry_lw", &
573           &  dim3_name="column", dim2_name="level", dim1_name="gpoint_lw", &
574           &  units_str="1", long_name="Clear-sky longwave asymmetry factor")
575      end if
576
577      if (config%do_clouds) then
578        call out_file%define_variable("od_lw_cloud", &
579             &  dim3_name="column", dim2_name="level", dim1_name="band_lw", &
580             &  units_str="1", long_name="In-cloud longwave optical depth")
581        if (config%do_lw_cloud_scattering) then
582          call out_file%define_variable("ssa_lw_cloud", &
583               &  dim3_name="column", dim2_name="level", dim1_name="band_lw", &
584               &  units_str="1", long_name="Cloud longwave single scattering albedo")
585          call out_file%define_variable("asymmetry_lw_cloud", &
586               &  dim3_name="column", dim2_name="level", dim1_name="band_lw", &
587               &  units_str="1", long_name="Cloud longwave asymmetry factor")
588        end if
589      end if ! do_clouds
590    end if ! do_lw
591   
592    if (config%do_sw) then
593      call out_file%define_variable("incoming_sw", &
594           &  dim2_name="column", dim1_name="gpoint_sw", &
595           &  units_str="W m-2", long_name="Incoming shortwave flux at top-of-atmosphere in direction of sun")
596
597      call out_file%define_variable("sw_albedo", &
598           &  dim2_name="column", dim1_name="gpoint_sw", &
599           &  units_str="1", long_name="Surface shortwave albedo to diffuse radiation")
600      call out_file%define_variable("sw_albedo_direct", &
601           &  dim2_name="column", dim1_name="gpoint_sw", &
602           &  units_str="1", long_name="Surface shortwave albedo to direct radiation")
603
604      call out_file%define_variable("od_sw", &
605           &  dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
606           &  units_str="1", long_name="Clear-sky shortwave optical depth")
607      call out_file%define_variable("ssa_sw", &
608           &  dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
609           &  units_str="1", long_name="Clear-sky shortwave single scattering albedo")
610      call out_file%define_variable("asymmetry_sw", &
611           &  dim3_name="column", dim2_name="level", dim1_name="gpoint_sw", &
612           &  units_str="1", long_name="Clear-sky shortwave asymmetry factor")
613
614      if (config%do_clouds) then
615        call out_file%define_variable("od_sw_cloud", &
616             &  dim3_name="column", dim2_name="level", dim1_name="band_sw", &
617             &  units_str="1", long_name="In-cloud shortwave optical depth")
618        call out_file%define_variable("ssa_sw_cloud", &
619             &  dim3_name="column", dim2_name="level", dim1_name="band_sw", &
620             &  units_str="1", long_name="Cloud shortwave single scattering albedo")
621        call out_file%define_variable("asymmetry_sw_cloud", &
622             &  dim3_name="column", dim2_name="level", dim1_name="band_sw", &
623             &  units_str="1", long_name="Cloud shortwave asymmetry factor")
624      end if
625    end if
626   
627    if (config%do_clouds) then
628      if (allocated(cloud%fractional_std)) then
629        call out_file%define_variable("fractional_std", &
630             &  dim2_name="column", dim1_name="level", units_str="1", &
631             &  long_name="Fractional standard deviation of cloud optical depth")
632      end if
633      if (allocated(cloud%inv_cloud_effective_size)) then
634        call out_file%define_variable("inv_cloud_effective_size", &
635             &  dim2_name="column", dim1_name="level", units_str="m-1", &
636             &  long_name="Inverse of cloud effective horizontal size")
637      end if
638      if (allocated(cloud%inv_inhom_effective_size)) then
639        call out_file%define_variable("inv_inhom_effective_size", &
640             &  dim2_name="column", dim1_name="level", units_str="m-1", &
641             &  long_name="Inverse of cloud inhomogeneity effective horizontal size")
642      end if
643   end if
644
645    ! Write variables
646    call out_file%put("pressure_hl", thermodynamics%pressure_hl(istartcol:iendcol,:))
647
648    if (allocated(thermodynamics%h2o_sat_liq) .and. config%use_aerosols) then
649      call out_file%put("q_sat_liquid", thermodynamics%h2o_sat_liq(istartcol:iendcol,:))
650    end if
651
652    if (config%do_clouds) then
653      call out_file%put("cloud_fraction", cloud%fraction(istartcol:iendcol,:))
654      call out_file%put("overlap_param", cloud%overlap_param(istartcol:iendcol,:))
655    end if
656
657    if (config%do_sw) then
658      call out_file%put("cos_solar_zenith_angle", single_level%cos_sza(istartcol:iendcol))
659      call out_file%put("sw_albedo", sw_albedo_diffuse, do_transp=.false.)
660      call out_file%put("sw_albedo_direct", sw_albedo_direct, do_transp=.false.)
661    end if
662
663    if (config%do_lw) then
664      call out_file%put("lw_emissivity", 1.0_jprb - lw_albedo, do_transp=.false.)
665      call out_file%put("planck_hl", planck_hl)
666      call out_file%put("lw_emission", lw_emission, do_transp=.false.)
667     
668      call out_file%put("od_lw", od_lw)
669      if (config%do_lw_aerosol_scattering) then
670        call out_file%put("ssa_lw", ssa_lw)
671        call out_file%put("asymmetry_lw", g_lw)
672      end if
673     
674      if (config%do_clouds) then
675        call out_file%put("od_lw_cloud", od_lw_cloud)
676        if (config%do_lw_cloud_scattering) then
677          call out_file%put("ssa_lw_cloud", ssa_lw_cloud)
678          call out_file%put("asymmetry_lw_cloud", g_lw_cloud)
679        end if
680      end if
681    end if
682   
683    if (config%do_sw) then
684      call out_file%put("incoming_sw", incoming_sw, do_transp=.false.)
685     
686      call out_file%put("od_sw", od_sw)
687      call out_file%put("ssa_sw", ssa_sw)
688      call out_file%put("asymmetry_sw", g_sw)
689     
690      if (config%do_clouds) then
691        call out_file%put("od_sw_cloud", od_sw_cloud)
692        call out_file%put("ssa_sw_cloud", ssa_sw_cloud)
693        call out_file%put("asymmetry_sw_cloud", g_sw_cloud)
694      end if
695    end if
696   
697    if (config%do_clouds) then
698      if (allocated(cloud%fractional_std)) then
699        call out_file%put("fractional_std", cloud%fractional_std(istartcol:iendcol,:))
700      end if
701      if (allocated(cloud%inv_cloud_effective_size)) then
702        call out_file%put("inv_cloud_effective_size", cloud%inv_cloud_effective_size(istartcol:iendcol,:))
703      end if
704      if (allocated(cloud%inv_inhom_effective_size)) then
705        call out_file%put("inv_inhom_effective_size", cloud%inv_inhom_effective_size(istartcol:iendcol,:))
706      end if
707    end if
708
709    ! Close the file
710    call out_file%close()
711
712    !$OMP END CRITICAL
713   
714  end subroutine save_radiative_properties
715 
716
717  !---------------------------------------------------------------------
718  ! Save inputs to the radiation scheme
719  subroutine save_inputs(file_name, config, single_level, thermodynamics, &
720       &                 gas, cloud, aerosol, lat, lon, iverbose)
721    use yomhook,                  only : lhook, dr_hook
722
723    use radiation_config,         only : config_type
724    use radiation_single_level,   only : single_level_type
725    use radiation_thermodynamics, only : thermodynamics_type
726    use radiation_gas
727    use radiation_cloud,          only : cloud_type
728    use radiation_aerosol,        only : aerosol_type
729    use easy_netcdf
730
731    character(len=*),             intent(in)   :: file_name
732    type(config_type),            intent(in)   :: config
733    type(single_level_type),      intent(in)   :: single_level
734    type(thermodynamics_type),    intent(in)   :: thermodynamics
735    type(gas_type),               intent(inout):: gas
736    type(cloud_type),             intent(in)   :: cloud
737    type(aerosol_type), optional, intent(in)   :: aerosol
738    real(jprb),         optional, intent(in)   :: lat(:), lon(:)
739    integer,            optional, intent(in)   :: iverbose
740
741    real(jprb), allocatable :: mixing_ratio(:,:)
742    real(jprb), allocatable :: aerosol_mmr(:,:,:)
743    real(jprb), allocatable :: seed(:)
744    integer       :: i_local_verbose
745    integer       :: ncol, nlev
746    integer       :: jgas
747    character(32) :: var_name, long_name
748
749    ! Object for output NetCDF file
750    type(netcdf_file) :: out_file
751
752    logical :: do_aerosol
753
754    real(jprb) :: hook_handle
755
756    if (lhook) call dr_hook('radiation_save:save_inputs',0,hook_handle)
757
758    if (present(iverbose)) then
759      i_local_verbose = iverbose
760    else
761      i_local_verbose = config%iverbose
762    end if
763
764    ! Work out array dimensions
765    ncol = size(thermodynamics%pressure_hl,1)
766    nlev = size(thermodynamics%pressure_hl,2)
767    nlev = nlev - 1
768   
769    do_aerosol = config%use_aerosols .and. present(aerosol)
770
771    ! Open the file
772    call out_file%create(trim(file_name), iverbose=i_local_verbose)
773
774    ! Variables stored internally with column varying fastest, but in
775    ! output file column varies most slowly so need to transpose
776    call out_file%transpose_matrices(.true.)
777
778    ! In the case of aerosols we convert dimensions (ncol,nlev,ntype)
779    ! in memory to (nlev,ntype,ncol) in file (in both cases the first
780    ! dimension varying fastest).
781    call out_file%permute_3d_arrays( (/ 2, 3, 1 /) ) ! For aerosols
782
783    ! Define dimensions
784    call out_file%define_dimension("column",     ncol)
785    call out_file%define_dimension("level",      nlev)
786    call out_file%define_dimension("half_level", nlev+1)
787    if (allocated(cloud%overlap_param)) then
788      call out_file%define_dimension("level_interface", nlev-1)
789    end if
790    call out_file%define_dimension("sw_albedo_band", &
791         &                         size(single_level%sw_albedo,2))
792    call out_file%define_dimension("lw_emissivity_band", &
793         &                         size(single_level%lw_emissivity,2))
794
795    if (do_aerosol) then
796      call out_file%define_dimension("aerosol_type", size(aerosol%mixing_ratio,3))
797    end if
798
799    ! Put global attributes
800    call out_file%put_global_attributes( &
801         &   title_str="Input profiles to the ecRad offline radiation model", &
802         &   references_str="Hogan, R. J., and A. Bozzo, 2018: A flexible and efficient radiation " &
803         &   //"scheme for the ECMWF model. J. Adv. Modeling Earth Sys., 10, 1990–2008", &
804         &   source_str="ecRad offline radiation model")
805
806    ! Define single-level variables
807    call out_file%define_variable("solar_irradiance", &
808         &   units_str="W m-2", long_name="Solar irradiance at Earth's orbit")
809    if (present(lat)) then
810      call out_file%define_variable("lat", &
811           &   dim1_name="column", units_str="degrees_north", long_name="Latitude")
812    end if
813    if (present(lon)) then
814      call out_file%define_variable("lon", &
815           &   dim1_name="column", units_str="degrees_east", long_name="Longitude")
816    end if
817    call out_file%define_variable("skin_temperature", &
818         &   dim1_name="column", units_str="K", long_name="Skin_temperature")
819    if (config%do_sw) then
820      call out_file%define_variable("cos_solar_zenith_angle", &
821           &   dim1_name="column", units_str="1", &
822           &   long_name="Cosine of the solar zenith angle")
823    end if
824
825    if (allocated(single_level%sw_albedo_direct)) then
826      call out_file%define_variable("sw_albedo", &
827           &   dim2_name="column", dim1_name="sw_albedo_band", &
828           &   units_str="1", long_name="Shortwave surface albedo to diffuse radiation")
829            call out_file%define_variable("sw_albedo_direct", &
830           &   dim2_name="column", dim1_name="sw_albedo_band", &
831           &   units_str="1", long_name="Shortwave surface albedo to direct radiation")
832    else
833      call out_file%define_variable("sw_albedo", &
834           &   dim2_name="column", dim1_name="sw_albedo_band", &
835           &   units_str="1", long_name="Shortwave surface albedo")
836     
837    end if
838    call out_file%define_variable("lw_emissivity", &
839         &   dim2_name="column", dim1_name="lw_emissivity_band", &
840         &   units_str="1", long_name="Longwave surface emissivity")
841
842    if (allocated(single_level%iseed)) then
843      call out_file%define_variable("iseed", &
844           &   dim1_name="column", units_str="1", is_double=.true., &
845           &   long_name="Seed for random-number generator")
846    end if
847
848    ! Define thermodynamic variables on half levels
849    call out_file%define_variable("pressure_hl", &
850         &   dim2_name="column", dim1_name="half_level", &
851         &   units_str="Pa", long_name="Pressure")
852    call out_file%define_variable("temperature_hl", &
853         &   dim2_name="column", dim1_name="half_level", &
854         &   units_str="K", long_name="Temperature")
855
856    ! Define gas mixing ratios on full levels
857    call out_file%define_variable("q", &
858         &   dim2_name="column", dim1_name="level", &
859         &   units_str="1", long_name="Specific humidity")
860    call out_file%define_variable("o3_mmr", &
861         &   dim2_name="column", dim1_name="level", &
862         &   units_str="1", long_name="Ozone mass mixing ratio")
863    do jgas = 1,NMaxGases
864      if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
865        write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
866        write(long_name,'(a,a)') trim(GasName(jgas)), ' volume mixing ratio'
867        call out_file%define_variable(trim(var_name), &
868             &   dim2_name="column", dim1_name="level", &
869             &   units_str="1", long_name=trim(long_name))
870      end if
871    end do
872
873    if (config%do_clouds) then
874      ! Define cloud variables on full levels
875      call out_file%define_variable("cloud_fraction", &
876           &   dim2_name="column", dim1_name="level", &
877           &   units_str="1", long_name="Cloud fraction")
878      call out_file%define_variable("q_liquid", &
879           &   dim2_name="column", dim1_name="level", &
880           &   units_str="1", long_name="Gridbox-mean liquid water mixing ratio")
881      call out_file%define_variable("q_ice", &
882           &   dim2_name="column", dim1_name="level", &
883           &   units_str="1", long_name="Gridbox-mean ice water mixing ratio")
884      call out_file%define_variable("re_liquid", &
885           &   dim2_name="column", dim1_name="level", &
886           &   units_str="m", long_name="Ice effective radius")
887      if (allocated(cloud%re_ice)) then
888        call out_file%define_variable("re_ice", &
889             &   dim2_name="column", dim1_name="level", &
890             &   units_str="m", long_name="Ice effective radius")
891      end if
892      if (allocated(cloud%overlap_param)) then
893        call out_file%define_variable("overlap_param", &
894             &  dim2_name="column", dim1_name="level_interface", &
895             &  units_str="1", long_name="Cloud overlap parameter")
896      end if
897      if (allocated(cloud%fractional_std)) then
898        call out_file%define_variable("fractional_std", &
899             &  dim2_name="column", dim1_name="level", units_str="1", &
900             &  long_name="Fractional standard deviation of cloud optical depth")
901      end if
902      if (allocated(cloud%inv_cloud_effective_size)) then
903        call out_file%define_variable("inv_cloud_effective_size", &
904             &   dim2_name="column", dim1_name="level", units_str="m-1", &
905             &   long_name="Inverse of cloud effective horizontal size")
906      end if
907      if (allocated(cloud%inv_inhom_effective_size)) then
908        call out_file%define_variable("inv_inhom_effective_size", &
909             &  dim2_name="column", dim1_name="level", units_str="m-1", &
910             &  long_name="Inverse of cloud inhomogeneity effective horizontal size")
911      end if
912    end if ! do_clouds
913
914    ! Define aerosol mass mixing ratio
915    if (do_aerosol) then
916      call out_file%define_variable("aerosol_mmr", &
917           &   dim3_name="column", dim2_name="aerosol_type", &
918           &   dim1_name="level", units_str="kg kg-1", &
919           &   long_name="Aerosol mass mixing ratio")
920    end if
921
922    ! Write variables
923    call out_file%put("solar_irradiance", single_level%solar_irradiance)
924    if (present(lat)) then
925      call out_file%put("lat", lat)
926    end if
927    if (present(lon)) then
928      call out_file%put("lon", lon)
929    end if
930    call out_file%put("skin_temperature", single_level%skin_temperature)
931    if (config%do_sw) then
932      call out_file%put("cos_solar_zenith_angle", single_level%cos_sza)
933    end if
934    call out_file%put("sw_albedo", single_level%sw_albedo)
935    if (allocated(single_level%sw_albedo_direct)) then
936      call out_file%put("sw_albedo_direct", single_level%sw_albedo_direct)
937    end if
938    call out_file%put("lw_emissivity", single_level%lw_emissivity)
939    if (config%do_clouds .and. allocated(single_level%iseed)) then
940      allocate(seed(ncol))
941      seed = single_level%iseed
942      call out_file%put("iseed", seed)
943      deallocate(seed)
944    end if
945
946    call out_file%put("pressure_hl", thermodynamics%pressure_hl)
947    call out_file%put("temperature_hl", thermodynamics%temperature_hl)
948
949    allocate(mixing_ratio(ncol,nlev))
950    call gas%get(IH2O, IMassMixingRatio, mixing_ratio)
951    call out_file%put("q", mixing_ratio)
952    call gas%get(IO3, IMassMixingRatio, mixing_ratio)
953    call out_file%put("o3_mmr", mixing_ratio)
954    do jgas = 1,NMaxGases
955      if (gas%is_present(jgas) .and. jgas /= IH2O .and. jgas /= IO3) then
956        write(var_name,'(a,a)') trim(GasLowerCaseName(jgas)), '_vmr'
957        call gas%get(jgas, IVolumeMixingRatio, mixing_ratio)
958        call out_file%put(trim(var_name), mixing_ratio)
959      end if
960    end do
961    deallocate(mixing_ratio)
962
963    if (config%do_clouds) then
964      call out_file%put("cloud_fraction", cloud%fraction)
965      call out_file%put("q_liquid", cloud%q_liq)
966      call out_file%put("q_ice", cloud%q_ice)
967      call out_file%put("re_liquid", cloud%re_liq)
968      if (allocated(cloud%re_ice)) then
969        call out_file%put("re_ice", cloud%re_ice)
970      end if
971      if (allocated(cloud%overlap_param)) then
972        call out_file%put("overlap_param", cloud%overlap_param)
973      end if
974      if (allocated(cloud%fractional_std)) then
975        call out_file%put("fractional_std", cloud%fractional_std)
976      end if
977      if (allocated(cloud%inv_cloud_effective_size)) then
978        call out_file%put("inv_cloud_effective_size", cloud%inv_cloud_effective_size)
979      end if
980      if (allocated(cloud%inv_inhom_effective_size)) then
981        call out_file%put("inv_inhom_effective_size", cloud%inv_inhom_effective_size)
982      end if
983    end if
984
985    if (do_aerosol) then
986      allocate(aerosol_mmr(ncol, nlev, size(aerosol%mixing_ratio,3)))
987      aerosol_mmr = 0.0_jprb
988      aerosol_mmr(:,aerosol%istartlev:aerosol%iendlev,:) = aerosol%mixing_ratio
989      call out_file%put("aerosol_mmr", aerosol_mmr)
990      deallocate(aerosol_mmr)
991    end if
992
993    ! Close the file
994    call out_file%close()
995
996    if (lhook) call dr_hook('radiation_save:save_inputs',1,hook_handle)
997   
998  end subroutine save_inputs
999
1000end module radiation_save
Note: See TracBrowser for help on using the repository browser.