! ecrad_driver.F90 - Driver for offline ECRAD radiation scheme ! ! (C) Copyright 2014- ECMWF. ! ! This software is licensed under the terms of the Apache Licence Version 2.0 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. ! ! In applying this licence, ECMWF does not waive the privileges and immunities ! granted to it by virtue of its status as an intergovernmental organisation ! nor does it submit to any jurisdiction. ! ! Author: Robin Hogan ! Email: r.j.hogan@ecmwf.int ! ! ECRAD is the radiation scheme used in the ECMWF Integrated ! Forecasting System in cycle 43R3 and later. Several solvers are ! available, including McICA, Tripleclouds and SPARTACUS (the Speedy ! Algorithm for Radiative Transfer through Cloud Sides, a modification ! of the two-stream formulation of shortwave and longwave radiative ! transfer to account for 3D radiative effects). Gas optical ! properties are provided by the RRTM-G gas optics scheme. ! This program takes three arguments: ! 1) Namelist file to configure the radiation calculation ! 2) Name of a NetCDF file containing one or more atmospheric profiles ! 3) Name of output NetCDF file program ecrad_driver ! -------------------------------------------------------- ! Section 1: Declarations ! -------------------------------------------------------- use parkind1, only : jprb, jprd ! Working/double precision use radiation_io, only : nulout use radiation_interface, only : setup_radiation, radiation, set_gas_units use radiation_config, only : config_type use radiation_single_level, only : single_level_type use radiation_thermodynamics, only : thermodynamics_type use radiation_gas, only : gas_type, & & IVolumeMixingRatio, IMassMixingRatio, & & IH2O, ICO2, IO3, IN2O, ICO, ICH4, IO2, ICFC11, ICFC12, & & IHCFC22, ICCl4, GasName, GasLowerCaseName, NMaxGases use radiation_cloud, only : cloud_type use radiation_aerosol, only : aerosol_type use radiation_flux, only : flux_type use radiation_save, only : save_fluxes, save_net_fluxes, & & save_inputs, save_sw_diagnostics use radiation_general_cloud_optics, only : save_general_cloud_optics use ecrad_driver_config, only : driver_config_type use ecrad_driver_read_input, only : read_input use easy_netcdf use print_matrix_mod, only : print_matrix implicit none ! The NetCDF file containing the input profiles type(netcdf_file) :: file ! Derived types for the inputs to the radiation scheme type(config_type) :: config type(single_level_type) :: single_level type(thermodynamics_type) :: thermodynamics type(gas_type) :: gas type(cloud_type) :: cloud type(aerosol_type) :: aerosol ! Configuration specific to this driver type(driver_config_type) :: driver_config ! Derived type containing outputs from the radiation scheme type(flux_type) :: flux integer :: ncol, nlev ! Number of columns and levels integer :: istartcol, iendcol ! Range of columns to process ! Name of file names specified on command line character(len=512) :: file_name integer :: istatus ! Result of command_argument_count ! For parallel processing of multiple blocks integer :: jblock, nblock ! Block loop index and number ! Mapping matrix for shortwave spectral diagnostics real(jprb), allocatable :: sw_diag_mapping(:,:) #ifndef NO_OPENMP ! OpenMP functions integer, external :: omp_get_thread_num real(kind=jprd), external :: omp_get_wtime ! Start/stop time in seconds real(kind=jprd) :: tstart, tstop #endif ! For demonstration of get_sw_weights later on ! Ultraviolet weightings !integer :: nweight_uv !integer :: iband_uv(100) !real(jprb) :: weight_uv(100) ! Photosynthetically active radiation weightings !integer :: nweight_par !integer :: iband_par(100) !real(jprb) :: weight_par(100) ! Loop index for repeats (for benchmarking) integer :: jrepeat ! Are any variables out of bounds? logical :: is_out_of_bounds ! integer :: iband(20), nweights ! real(jprb) :: weight(20) ! -------------------------------------------------------- ! Section 2: Configure ! -------------------------------------------------------- ! Check program called with correct number of arguments if (command_argument_count() < 3) then stop 'Usage: ecrad config.nam input_file.nc output_file.nc' end if ! Use namelist to configure the radiation calculation call get_command_argument(1, file_name, status=istatus) if (istatus /= 0) then stop 'Failed to read name of namelist file as string of length < 512' end if ! Read "radiation" namelist into radiation configuration type call config%read(file_name=file_name) ! Read "radiation_driver" namelist into radiation driver config type call driver_config%read(file_name) if (driver_config%iverbose >= 2) then write(nulout,'(a)') '-------------------------- OFFLINE ECRAD RADIATION SCHEME --------------------------' write(nulout,'(a)') 'Copyright (C) 2014- ECMWF' write(nulout,'(a)') 'Contact: Robin Hogan (r.j.hogan@ecmwf.int)' #ifdef PARKIND1_SINGLE write(nulout,'(a)') 'Floating-point precision: single' #else write(nulout,'(a)') 'Floating-point precision: double' #endif call config%print(driver_config%iverbose) end if ! Albedo/emissivity intervals may be specified like this !call config%define_sw_albedo_intervals(6, & ! & [0.25e-6_jprb, 0.44e-6_jprb, 0.69e-6_jprb, & ! & 1.19_jprb, 2.38e-6_jprb], [1,2,3,4,5,6], & ! & do_nearest=.false.) !call config%define_lw_emiss_intervals(3, & ! & [8.0e-6_jprb, 13.0e-6_jprb], [1,2,1], & ! & do_nearest=.false.) ! If monochromatic aerosol properties are required, then the ! wavelengths can be specified (in metres) as follows - these can be ! whatever you like for the general aerosol optics, but must match ! the monochromatic values in the aerosol input file for the older ! aerosol optics !call config%set_aerosol_wavelength_mono( & ! & [3.4e-07_jprb, 3.55e-07_jprb, 3.8e-07_jprb, 4.0e-07_jprb, 4.4e-07_jprb, & ! & 4.69e-07_jprb, 5.0e-07_jprb, 5.32e-07_jprb, 5.5e-07_jprb, 6.45e-07_jprb, & ! & 6.7e-07_jprb, 8.0e-07_jprb, 8.58e-07_jprb, 8.65e-07_jprb, 1.02e-06_jprb, & ! & 1.064e-06_jprb, 1.24e-06_jprb, 1.64e-06_jprb, 2.13e-06_jprb, 1.0e-05_jprb]) ! Setup the radiation scheme: load the coefficients for gas and ! cloud optics, currently from RRTMG call setup_radiation(config) ! Demonstration of how to get weights for UV and PAR fluxes !if (config%do_sw) then ! call config%get_sw_weights(0.2e-6_jprb, 0.4415e-6_jprb,& ! & nweight_uv, iband_uv, weight_uv,& ! & 'ultraviolet') ! call config%get_sw_weights(0.4e-6_jprb, 0.7e-6_jprb,& ! & nweight_par, iband_par, weight_par,& ! & 'photosynthetically active radiation, PAR') !end if ! Optionally compute shortwave spectral diagnostics in ! user-specified wavlength intervals if (driver_config%n_sw_diag > 0) then if (.not. config%do_surface_sw_spectral_flux) then stop 'Error: shortwave spectral diagnostics require do_surface_sw_spectral_flux=true' end if call config%get_sw_mapping(driver_config%sw_diag_wavelength_bound(1:driver_config%n_sw_diag+1), & & sw_diag_mapping, 'user-specified diagnostic intervals') !if (driver_config%iverbose >= 3) then ! call print_matrix(sw_diag_mapping, 'Shortwave diagnostic mapping', nulout) !end if end if if (driver_config%do_save_aerosol_optics) then call config%aerosol_optics%save('aerosol_optics.nc', iverbose=driver_config%iverbose) end if if (driver_config%do_save_cloud_optics .and. config%use_general_cloud_optics) then call save_general_cloud_optics(config, 'hydrometeor_optics', iverbose=driver_config%iverbose) end if ! -------------------------------------------------------- ! Section 3: Read input data file ! -------------------------------------------------------- ! Get NetCDF input file name call get_command_argument(2, file_name, status=istatus) if (istatus /= 0) then stop 'Failed to read name of input NetCDF file as string of length < 512' end if ! Open the file and configure the way it is read call file%open(trim(file_name), iverbose=driver_config%iverbose) ! Get NetCDF output file name call get_command_argument(3, file_name, status=istatus) if (istatus /= 0) then stop 'Failed to read name of output NetCDF file as string of length < 512' end if ! 2D arrays are assumed to be stored in the file with height varying ! more rapidly than column index. Specifying "true" here transposes ! all 2D arrays so that the column index varies fastest within the ! program. call file%transpose_matrices(.true.) ! Read input variables from NetCDF file call read_input(file, config, driver_config, ncol, nlev, & & single_level, thermodynamics, & & gas, cloud, aerosol) ! Close input file call file%close() ! Compute seed from skin temperature residual ! single_level%iseed = int(1.0e9*(single_level%skin_temperature & ! & -int(single_level%skin_temperature))) ! Set first and last columns to process if (driver_config%iendcol < 1 .or. driver_config%iendcol > ncol) then driver_config%iendcol = ncol end if if (driver_config%istartcol > driver_config%iendcol) then write(nulout,'(a,i0,a,i0,a,i0,a)') '*** Error: requested column range (', & & driver_config%istartcol, & & ' to ', driver_config%iendcol, ') is out of the range in the data (1 to ', & & ncol, ')' stop 1 end if ! Store inputs if (driver_config%do_save_inputs) then call save_inputs('inputs.nc', config, single_level, thermodynamics, & & gas, cloud, aerosol, & & lat=spread(0.0_jprb,1,ncol), & & lon=spread(0.0_jprb,1,ncol), & & iverbose=driver_config%iverbose) end if ! -------------------------------------------------------- ! Section 4: Call radiation scheme ! -------------------------------------------------------- ! Ensure the units of the gas mixing ratios are what is required ! by the gas absorption model call set_gas_units(config, gas) ! Compute saturation with respect to liquid (needed for aerosol ! hydration) call call thermodynamics%calc_saturation_wrt_liquid(driver_config%istartcol,driver_config%iendcol) ! Check inputs are within physical bounds, printing message if not is_out_of_bounds = gas%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & & driver_config%do_correct_unphysical_inputs) & & .or. single_level%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & & driver_config%do_correct_unphysical_inputs) & & .or. thermodynamics%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & & driver_config%do_correct_unphysical_inputs) & & .or. cloud%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & & driver_config%do_correct_unphysical_inputs) & & .or. aerosol%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & & driver_config%do_correct_unphysical_inputs) ! Allocate memory for the flux profiles, which may include arrays ! of dimension n_bands_sw/n_bands_lw, so must be called after ! setup_radiation call flux%allocate(config, 1, ncol, nlev) if (driver_config%iverbose >= 2) then write(nulout,'(a)') 'Performing radiative transfer calculations' end if ! Option of repeating calculation multiple time for more accurate ! profiling #ifndef NO_OPENMP tstart = omp_get_wtime() #endif do jrepeat = 1,driver_config%nrepeat if (driver_config%do_parallel) then ! Run radiation scheme over blocks of columns in parallel ! Compute number of blocks to process nblock = (driver_config%iendcol - driver_config%istartcol & & + driver_config%nblocksize) / driver_config%nblocksize !$OMP PARALLEL DO PRIVATE(istartcol, iendcol) SCHEDULE(RUNTIME) do jblock = 1, nblock ! Specify the range of columns to process. istartcol = (jblock-1) * driver_config%nblocksize & & + driver_config%istartcol iendcol = min(istartcol + driver_config%nblocksize - 1, & & driver_config%iendcol) if (driver_config%iverbose >= 3) then #ifndef NO_OPENMP write(nulout,'(a,i0,a,i0,a,i0)') 'Thread ', omp_get_thread_num(), & & ' processing columns ', istartcol, '-', iendcol #else write(nulout,'(a,i0,a,i0)') 'Processing columns ', istartcol, '-', iendcol #endif end if ! Call the ECRAD radiation scheme call radiation(ncol, nlev, istartcol, iendcol, config, & & single_level, thermodynamics, gas, cloud, aerosol, flux) end do !$OMP END PARALLEL DO else ! Run radiation scheme serially if (driver_config%iverbose >= 3) then write(nulout,'(a,i0,a)') 'Processing ', ncol, ' columns' end if ! Call the ECRAD radiation scheme call radiation(ncol, nlev, driver_config%istartcol, driver_config%iendcol, & & config, single_level, thermodynamics, gas, cloud, aerosol, flux) end if end do #ifndef NO_OPENMP tstop = omp_get_wtime() write(nulout, '(a,g12.5,a)') 'Time elapsed in radiative transfer: ', tstop-tstart, ' seconds' #endif ! -------------------------------------------------------- ! Section 5: Check and save output ! -------------------------------------------------------- is_out_of_bounds = flux%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol) ! Store the fluxes in the output file if (.not. driver_config%do_save_net_fluxes) then call save_fluxes(file_name, config, thermodynamics, flux, & & iverbose=driver_config%iverbose, is_hdf5_file=driver_config%do_write_hdf5, & & experiment_name=driver_config%experiment_name, & & is_double_precision=driver_config%do_write_double_precision) else call save_net_fluxes(file_name, config, thermodynamics, flux, & & iverbose=driver_config%iverbose, is_hdf5_file=driver_config%do_write_hdf5, & & experiment_name=driver_config%experiment_name, & & is_double_precision=driver_config%do_write_double_precision) end if if (driver_config%n_sw_diag > 0) then ! Store spectral fluxes in user-defined intervals in a second ! output file call save_sw_diagnostics(driver_config%sw_diag_file_name, config, & & driver_config%sw_diag_wavelength_bound(1:driver_config%n_sw_diag+1), & & sw_diag_mapping, flux, iverbose=driver_config%iverbose, & & is_hdf5_file=driver_config%do_write_hdf5, & & experiment_name=driver_config%experiment_name, & & is_double_precision=driver_config%do_write_double_precision) end if if (driver_config%iverbose >= 2) then write(nulout,'(a)') '------------------------------------------------------------------------------------' end if end program ecrad_driver