| 1 | ! ecrad_driver.F90 - Driver for offline ECRAD 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 | ! ECRAD is the radiation scheme used in the ECMWF Integrated |
|---|
| 16 | ! Forecasting System in cycle 43R3 and later. Several solvers are |
|---|
| 17 | ! available, including McICA, Tripleclouds and SPARTACUS (the Speedy |
|---|
| 18 | ! Algorithm for Radiative Transfer through Cloud Sides, a modification |
|---|
| 19 | ! of the two-stream formulation of shortwave and longwave radiative |
|---|
| 20 | ! transfer to account for 3D radiative effects). Gas optical |
|---|
| 21 | ! properties are provided by the RRTM-G gas optics scheme. |
|---|
| 22 | |
|---|
| 23 | ! This program takes three arguments: |
|---|
| 24 | ! 1) Namelist file to configure the radiation calculation |
|---|
| 25 | ! 2) Name of a NetCDF file containing one or more atmospheric profiles |
|---|
| 26 | ! 3) Name of output NetCDF file |
|---|
| 27 | |
|---|
| 28 | program ecrad_driver |
|---|
| 29 | |
|---|
| 30 | ! -------------------------------------------------------- |
|---|
| 31 | ! Section 1: Declarations |
|---|
| 32 | ! -------------------------------------------------------- |
|---|
| 33 | use parkind1, only : jprb, jprd ! Working/double precision |
|---|
| 34 | |
|---|
| 35 | use radiation_io, only : nulout |
|---|
| 36 | use radiation_interface, only : setup_radiation, radiation, set_gas_units |
|---|
| 37 | use radiation_config, only : config_type |
|---|
| 38 | use radiation_single_level, only : single_level_type |
|---|
| 39 | use radiation_thermodynamics, only : thermodynamics_type |
|---|
| 40 | use radiation_gas, only : gas_type, & |
|---|
| 41 | & IVolumeMixingRatio, IMassMixingRatio, & |
|---|
| 42 | & IH2O, ICO2, IO3, IN2O, ICO, ICH4, IO2, ICFC11, ICFC12, & |
|---|
| 43 | & IHCFC22, ICCl4, GasName, GasLowerCaseName, NMaxGases |
|---|
| 44 | use radiation_cloud, only : cloud_type |
|---|
| 45 | use radiation_aerosol, only : aerosol_type |
|---|
| 46 | use radiation_flux, only : flux_type |
|---|
| 47 | use radiation_save, only : save_fluxes, save_net_fluxes, & |
|---|
| 48 | & save_inputs, save_sw_diagnostics |
|---|
| 49 | use radiation_general_cloud_optics, only : save_general_cloud_optics |
|---|
| 50 | use ecrad_driver_config, only : driver_config_type |
|---|
| 51 | use ecrad_driver_read_input, only : read_input |
|---|
| 52 | use easy_netcdf |
|---|
| 53 | use print_matrix_mod, only : print_matrix |
|---|
| 54 | |
|---|
| 55 | implicit none |
|---|
| 56 | |
|---|
| 57 | ! The NetCDF file containing the input profiles |
|---|
| 58 | type(netcdf_file) :: file |
|---|
| 59 | |
|---|
| 60 | ! Derived types for the inputs to the radiation scheme |
|---|
| 61 | type(config_type) :: config |
|---|
| 62 | type(single_level_type) :: single_level |
|---|
| 63 | type(thermodynamics_type) :: thermodynamics |
|---|
| 64 | type(gas_type) :: gas |
|---|
| 65 | type(cloud_type) :: cloud |
|---|
| 66 | type(aerosol_type) :: aerosol |
|---|
| 67 | |
|---|
| 68 | ! Configuration specific to this driver |
|---|
| 69 | type(driver_config_type) :: driver_config |
|---|
| 70 | |
|---|
| 71 | ! Derived type containing outputs from the radiation scheme |
|---|
| 72 | type(flux_type) :: flux |
|---|
| 73 | |
|---|
| 74 | integer :: ncol, nlev ! Number of columns and levels |
|---|
| 75 | integer :: istartcol, iendcol ! Range of columns to process |
|---|
| 76 | |
|---|
| 77 | ! Name of file names specified on command line |
|---|
| 78 | character(len=512) :: file_name |
|---|
| 79 | integer :: istatus ! Result of command_argument_count |
|---|
| 80 | |
|---|
| 81 | ! For parallel processing of multiple blocks |
|---|
| 82 | integer :: jblock, nblock ! Block loop index and number |
|---|
| 83 | |
|---|
| 84 | ! Mapping matrix for shortwave spectral diagnostics |
|---|
| 85 | real(jprb), allocatable :: sw_diag_mapping(:,:) |
|---|
| 86 | |
|---|
| 87 | #ifndef NO_OPENMP |
|---|
| 88 | ! OpenMP functions |
|---|
| 89 | integer, external :: omp_get_thread_num |
|---|
| 90 | real(kind=jprd), external :: omp_get_wtime |
|---|
| 91 | ! Start/stop time in seconds |
|---|
| 92 | real(kind=jprd) :: tstart, tstop |
|---|
| 93 | #endif |
|---|
| 94 | |
|---|
| 95 | ! For demonstration of get_sw_weights later on |
|---|
| 96 | ! Ultraviolet weightings |
|---|
| 97 | !integer :: nweight_uv |
|---|
| 98 | !integer :: iband_uv(100) |
|---|
| 99 | !real(jprb) :: weight_uv(100) |
|---|
| 100 | ! Photosynthetically active radiation weightings |
|---|
| 101 | !integer :: nweight_par |
|---|
| 102 | !integer :: iband_par(100) |
|---|
| 103 | !real(jprb) :: weight_par(100) |
|---|
| 104 | |
|---|
| 105 | ! Loop index for repeats (for benchmarking) |
|---|
| 106 | integer :: jrepeat |
|---|
| 107 | |
|---|
| 108 | ! Are any variables out of bounds? |
|---|
| 109 | logical :: is_out_of_bounds |
|---|
| 110 | |
|---|
| 111 | ! integer :: iband(20), nweights |
|---|
| 112 | ! real(jprb) :: weight(20) |
|---|
| 113 | |
|---|
| 114 | |
|---|
| 115 | ! -------------------------------------------------------- |
|---|
| 116 | ! Section 2: Configure |
|---|
| 117 | ! -------------------------------------------------------- |
|---|
| 118 | |
|---|
| 119 | ! Check program called with correct number of arguments |
|---|
| 120 | if (command_argument_count() < 3) then |
|---|
| 121 | stop 'Usage: ecrad config.nam input_file.nc output_file.nc' |
|---|
| 122 | end if |
|---|
| 123 | |
|---|
| 124 | ! Use namelist to configure the radiation calculation |
|---|
| 125 | call get_command_argument(1, file_name, status=istatus) |
|---|
| 126 | if (istatus /= 0) then |
|---|
| 127 | stop 'Failed to read name of namelist file as string of length < 512' |
|---|
| 128 | end if |
|---|
| 129 | |
|---|
| 130 | ! Read "radiation" namelist into radiation configuration type |
|---|
| 131 | call config%read(file_name=file_name) |
|---|
| 132 | |
|---|
| 133 | ! Read "radiation_driver" namelist into radiation driver config type |
|---|
| 134 | call driver_config%read(file_name) |
|---|
| 135 | |
|---|
| 136 | if (driver_config%iverbose >= 2) then |
|---|
| 137 | write(nulout,'(a)') '-------------------------- OFFLINE ECRAD RADIATION SCHEME --------------------------' |
|---|
| 138 | write(nulout,'(a)') 'Copyright (C) 2014- ECMWF' |
|---|
| 139 | write(nulout,'(a)') 'Contact: Robin Hogan (r.j.hogan@ecmwf.int)' |
|---|
| 140 | #ifdef PARKIND1_SINGLE |
|---|
| 141 | write(nulout,'(a)') 'Floating-point precision: single' |
|---|
| 142 | #else |
|---|
| 143 | write(nulout,'(a)') 'Floating-point precision: double' |
|---|
| 144 | #endif |
|---|
| 145 | call config%print(driver_config%iverbose) |
|---|
| 146 | end if |
|---|
| 147 | |
|---|
| 148 | ! Albedo/emissivity intervals may be specified like this |
|---|
| 149 | !call config%define_sw_albedo_intervals(6, & |
|---|
| 150 | ! & [0.25e-6_jprb, 0.44e-6_jprb, 0.69e-6_jprb, & |
|---|
| 151 | ! & 1.19_jprb, 2.38e-6_jprb], [1,2,3,4,5,6], & |
|---|
| 152 | ! & do_nearest=.false.) |
|---|
| 153 | !call config%define_lw_emiss_intervals(3, & |
|---|
| 154 | ! & [8.0e-6_jprb, 13.0e-6_jprb], [1,2,1], & |
|---|
| 155 | ! & do_nearest=.false.) |
|---|
| 156 | |
|---|
| 157 | ! If monochromatic aerosol properties are required, then the |
|---|
| 158 | ! wavelengths can be specified (in metres) as follows - these can be |
|---|
| 159 | ! whatever you like for the general aerosol optics, but must match |
|---|
| 160 | ! the monochromatic values in the aerosol input file for the older |
|---|
| 161 | ! aerosol optics |
|---|
| 162 | !call config%set_aerosol_wavelength_mono( & |
|---|
| 163 | ! & [3.4e-07_jprb, 3.55e-07_jprb, 3.8e-07_jprb, 4.0e-07_jprb, 4.4e-07_jprb, & |
|---|
| 164 | ! & 4.69e-07_jprb, 5.0e-07_jprb, 5.32e-07_jprb, 5.5e-07_jprb, 6.45e-07_jprb, & |
|---|
| 165 | ! & 6.7e-07_jprb, 8.0e-07_jprb, 8.58e-07_jprb, 8.65e-07_jprb, 1.02e-06_jprb, & |
|---|
| 166 | ! & 1.064e-06_jprb, 1.24e-06_jprb, 1.64e-06_jprb, 2.13e-06_jprb, 1.0e-05_jprb]) |
|---|
| 167 | |
|---|
| 168 | ! Setup the radiation scheme: load the coefficients for gas and |
|---|
| 169 | ! cloud optics, currently from RRTMG |
|---|
| 170 | call setup_radiation(config) |
|---|
| 171 | |
|---|
| 172 | ! Demonstration of how to get weights for UV and PAR fluxes |
|---|
| 173 | !if (config%do_sw) then |
|---|
| 174 | ! call config%get_sw_weights(0.2e-6_jprb, 0.4415e-6_jprb,& |
|---|
| 175 | ! & nweight_uv, iband_uv, weight_uv,& |
|---|
| 176 | ! & 'ultraviolet') |
|---|
| 177 | ! call config%get_sw_weights(0.4e-6_jprb, 0.7e-6_jprb,& |
|---|
| 178 | ! & nweight_par, iband_par, weight_par,& |
|---|
| 179 | ! & 'photosynthetically active radiation, PAR') |
|---|
| 180 | !end if |
|---|
| 181 | |
|---|
| 182 | ! Optionally compute shortwave spectral diagnostics in |
|---|
| 183 | ! user-specified wavlength intervals |
|---|
| 184 | if (driver_config%n_sw_diag > 0) then |
|---|
| 185 | if (.not. config%do_surface_sw_spectral_flux) then |
|---|
| 186 | stop 'Error: shortwave spectral diagnostics require do_surface_sw_spectral_flux=true' |
|---|
| 187 | end if |
|---|
| 188 | call config%get_sw_mapping(driver_config%sw_diag_wavelength_bound(1:driver_config%n_sw_diag+1), & |
|---|
| 189 | & sw_diag_mapping, 'user-specified diagnostic intervals') |
|---|
| 190 | !if (driver_config%iverbose >= 3) then |
|---|
| 191 | ! call print_matrix(sw_diag_mapping, 'Shortwave diagnostic mapping', nulout) |
|---|
| 192 | !end if |
|---|
| 193 | end if |
|---|
| 194 | |
|---|
| 195 | if (driver_config%do_save_aerosol_optics) then |
|---|
| 196 | call config%aerosol_optics%save('aerosol_optics.nc', iverbose=driver_config%iverbose) |
|---|
| 197 | end if |
|---|
| 198 | |
|---|
| 199 | if (driver_config%do_save_cloud_optics .and. config%use_general_cloud_optics) then |
|---|
| 200 | call save_general_cloud_optics(config, 'hydrometeor_optics', iverbose=driver_config%iverbose) |
|---|
| 201 | end if |
|---|
| 202 | |
|---|
| 203 | ! -------------------------------------------------------- |
|---|
| 204 | ! Section 3: Read input data file |
|---|
| 205 | ! -------------------------------------------------------- |
|---|
| 206 | |
|---|
| 207 | ! Get NetCDF input file name |
|---|
| 208 | call get_command_argument(2, file_name, status=istatus) |
|---|
| 209 | if (istatus /= 0) then |
|---|
| 210 | stop 'Failed to read name of input NetCDF file as string of length < 512' |
|---|
| 211 | end if |
|---|
| 212 | |
|---|
| 213 | ! Open the file and configure the way it is read |
|---|
| 214 | call file%open(trim(file_name), iverbose=driver_config%iverbose) |
|---|
| 215 | |
|---|
| 216 | ! Get NetCDF output file name |
|---|
| 217 | call get_command_argument(3, file_name, status=istatus) |
|---|
| 218 | if (istatus /= 0) then |
|---|
| 219 | stop 'Failed to read name of output NetCDF file as string of length < 512' |
|---|
| 220 | end if |
|---|
| 221 | |
|---|
| 222 | ! 2D arrays are assumed to be stored in the file with height varying |
|---|
| 223 | ! more rapidly than column index. Specifying "true" here transposes |
|---|
| 224 | ! all 2D arrays so that the column index varies fastest within the |
|---|
| 225 | ! program. |
|---|
| 226 | call file%transpose_matrices(.true.) |
|---|
| 227 | |
|---|
| 228 | ! Read input variables from NetCDF file |
|---|
| 229 | call read_input(file, config, driver_config, ncol, nlev, & |
|---|
| 230 | & single_level, thermodynamics, & |
|---|
| 231 | & gas, cloud, aerosol) |
|---|
| 232 | |
|---|
| 233 | ! Close input file |
|---|
| 234 | call file%close() |
|---|
| 235 | |
|---|
| 236 | ! Compute seed from skin temperature residual |
|---|
| 237 | ! single_level%iseed = int(1.0e9*(single_level%skin_temperature & |
|---|
| 238 | ! & -int(single_level%skin_temperature))) |
|---|
| 239 | |
|---|
| 240 | ! Set first and last columns to process |
|---|
| 241 | if (driver_config%iendcol < 1 .or. driver_config%iendcol > ncol) then |
|---|
| 242 | driver_config%iendcol = ncol |
|---|
| 243 | end if |
|---|
| 244 | |
|---|
| 245 | if (driver_config%istartcol > driver_config%iendcol) then |
|---|
| 246 | write(nulout,'(a,i0,a,i0,a,i0,a)') '*** Error: requested column range (', & |
|---|
| 247 | & driver_config%istartcol, & |
|---|
| 248 | & ' to ', driver_config%iendcol, ') is out of the range in the data (1 to ', & |
|---|
| 249 | & ncol, ')' |
|---|
| 250 | stop 1 |
|---|
| 251 | end if |
|---|
| 252 | |
|---|
| 253 | ! Store inputs |
|---|
| 254 | if (driver_config%do_save_inputs) then |
|---|
| 255 | call save_inputs('inputs.nc', config, single_level, thermodynamics, & |
|---|
| 256 | & gas, cloud, aerosol, & |
|---|
| 257 | & lat=spread(0.0_jprb,1,ncol), & |
|---|
| 258 | & lon=spread(0.0_jprb,1,ncol), & |
|---|
| 259 | & iverbose=driver_config%iverbose) |
|---|
| 260 | end if |
|---|
| 261 | |
|---|
| 262 | ! -------------------------------------------------------- |
|---|
| 263 | ! Section 4: Call radiation scheme |
|---|
| 264 | ! -------------------------------------------------------- |
|---|
| 265 | |
|---|
| 266 | ! Ensure the units of the gas mixing ratios are what is required |
|---|
| 267 | ! by the gas absorption model |
|---|
| 268 | call set_gas_units(config, gas) |
|---|
| 269 | |
|---|
| 270 | ! Compute saturation with respect to liquid (needed for aerosol |
|---|
| 271 | ! hydration) call |
|---|
| 272 | call thermodynamics%calc_saturation_wrt_liquid(driver_config%istartcol,driver_config%iendcol) |
|---|
| 273 | |
|---|
| 274 | ! Check inputs are within physical bounds, printing message if not |
|---|
| 275 | is_out_of_bounds = gas%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 276 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 277 | & .or. single_level%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 278 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 279 | & .or. thermodynamics%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 280 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 281 | & .or. cloud%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 282 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 283 | & .or. aerosol%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 284 | & driver_config%do_correct_unphysical_inputs) |
|---|
| 285 | |
|---|
| 286 | ! Allocate memory for the flux profiles, which may include arrays |
|---|
| 287 | ! of dimension n_bands_sw/n_bands_lw, so must be called after |
|---|
| 288 | ! setup_radiation |
|---|
| 289 | call flux%allocate(config, 1, ncol, nlev) |
|---|
| 290 | |
|---|
| 291 | if (driver_config%iverbose >= 2) then |
|---|
| 292 | write(nulout,'(a)') 'Performing radiative transfer calculations' |
|---|
| 293 | end if |
|---|
| 294 | |
|---|
| 295 | ! Option of repeating calculation multiple time for more accurate |
|---|
| 296 | ! profiling |
|---|
| 297 | #ifndef NO_OPENMP |
|---|
| 298 | tstart = omp_get_wtime() |
|---|
| 299 | #endif |
|---|
| 300 | do jrepeat = 1,driver_config%nrepeat |
|---|
| 301 | |
|---|
| 302 | if (driver_config%do_parallel) then |
|---|
| 303 | ! Run radiation scheme over blocks of columns in parallel |
|---|
| 304 | |
|---|
| 305 | ! Compute number of blocks to process |
|---|
| 306 | nblock = (driver_config%iendcol - driver_config%istartcol & |
|---|
| 307 | & + driver_config%nblocksize) / driver_config%nblocksize |
|---|
| 308 | |
|---|
| 309 | !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME) |
|---|
| 310 | do jblock = 1, nblock |
|---|
| 311 | ! Specify the range of columns to process. |
|---|
| 312 | istartcol = (jblock-1) * driver_config%nblocksize & |
|---|
| 313 | & + driver_config%istartcol |
|---|
| 314 | iendcol = min(istartcol + driver_config%nblocksize - 1, & |
|---|
| 315 | & driver_config%iendcol) |
|---|
| 316 | |
|---|
| 317 | if (driver_config%iverbose >= 3) then |
|---|
| 318 | #ifndef NO_OPENMP |
|---|
| 319 | write(nulout,'(a,i0,a,i0,a,i0)') 'Thread ', omp_get_thread_num(), & |
|---|
| 320 | & ' processing columns ', istartcol, '-', iendcol |
|---|
| 321 | #else |
|---|
| 322 | write(nulout,'(a,i0,a,i0)') 'Processing columns ', istartcol, '-', iendcol |
|---|
| 323 | #endif |
|---|
| 324 | end if |
|---|
| 325 | |
|---|
| 326 | ! Call the ECRAD radiation scheme |
|---|
| 327 | call radiation(ncol, nlev, istartcol, iendcol, config, & |
|---|
| 328 | & single_level, thermodynamics, gas, cloud, aerosol, flux) |
|---|
| 329 | |
|---|
| 330 | end do |
|---|
| 331 | !$OMP END PARALLEL DO |
|---|
| 332 | |
|---|
| 333 | else |
|---|
| 334 | ! Run radiation scheme serially |
|---|
| 335 | if (driver_config%iverbose >= 3) then |
|---|
| 336 | write(nulout,'(a,i0,a)') 'Processing ', ncol, ' columns' |
|---|
| 337 | end if |
|---|
| 338 | |
|---|
| 339 | ! Call the ECRAD radiation scheme |
|---|
| 340 | call radiation(ncol, nlev, driver_config%istartcol, driver_config%iendcol, & |
|---|
| 341 | & config, single_level, thermodynamics, gas, cloud, aerosol, flux) |
|---|
| 342 | |
|---|
| 343 | end if |
|---|
| 344 | |
|---|
| 345 | end do |
|---|
| 346 | |
|---|
| 347 | #ifndef NO_OPENMP |
|---|
| 348 | tstop = omp_get_wtime() |
|---|
| 349 | write(nulout, '(a,g12.5,a)') 'Time elapsed in radiative transfer: ', tstop-tstart, ' seconds' |
|---|
| 350 | #endif |
|---|
| 351 | |
|---|
| 352 | ! -------------------------------------------------------- |
|---|
| 353 | ! Section 5: Check and save output |
|---|
| 354 | ! -------------------------------------------------------- |
|---|
| 355 | |
|---|
| 356 | is_out_of_bounds = flux%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol) |
|---|
| 357 | |
|---|
| 358 | ! Store the fluxes in the output file |
|---|
| 359 | if (.not. driver_config%do_save_net_fluxes) then |
|---|
| 360 | call save_fluxes(file_name, config, thermodynamics, flux, & |
|---|
| 361 | & iverbose=driver_config%iverbose, is_hdf5_file=driver_config%do_write_hdf5, & |
|---|
| 362 | & experiment_name=driver_config%experiment_name, & |
|---|
| 363 | & is_double_precision=driver_config%do_write_double_precision) |
|---|
| 364 | else |
|---|
| 365 | call save_net_fluxes(file_name, config, thermodynamics, flux, & |
|---|
| 366 | & iverbose=driver_config%iverbose, is_hdf5_file=driver_config%do_write_hdf5, & |
|---|
| 367 | & experiment_name=driver_config%experiment_name, & |
|---|
| 368 | & is_double_precision=driver_config%do_write_double_precision) |
|---|
| 369 | end if |
|---|
| 370 | |
|---|
| 371 | if (driver_config%n_sw_diag > 0) then |
|---|
| 372 | ! Store spectral fluxes in user-defined intervals in a second |
|---|
| 373 | ! output file |
|---|
| 374 | call save_sw_diagnostics(driver_config%sw_diag_file_name, config, & |
|---|
| 375 | & driver_config%sw_diag_wavelength_bound(1:driver_config%n_sw_diag+1), & |
|---|
| 376 | & sw_diag_mapping, flux, iverbose=driver_config%iverbose, & |
|---|
| 377 | & is_hdf5_file=driver_config%do_write_hdf5, & |
|---|
| 378 | & experiment_name=driver_config%experiment_name, & |
|---|
| 379 | & is_double_precision=driver_config%do_write_double_precision) |
|---|
| 380 | end if |
|---|
| 381 | |
|---|
| 382 | if (driver_config%iverbose >= 2) then |
|---|
| 383 | write(nulout,'(a)') '------------------------------------------------------------------------------------' |
|---|
| 384 | end if |
|---|
| 385 | |
|---|
| 386 | end program ecrad_driver |
|---|
