| 1 | ! ecrad_ifs_driver_blocked.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, but note |
|---|
| 25 | ! that only the radiation_config group is read |
|---|
| 26 | ! 2) Name of a NetCDF file containing one or more atmospheric profiles |
|---|
| 27 | ! 3) Name of output NetCDF file |
|---|
| 28 | ! |
|---|
| 29 | ! This version uses the infrastructure of the IFS, such as computing |
|---|
| 30 | ! effective radius and cloud overlap from latitude and other |
|---|
| 31 | ! variables. To configure ecRad in this version you need to edit |
|---|
| 32 | ! ifs/yoerad.F90 in the ecRad package, but these options can be |
|---|
| 33 | ! overridden with the "radiation" namelist. This file requires the |
|---|
| 34 | ! input data to have compatible settings, e.g. the right number of |
|---|
| 35 | ! aerosol variables, and surface albedo/emissivity bands; a test file |
|---|
| 36 | ! satisfying this requirement is test/ifs/ecrad_meridian.nc in the |
|---|
| 37 | ! ecRad package. |
|---|
| 38 | ! |
|---|
| 39 | ! Note that the purpose of this file is simply to demonstrate the use |
|---|
| 40 | ! of the setup_radiation_scheme and radiation_scheme routines as well |
|---|
| 41 | ! as the use of a blocked memory layout to improve cache efficiency; |
|---|
| 42 | ! all the rest is using the offline ecRad driver containers to read |
|---|
| 43 | ! a NetCDF file to memory and pass it into these routines. |
|---|
| 44 | |
|---|
| 45 | program ecrad_ifs_driver |
|---|
| 46 | |
|---|
| 47 | ! -------------------------------------------------------- |
|---|
| 48 | ! Section 1: Declarations |
|---|
| 49 | ! -------------------------------------------------------- |
|---|
| 50 | use parkind1, only : jprb, jprd ! Working/double precision |
|---|
| 51 | |
|---|
| 52 | use radiation_io, only : nulout |
|---|
| 53 | use radiation_single_level, only : single_level_type |
|---|
| 54 | use radiation_thermodynamics, only : thermodynamics_type |
|---|
| 55 | use radiation_gas, only : gas_type, IMassMixingRatio, & |
|---|
| 56 | & IH2O, ICO2, IO3, IN2O, INO2, ICO, ICH4, IO2, ICFC11, ICFC12, & |
|---|
| 57 | & IHCFC22, ICCl4 |
|---|
| 58 | use radiation_cloud, only : cloud_type |
|---|
| 59 | use radiation_aerosol, only : aerosol_type |
|---|
| 60 | use radiation_flux, only : flux_type |
|---|
| 61 | use radiation_save, only : save_net_fluxes |
|---|
| 62 | use radiation_setup, only : tradiation, setup_radiation_scheme |
|---|
| 63 | use radiation_constants, only : Pi |
|---|
| 64 | use ecrad_driver_config, only : driver_config_type |
|---|
| 65 | use ecrad_driver_read_input, only : read_input |
|---|
| 66 | use easy_netcdf |
|---|
| 67 | use ifs_blocking |
|---|
| 68 | |
|---|
| 69 | implicit none |
|---|
| 70 | |
|---|
| 71 | #include "radiation_scheme.intfb.h" |
|---|
| 72 | |
|---|
| 73 | ! The NetCDF file containing the input profiles |
|---|
| 74 | type(netcdf_file) :: file |
|---|
| 75 | |
|---|
| 76 | ! Configuration for the radiation scheme, IFS style |
|---|
| 77 | type(tradiation) :: yradiation |
|---|
| 78 | |
|---|
| 79 | ! Derived types for the inputs to the radiation scheme |
|---|
| 80 | type(single_level_type) :: single_level |
|---|
| 81 | type(thermodynamics_type) :: thermodynamics |
|---|
| 82 | type(gas_type) :: gas |
|---|
| 83 | type(cloud_type) :: cloud |
|---|
| 84 | type(aerosol_type) :: aerosol |
|---|
| 85 | |
|---|
| 86 | ! Configuration specific to this driver |
|---|
| 87 | type(driver_config_type) :: driver_config |
|---|
| 88 | |
|---|
| 89 | ! Derived type containing outputs from the radiation scheme |
|---|
| 90 | type(flux_type) :: flux |
|---|
| 91 | |
|---|
| 92 | ! Additional arrays passed to radiation_scheme |
|---|
| 93 | real(jprb), allocatable, dimension(:) :: ccn_land, ccn_sea, sin_latitude, longitude_rad, land_frac |
|---|
| 94 | real(jprb), allocatable, dimension(:,:) :: pressure_fl, temperature_fl |
|---|
| 95 | real(jprb), allocatable, dimension(:) :: flux_sw_direct_normal, flux_uv, flux_par, flux_par_clear, & |
|---|
| 96 | & emissivity_out |
|---|
| 97 | real(jprb), allocatable, dimension(:,:) :: flux_diffuse_band, flux_direct_band |
|---|
| 98 | |
|---|
| 99 | ! Bespoke data types to set-up the blocked memory layout |
|---|
| 100 | type(ifs_config_type) :: ifs_config |
|---|
| 101 | real(kind=jprb), allocatable :: zrgp(:,:,:) ! monolithic IFS data structure |
|---|
| 102 | integer, allocatable :: iseed(:,:) ! Seed for random number generator |
|---|
| 103 | |
|---|
| 104 | integer :: ncol, nlev ! Number of columns and levels |
|---|
| 105 | integer :: nproma ! block size |
|---|
| 106 | |
|---|
| 107 | ! Name of file names specified on command line |
|---|
| 108 | character(len=512) :: file_name |
|---|
| 109 | integer :: istatus ! Result of command_argument_count |
|---|
| 110 | |
|---|
| 111 | #ifndef NO_OPENMP |
|---|
| 112 | ! OpenMP functions |
|---|
| 113 | integer, external :: omp_get_thread_num |
|---|
| 114 | real(kind=jprd), external :: omp_get_wtime |
|---|
| 115 | ! Start/stop time in seconds |
|---|
| 116 | real(kind=jprd) :: tstart, tstop |
|---|
| 117 | #endif |
|---|
| 118 | |
|---|
| 119 | ! For demonstration of get_sw_weights later on |
|---|
| 120 | ! Ultraviolet weightings |
|---|
| 121 | !integer :: nweight_uv |
|---|
| 122 | !integer :: iband_uv(100) |
|---|
| 123 | !real(jprb) :: weight_uv(100) |
|---|
| 124 | ! Photosynthetically active radiation weightings |
|---|
| 125 | !integer :: nweight_par |
|---|
| 126 | !integer :: iband_par(100) |
|---|
| 127 | !real(jprb) :: weight_par(100) |
|---|
| 128 | |
|---|
| 129 | ! Loop index for repeats (for benchmarking) |
|---|
| 130 | integer :: jrepeat |
|---|
| 131 | |
|---|
| 132 | ! Loop index |
|---|
| 133 | integer :: jrl, ibeg, iend, il, ib |
|---|
| 134 | |
|---|
| 135 | ! Are any variables out of bounds? |
|---|
| 136 | logical :: is_out_of_bounds |
|---|
| 137 | |
|---|
| 138 | ! integer :: iband(20), nweights |
|---|
| 139 | ! real(jprb) :: weight(20) |
|---|
| 140 | |
|---|
| 141 | |
|---|
| 142 | ! -------------------------------------------------------- |
|---|
| 143 | ! Section 2: Configure |
|---|
| 144 | ! -------------------------------------------------------- |
|---|
| 145 | |
|---|
| 146 | ! Check program called with correct number of arguments |
|---|
| 147 | if (command_argument_count() < 3) then |
|---|
| 148 | stop 'Usage: ecrad config.nam input_file.nc output_file.nc' |
|---|
| 149 | end if |
|---|
| 150 | |
|---|
| 151 | ! Use namelist to configure the radiation calculation |
|---|
| 152 | call get_command_argument(1, file_name, status=istatus) |
|---|
| 153 | if (istatus /= 0) then |
|---|
| 154 | stop 'Failed to read name of namelist file as string of length < 512' |
|---|
| 155 | end if |
|---|
| 156 | |
|---|
| 157 | ! Read "radiation_driver" namelist into radiation driver config type |
|---|
| 158 | call driver_config%read(file_name) |
|---|
| 159 | nproma = driver_config%nblocksize |
|---|
| 160 | |
|---|
| 161 | if (driver_config%iverbose >= 2) then |
|---|
| 162 | write(nulout,'(a)') '-------------------------- OFFLINE ECRAD RADIATION SCHEME --------------------------' |
|---|
| 163 | write(nulout,'(a)') 'Copyright (C) 2014- ECMWF' |
|---|
| 164 | write(nulout,'(a)') 'Contact: Robin Hogan (r.j.hogan@ecmwf.int)' |
|---|
| 165 | #ifdef PARKIND1_SINGLE |
|---|
| 166 | write(nulout,'(a)') 'Floating-point precision: single' |
|---|
| 167 | #else |
|---|
| 168 | write(nulout,'(a)') 'Floating-point precision: double' |
|---|
| 169 | #endif |
|---|
| 170 | end if |
|---|
| 171 | |
|---|
| 172 | ! Albedo/emissivity intervals may be specified like this |
|---|
| 173 | !call config%define_sw_albedo_intervals(6, & |
|---|
| 174 | ! & [0.25e-6_jprb, 0.44e-6_jprb, 0.69e-6_jprb, & |
|---|
| 175 | ! & 1.19_jprb, 2.38e-6_jprb], [1,2,3,4,5,6], & |
|---|
| 176 | ! & do_nearest=.false.) |
|---|
| 177 | !call config%define_lw_emiss_intervals(3, & |
|---|
| 178 | ! & [8.0e-6_jprb, 13.0e-6_jprb], [1,2,1], & |
|---|
| 179 | ! & do_nearest=.false.) |
|---|
| 180 | |
|---|
| 181 | ! If monochromatic aerosol properties are required, then the |
|---|
| 182 | ! wavelengths can be specified (in metres) as follows - these can be |
|---|
| 183 | ! whatever you like for the general aerosol optics, but must match |
|---|
| 184 | ! the monochromatic values in the aerosol input file for the older |
|---|
| 185 | ! aerosol optics |
|---|
| 186 | !call config%set_aerosol_wavelength_mono( & |
|---|
| 187 | ! & [3.4e-07_jprb, 3.55e-07_jprb, 3.8e-07_jprb, 4.0e-07_jprb, 4.4e-07_jprb, & |
|---|
| 188 | ! & 4.69e-07_jprb, 5.0e-07_jprb, 5.32e-07_jprb, 5.5e-07_jprb, 6.45e-07_jprb, & |
|---|
| 189 | ! & 6.7e-07_jprb, 8.0e-07_jprb, 8.58e-07_jprb, 8.65e-07_jprb, 1.02e-06_jprb, & |
|---|
| 190 | ! & 1.064e-06_jprb, 1.24e-06_jprb, 1.64e-06_jprb, 2.13e-06_jprb, 1.0e-05_jprb]) |
|---|
| 191 | |
|---|
| 192 | call yradiation%rad_config%read(file_name=file_name) |
|---|
| 193 | |
|---|
| 194 | ! Setup aerosols |
|---|
| 195 | if (yradiation%rad_config%use_aerosols) then |
|---|
| 196 | yradiation%yrerad%naermacc = 1 ! MACC-derived aerosol climatology on a NMCLAT x NMCLON grid |
|---|
| 197 | else |
|---|
| 198 | yradiation%yrerad%naermacc = 0 |
|---|
| 199 | endif |
|---|
| 200 | |
|---|
| 201 | ! Setup the radiation scheme: load the coefficients for gas and |
|---|
| 202 | ! cloud optics, currently from RRTMG |
|---|
| 203 | call setup_radiation_scheme(yradiation, .true., file_name=file_name) |
|---|
| 204 | ! Or call without specifying the namelist filename, in which case |
|---|
| 205 | ! the default settings are from yoerad.F90 |
|---|
| 206 | !call setup_radiation_scheme(yradiation, .true.) |
|---|
| 207 | |
|---|
| 208 | ! Demonstration of how to get weights for UV and PAR fluxes |
|---|
| 209 | !if (config%do_sw) then |
|---|
| 210 | ! call config%get_sw_weights(0.2e-6_jprb, 0.4415e-6_jprb,& |
|---|
| 211 | ! & nweight_uv, iband_uv, weight_uv,& |
|---|
| 212 | ! & 'ultraviolet') |
|---|
| 213 | ! call config%get_sw_weights(0.4e-6_jprb, 0.7e-6_jprb,& |
|---|
| 214 | ! & nweight_par, iband_par, weight_par,& |
|---|
| 215 | ! & 'photosynthetically active radiation, PAR') |
|---|
| 216 | !end if |
|---|
| 217 | |
|---|
| 218 | ! -------------------------------------------------------- |
|---|
| 219 | ! Section 3: Read input data file |
|---|
| 220 | ! -------------------------------------------------------- |
|---|
| 221 | |
|---|
| 222 | ! Get NetCDF input file name |
|---|
| 223 | call get_command_argument(2, file_name, status=istatus) |
|---|
| 224 | if (istatus /= 0) then |
|---|
| 225 | stop 'Failed to read name of input NetCDF file as string of length < 512' |
|---|
| 226 | end if |
|---|
| 227 | |
|---|
| 228 | ! Open the file and configure the way it is read |
|---|
| 229 | call file%open(trim(file_name), iverbose=driver_config%iverbose) |
|---|
| 230 | |
|---|
| 231 | ! Get NetCDF output file name |
|---|
| 232 | call get_command_argument(3, file_name, status=istatus) |
|---|
| 233 | if (istatus /= 0) then |
|---|
| 234 | stop 'Failed to read name of output NetCDF file as string of length < 512' |
|---|
| 235 | end if |
|---|
| 236 | |
|---|
| 237 | ! 2D arrays are assumed to be stored in the file with height varying |
|---|
| 238 | ! more rapidly than column index. Specifying "true" here transposes |
|---|
| 239 | ! all 2D arrays so that the column index varies fastest within the |
|---|
| 240 | ! program. |
|---|
| 241 | call file%transpose_matrices(.true.) |
|---|
| 242 | |
|---|
| 243 | ! Read input variables from NetCDF file, noting that cloud overlap |
|---|
| 244 | ! and effective radius are ignored |
|---|
| 245 | call read_input(file, yradiation%rad_config, driver_config, ncol, nlev, & |
|---|
| 246 | & single_level, thermodynamics, & |
|---|
| 247 | & gas, cloud, aerosol) |
|---|
| 248 | |
|---|
| 249 | ! Latitude is used for cloud overlap and ice effective radius |
|---|
| 250 | if (file%exists('lat')) then |
|---|
| 251 | call file%get('lat', sin_latitude) |
|---|
| 252 | sin_latitude = sin(sin_latitude * Pi/180.0_jprb) |
|---|
| 253 | else |
|---|
| 254 | allocate(sin_latitude(ncol)) |
|---|
| 255 | sin_latitude = 0.0_jprb |
|---|
| 256 | end if |
|---|
| 257 | |
|---|
| 258 | if (file%exists('lon')) then |
|---|
| 259 | call file%get('lon', longitude_rad) |
|---|
| 260 | longitude_rad = longitude_rad * Pi/180.0_jprb |
|---|
| 261 | else |
|---|
| 262 | allocate(longitude_rad(ncol)) |
|---|
| 263 | longitude_rad = 0.0_jprb |
|---|
| 264 | end if |
|---|
| 265 | |
|---|
| 266 | ! Close input file |
|---|
| 267 | call file%close() |
|---|
| 268 | |
|---|
| 269 | ! Convert gas units to mass-mixing ratio |
|---|
| 270 | call gas%set_units(IMassMixingRatio) |
|---|
| 271 | |
|---|
| 272 | ! Compute seed from skin temperature residual |
|---|
| 273 | ! single_level%iseed = int(1.0e9*(single_level%skin_temperature & |
|---|
| 274 | ! & -int(single_level%skin_temperature))) |
|---|
| 275 | |
|---|
| 276 | ! Set first and last columns to process |
|---|
| 277 | if (driver_config%iendcol < 1 .or. driver_config%iendcol > ncol) then |
|---|
| 278 | driver_config%iendcol = ncol |
|---|
| 279 | end if |
|---|
| 280 | |
|---|
| 281 | if (driver_config%istartcol > driver_config%iendcol) then |
|---|
| 282 | write(nulout,'(a,i0,a,i0,a,i0,a)') '*** Error: requested column range (', & |
|---|
| 283 | & driver_config%istartcol, & |
|---|
| 284 | & ' to ', driver_config%iendcol, ') is out of the range in the data (1 to ', & |
|---|
| 285 | & ncol, ')' |
|---|
| 286 | stop 1 |
|---|
| 287 | end if |
|---|
| 288 | |
|---|
| 289 | ! -------------------------------------------------------- |
|---|
| 290 | ! Section 4: Call radiation scheme |
|---|
| 291 | ! -------------------------------------------------------- |
|---|
| 292 | |
|---|
| 293 | ! Compute saturation with respect to liquid (needed for aerosol |
|---|
| 294 | ! hydration) call |
|---|
| 295 | ! call thermodynamics%calc_saturation_wrt_liquid(driver_config%istartcol,driver_config%iendcol) |
|---|
| 296 | |
|---|
| 297 | ! Check inputs are within physical bounds, printing message if not |
|---|
| 298 | is_out_of_bounds = gas%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 299 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 300 | & .or. single_level%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 301 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 302 | & .or. thermodynamics%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 303 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 304 | & .or. cloud%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 305 | & driver_config%do_correct_unphysical_inputs) & |
|---|
| 306 | & .or. aerosol%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol, & |
|---|
| 307 | & driver_config%do_correct_unphysical_inputs) |
|---|
| 308 | |
|---|
| 309 | ! Allocate memory for the flux profiles, which may include arrays |
|---|
| 310 | ! of dimension n_bands_sw/n_bands_lw, so must be called after |
|---|
| 311 | ! setup_radiation |
|---|
| 312 | call flux%allocate(yradiation%rad_config, 1, ncol, nlev) |
|---|
| 313 | |
|---|
| 314 | ! set relevant fluxes to zero |
|---|
| 315 | flux%lw_up(:,:) = 0._jprb |
|---|
| 316 | flux%lw_dn(:,:) = 0._jprb |
|---|
| 317 | flux%sw_up(:,:) = 0._jprb |
|---|
| 318 | flux%sw_dn(:,:) = 0._jprb |
|---|
| 319 | flux%sw_dn_direct(:,:) = 0._jprb |
|---|
| 320 | flux%lw_up_clear(:,:) = 0._jprb |
|---|
| 321 | flux%lw_dn_clear(:,:) = 0._jprb |
|---|
| 322 | flux%sw_up_clear(:,:) = 0._jprb |
|---|
| 323 | flux%sw_dn_clear(:,:) = 0._jprb |
|---|
| 324 | flux%sw_dn_direct_clear(:,:) = 0._jprb |
|---|
| 325 | |
|---|
| 326 | flux%lw_dn_surf_canopy(:,:) = 0._jprb |
|---|
| 327 | flux%sw_dn_diffuse_surf_canopy(:,:) = 0._jprb |
|---|
| 328 | flux%sw_dn_direct_surf_canopy(:,:) = 0._jprb |
|---|
| 329 | flux%lw_derivatives(:,:) = 0._jprb |
|---|
| 330 | |
|---|
| 331 | ! Allocate memory for additional arrays |
|---|
| 332 | allocate(ccn_land(ncol)) |
|---|
| 333 | allocate(ccn_sea(ncol)) |
|---|
| 334 | allocate(land_frac(ncol)) |
|---|
| 335 | allocate(pressure_fl(ncol,nlev)) |
|---|
| 336 | allocate(temperature_fl(ncol,nlev)) |
|---|
| 337 | allocate(flux_sw_direct_normal(ncol)) |
|---|
| 338 | allocate(flux_uv(ncol)) |
|---|
| 339 | allocate(flux_par(ncol)) |
|---|
| 340 | allocate(flux_par_clear(ncol)) |
|---|
| 341 | allocate(emissivity_out(ncol)) |
|---|
| 342 | allocate(flux_diffuse_band(ncol,yradiation%yrerad%nsw)) |
|---|
| 343 | allocate(flux_direct_band(ncol,yradiation%yrerad%nsw)) |
|---|
| 344 | |
|---|
| 345 | ccn_land = yradiation%yrerad%rccnlnd |
|---|
| 346 | ccn_sea = yradiation%yrerad%rccnsea |
|---|
| 347 | pressure_fl = 0.5_jprb * (thermodynamics%pressure_hl(:,1:nlev)+thermodynamics%pressure_hl(:,2:nlev+1)) |
|---|
| 348 | temperature_fl = 0.5_jprb * (thermodynamics%temperature_hl(:,1:nlev)+thermodynamics%temperature_hl(:,2:nlev+1)) |
|---|
| 349 | |
|---|
| 350 | ! -------------------------------------------------------- |
|---|
| 351 | ! Section 4a: Reshuffle into blocked memory layout |
|---|
| 352 | ! -------------------------------------------------------- |
|---|
| 353 | |
|---|
| 354 | call ifs_setup_indices(driver_config, ifs_config, yradiation, nlev) |
|---|
| 355 | call ifs_copy_inputs_to_blocked(driver_config, ifs_config, yradiation,& |
|---|
| 356 | & ncol, nlev, single_level, thermodynamics, gas, cloud, aerosol,& |
|---|
| 357 | & sin_latitude, longitude_rad, land_frac, pressure_fl, temperature_fl,& |
|---|
| 358 | & zrgp, iseed=iseed) |
|---|
| 359 | |
|---|
| 360 | ! -------------------------------------------------------- |
|---|
| 361 | ! Section 4b: Call radiation_scheme with blocked memory data |
|---|
| 362 | ! -------------------------------------------------------- |
|---|
| 363 | |
|---|
| 364 | if (driver_config%iverbose >= 2) then |
|---|
| 365 | write(nulout,'(a)') 'Performing radiative transfer calculations' |
|---|
| 366 | end if |
|---|
| 367 | |
|---|
| 368 | ! Option of repeating calculation multiple time for more accurate |
|---|
| 369 | ! profiling |
|---|
| 370 | #ifndef NO_OPENMP |
|---|
| 371 | tstart = omp_get_wtime() |
|---|
| 372 | #endif |
|---|
| 373 | do jrepeat = 1,driver_config%nrepeat |
|---|
| 374 | |
|---|
| 375 | ! if (driver_config%do_parallel) then |
|---|
| 376 | ! Run radiation scheme over blocks of columns in parallel |
|---|
| 377 | |
|---|
| 378 | !$OMP PARALLEL DO SCHEDULE(DYNAMIC,1)& |
|---|
| 379 | !$OMP&PRIVATE(JRL,IBEG,IEND,IL,IB) |
|---|
| 380 | do jrl=1,ncol,nproma |
|---|
| 381 | ibeg=jrl |
|---|
| 382 | iend=min(ibeg+nproma-1,ncol) |
|---|
| 383 | il=iend-ibeg+1 |
|---|
| 384 | ib=(jrl-1)/nproma+1 |
|---|
| 385 | |
|---|
| 386 | if (driver_config%iverbose >= 3) then |
|---|
| 387 | #ifndef NO_OPENMP |
|---|
| 388 | write(nulout,'(a,i0,a,i0,a,i0)') 'Thread ', omp_get_thread_num(), & |
|---|
| 389 | & ' processing columns ', ibeg, '-', iend |
|---|
| 390 | #else |
|---|
| 391 | write(nulout,'(a,i0,a,i0)') 'Processing columns ', ibeg, '-', iend |
|---|
| 392 | #endif |
|---|
| 393 | end if |
|---|
| 394 | |
|---|
| 395 | ! Call the ECRAD radiation scheme |
|---|
| 396 | call radiation_scheme & |
|---|
| 397 | & (yradiation, & |
|---|
| 398 | & 1, il, nproma, & ! startcol, endcol, ncol |
|---|
| 399 | & nlev, size(aerosol%mixing_ratio,3), & ! nlev, naerosols |
|---|
| 400 | & single_level%solar_irradiance, & ! solar_irrad |
|---|
| 401 | ! array inputs |
|---|
| 402 | & zrgp(1,ifs_config%iamu0,ib), zrgp(1,ifs_config%its,ib), & ! mu0, skintemp |
|---|
| 403 | & zrgp(1,ifs_config%iald,ib) , zrgp(1,ifs_config%ialp,ib), & ! albedo_dif, albedo_dir |
|---|
| 404 | & zrgp(1,ifs_config%iemiss,ib), & ! spectral emissivity |
|---|
| 405 | & zrgp(1,ifs_config%iccnl,ib), zrgp(1,ifs_config%iccno,ib) ,& ! CCN concentration, land and sea |
|---|
| 406 | & zrgp(1,ifs_config%igelam,ib),zrgp(1,ifs_config%igemu,ib), & ! longitude, sine of latitude |
|---|
| 407 | & zrgp(1,ifs_config%islm,ib), & ! land sea mask |
|---|
| 408 | & zrgp(1,ifs_config%ipr,ib), zrgp(1,ifs_config%iti,ib), & ! full level pressure and temperature |
|---|
| 409 | & zrgp(1,ifs_config%iaprs,ib), zrgp(1,ifs_config%ihti,ib), & ! half-level pressure and temperature |
|---|
| 410 | & zrgp(1,ifs_config%iwv,ib), zrgp(1,ifs_config%iico2,ib), & |
|---|
| 411 | & zrgp(1,ifs_config%iich4,ib), zrgp(1,ifs_config%iin2o,ib), & |
|---|
| 412 | & zrgp(1,ifs_config%ino2,ib), zrgp(1,ifs_config%ic11,ib), & |
|---|
| 413 | & zrgp(1,ifs_config%ic12,ib), zrgp(1,ifs_config%ic22,ib), & |
|---|
| 414 | & zrgp(1,ifs_config%icl4,ib), zrgp(1,ifs_config%ioz,ib), & |
|---|
| 415 | & zrgp(1,ifs_config%iclc,ib), zrgp(1,ifs_config%ilwa,ib), & |
|---|
| 416 | & zrgp(1,ifs_config%iiwa,ib), zrgp(1,ifs_config%irwa,ib), & |
|---|
| 417 | & zrgp(1,ifs_config%iswa,ib), & |
|---|
| 418 | & zrgp(1,ifs_config%iaer,ib), zrgp(1,ifs_config%iaero,ib), & |
|---|
| 419 | ! flux outputs |
|---|
| 420 | & zrgp(1,ifs_config%ifrso,ib), zrgp(1,ifs_config%ifrth,ib), & |
|---|
| 421 | & zrgp(1,ifs_config%iswfc,ib), zrgp(1,ifs_config%ilwfc,ib),& |
|---|
| 422 | & zrgp(1,ifs_config%ifrsod,ib),zrgp(1,ifs_config%ifrted,ib), & |
|---|
| 423 | & zrgp(1,ifs_config%ifrsodc,ib),zrgp(1,ifs_config%ifrtedc,ib),& |
|---|
| 424 | & zrgp(1,ifs_config%ifdir,ib), zrgp(1,ifs_config%icdir,ib), & |
|---|
| 425 | & zrgp(1,ifs_config%isudu,ib), & |
|---|
| 426 | & zrgp(1,ifs_config%iuvdf,ib), zrgp(1,ifs_config%iparf,ib), & |
|---|
| 427 | & zrgp(1,ifs_config%iparcf,ib),zrgp(1,ifs_config%itincf,ib), & |
|---|
| 428 | & zrgp(1,ifs_config%iemit,ib) ,zrgp(1,ifs_config%ilwderivative,ib), & |
|---|
| 429 | & zrgp(1,ifs_config%iswdiffuseband,ib), zrgp(1,ifs_config%iswdirectband,ib)& |
|---|
| 430 | & ) |
|---|
| 431 | end do |
|---|
| 432 | !$OMP END PARALLEL DO |
|---|
| 433 | |
|---|
| 434 | ! else |
|---|
| 435 | ! Run radiation scheme serially |
|---|
| 436 | ! if (driver_config%iverbose >= 3) then |
|---|
| 437 | ! write(nulout,'(a,i0,a)') 'Processing ', ncol, ' columns' |
|---|
| 438 | ! end if |
|---|
| 439 | |
|---|
| 440 | ! Call the ECRAD radiation scheme |
|---|
| 441 | ! call radiation_scheme(ncol, nlev, driver_config%istartcol, driver_config%iendcol, & |
|---|
| 442 | ! & config, single_level, thermodynamics, gas, cloud, aerosol, flux) |
|---|
| 443 | |
|---|
| 444 | ! end if |
|---|
| 445 | |
|---|
| 446 | end do |
|---|
| 447 | |
|---|
| 448 | #ifndef NO_OPENMP |
|---|
| 449 | tstop = omp_get_wtime() |
|---|
| 450 | write(nulout, '(a,g12.5,a)') 'Time elapsed in radiative transfer: ', tstop-tstart, ' seconds' |
|---|
| 451 | #endif |
|---|
| 452 | |
|---|
| 453 | ! -------------------------------------------------------- |
|---|
| 454 | ! Section 4c: Copy fluxes from blocked memory data |
|---|
| 455 | ! -------------------------------------------------------- |
|---|
| 456 | |
|---|
| 457 | call ifs_copy_fluxes_from_blocked(driver_config, ifs_config, yradiation, ncol, nlev,& |
|---|
| 458 | & zrgp, flux, flux_sw_direct_normal, flux_uv, flux_par, flux_par_clear, & |
|---|
| 459 | & emissivity_out, flux_diffuse_band, flux_direct_band) |
|---|
| 460 | |
|---|
| 461 | ! "up" fluxes are actually net fluxes at this point - we modify the |
|---|
| 462 | ! upwelling flux so that net=dn-up, while the TOA and surface |
|---|
| 463 | ! downwelling fluxes are correct. |
|---|
| 464 | flux%sw_up = -flux%sw_up |
|---|
| 465 | flux%sw_up(:,1) = flux%sw_up(:,1)+flux%sw_dn(:,1) |
|---|
| 466 | flux%sw_up(:,nlev+1) = flux%sw_up(:,nlev+1)+flux%sw_dn(:,nlev+1) |
|---|
| 467 | |
|---|
| 468 | flux%lw_up = -flux%lw_up |
|---|
| 469 | flux%lw_up(:,1) = flux%lw_up(:,1)+flux%lw_dn(:,1) |
|---|
| 470 | flux%lw_up(:,nlev+1) = flux%lw_up(:,nlev+1)+flux%lw_dn(:,nlev+1) |
|---|
| 471 | |
|---|
| 472 | flux%sw_up_clear = -flux%sw_up_clear |
|---|
| 473 | flux%sw_up_clear(:,1) = flux%sw_up_clear(:,1)+flux%sw_dn_clear(:,1) |
|---|
| 474 | flux%sw_up_clear(:,nlev+1) = flux%sw_up_clear(:,nlev+1)+flux%sw_dn_clear(:,nlev+1) |
|---|
| 475 | |
|---|
| 476 | flux%lw_up_clear = -flux%lw_up_clear |
|---|
| 477 | flux%lw_up_clear(:,1) = flux%lw_up_clear(:,1)+flux%lw_dn_clear(:,1) |
|---|
| 478 | flux%lw_up_clear(:,nlev+1) = flux%lw_up_clear(:,nlev+1)+flux%lw_dn_clear(:,nlev+1) |
|---|
| 479 | |
|---|
| 480 | ! -------------------------------------------------------- |
|---|
| 481 | ! Section 5: Check and save output |
|---|
| 482 | ! -------------------------------------------------------- |
|---|
| 483 | |
|---|
| 484 | ! This is unreliable because only the net fluxes are valid: |
|---|
| 485 | !is_out_of_bounds = flux%out_of_physical_bounds(driver_config%istartcol, driver_config%iendcol) |
|---|
| 486 | |
|---|
| 487 | ! Store the fluxes in the output file |
|---|
| 488 | yradiation%rad_config%do_surface_sw_spectral_flux = .false. |
|---|
| 489 | yradiation%rad_config%do_canopy_fluxes_sw = .false. |
|---|
| 490 | yradiation%rad_config%do_canopy_fluxes_lw = .false. |
|---|
| 491 | |
|---|
| 492 | call save_net_fluxes(file_name, yradiation%rad_config, thermodynamics, flux, & |
|---|
| 493 | & iverbose=driver_config%iverbose, is_hdf5_file=driver_config%do_write_hdf5, & |
|---|
| 494 | & experiment_name=driver_config%experiment_name, & |
|---|
| 495 | & is_double_precision=driver_config%do_write_double_precision) |
|---|
| 496 | |
|---|
| 497 | if (driver_config%iverbose >= 2) then |
|---|
| 498 | write(nulout,'(a)') '------------------------------------------------------------------------------------' |
|---|
| 499 | end if |
|---|
| 500 | |
|---|
| 501 | end program ecrad_ifs_driver |
|---|
