source: LMDZ6/branches/cirrus/libf/phylmd/ecrad/driver/ecrad_driver.F90

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

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

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