source: LMDZ6/trunk/libf/phylmd/ecrad/radiation_ifs_rrtm.F90 @ 4660

Last change on this file since 4660 was 4489, checked in by lguez, 21 months ago

Merge LMDZ_ECRad branch back into trunk!

File size: 35.8 KB
Line 
1! radiation_ifs_rrtm.F90 - Interface to IFS implementation of RRTM-G
2!
3! (C) Copyright 2015- 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-11  R. Hogan  Receive "surface" dummy argument
17!   2017-09-08  R. Hogan  Reverted some changes
18!   2017-10-18  R. Hogan  Added planck_function public function
19!   2018-01-11  R. Hogan  Added optional spectral scaling of incoming solar radiation
20!   2018-02-22  R. Hogan  Optimized reverse indexing of heights
21!   2018-05-05  R. Hogan  gas_optics can be called for reduced number of levels
22!   2019-01-02  R. Hogan  Initialize shortwave props to zero in case sun below horizon
23
24module radiation_ifs_rrtm
25
26  implicit none
27
28  public  :: setup_gas_optics, gas_optics, planck_function, set_gas_units
29
30contains
31
32  !---------------------------------------------------------------------
33  ! Setup the IFS implementation of RRTM-G gas absorption model
34  subroutine setup_gas_optics(config, directory)
35
36    use yoerrtm,   only : jpglw
37    use yoesrtm,   only : jpgsw
38    use yoerrtftr, only : ngb_lw => ngb
39    use yoesrtm,   only : ngb_sw => ngbsw
40    use yomhook,   only : lhook, dr_hook
41
42    use radiation_config
43
44    type(config_type), intent(inout), target :: config
45    character(len=*), intent(in)     :: directory
46
47    integer :: irep ! For implied do
48
49    integer, parameter :: RRTM_GPOINT_REORDERING_LW(140) = (/ &
50          &   89, 90, 139, 77, 137, 69, 131, 97, 91, 70, 78, 71, 53, 72, 123, 54, 79, 98,  &
51          &   92, 55, 80, 132, 124, 81, 73, 56, 99, 82, 57, 23, 125, 100, 24, 74, 93, 58, 25,  &
52          &   83, 126, 75, 26, 11, 101, 133, 59, 27, 76, 140, 12, 84, 102, 94, 28, 127, 85,  &
53          &   13, 39, 60, 86, 103, 87, 109, 14, 29, 115, 40, 95, 15, 61, 88, 41, 110, 104, 1,  &
54          &   116, 42, 30, 134, 128, 138, 96, 62, 16, 43, 117, 63, 111, 44, 2, 64, 31, 65,  &
55          &   105, 17, 45, 66, 118, 32, 3, 33, 67, 18, 129, 135, 46, 112, 34, 106, 68, 35, 4,  &
56          &   119, 36, 47, 107, 19, 37, 38, 113, 48, 130, 5, 120, 49, 108, 20, 50, 51, 114,  &
57          &   21, 121, 52, 136, 122, 6, 22, 7, 8, 9, 10 &
58          & /)
59    integer, parameter :: RRTM_GPOINT_REORDERING_SW(112) = (/ &
60          &   35, 45, 19, 27, 36, 57, 20, 46, 58, 21, 28, 67, 55, 68, 37, 1, 69, 22, 29, 59,  &
61          &   78, 101, 79, 77, 70, 76, 47, 75, 30, 81, 60, 102, 80, 82, 23, 2, 83, 84, 85,  &
62          &   86, 103, 61, 31, 87, 56, 38, 71, 48, 88, 3, 62, 89, 24, 7, 49, 32, 104, 72, 90,  &
63          &   63, 39, 4, 8, 50, 91, 64, 40, 33, 25, 51, 95, 96, 73, 65, 9, 41, 97, 92, 105,  &
64          &   52, 5, 98, 10, 42, 99, 100, 66, 11, 74, 34, 53, 26, 6, 106, 12, 43, 13, 54, 93,  &
65          &   44, 107, 94, 14, 108, 15, 16, 109, 17, 18, 110, 111, 112 &
66          & /)
67   
68    real(jprb) :: hook_handle
69
70!#include "surdi.intfb.h"
71#include "surrtab.intfb.h"
72#include "surrtpk.intfb.h"
73#include "surrtrf.intfb.h"
74#include "rrtm_init_140gp.intfb.h"
75#include "srtm_init.intfb.h"
76
77    if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',0,hook_handle)
78
79    ! The IFS implementation of RRTMG uses many global variables.  In
80    ! the IFS these will have been set up already; otherwise set them
81    ! up now.
82    if (config%do_setup_ifsrrtm) then
83      !call SURDI
84      call SURRTAB
85      call SURRTPK
86      call SURRTRF
87      call RRTM_INIT_140GP(directory)
88      call SRTM_INIT(directory)
89    end if
90
91    ! Cloud and aerosol properties can only be defined per band
92    config%do_cloud_aerosol_per_sw_g_point = .false.
93    config%do_cloud_aerosol_per_lw_g_point = .false.
94
95    config%n_g_sw = jpgsw
96    config%n_g_lw = jpglw
97    config%n_bands_sw = 14
98    config%n_bands_lw = 16
99
100    ! Wavenumber ranges of each band may be needed so that the user
101    ! can compute UV and photosynthetically active radiation for a
102    ! particular wavelength range
103    call config%gas_optics_sw%spectral_def%allocate_bands_only( &
104         &  [2600.0_jprb, 3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, &
105         &   8050.0_jprb, 12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 820.0_jprb], &
106         &  [3250.0_jprb, 4000.0_jprb, 4650.0_jprb, 5150.0_jprb, 6150.0_jprb, 7700.0_jprb, 8050.0_jprb, &
107         &   12850.0_jprb, 16000.0_jprb, 22650.0_jprb, 29000.0_jprb, 38000.0_jprb, 50000.0_jprb, 2600.0_jprb])
108    call config%gas_optics_lw%spectral_def%allocate_bands_only( &
109         &  [10.0_jprb, 350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, &
110         &   1180.0_jprb, 1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb], &
111         &  [350.0_jprb, 500.0_jprb, 630.0_jprb, 700.0_jprb, 820.0_jprb, 980.0_jprb, 1080.0_jprb, 1180.0_jprb, &
112         &   1390.0_jprb, 1480.0_jprb, 1800.0_jprb, 2080.0_jprb, 2250.0_jprb, 2380.0_jprb, 2600.0_jprb, 3250.0_jprb])
113
114    allocate(config%i_band_from_g_sw          (config%n_g_sw))
115    allocate(config%i_band_from_g_lw          (config%n_g_lw))
116    allocate(config%i_band_from_reordered_g_sw(config%n_g_sw))
117    allocate(config%i_band_from_reordered_g_lw(config%n_g_lw))
118    allocate(config%i_g_from_reordered_g_sw(config%n_g_sw))
119    allocate(config%i_g_from_reordered_g_lw(config%n_g_lw))
120
121    ! Shortwave starts at 16: need to start at 1
122    config%i_band_from_g_sw = ngb_sw - ngb_sw(1)+1
123    config%i_band_from_g_lw = ngb_lw
124
125    if (config%i_solver_sw == ISolverSpartacus) then
126      ! SPARTACUS requires g points ordered in approximately
127      ! increasing order of optical depth
128      config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
129    else
130      ! Implied-do for no reordering
131!      config%i_g_from_reordered_g_sw = RRTM_GPOINT_REORDERING_SW
132      config%i_g_from_reordered_g_sw = (/ (irep, irep=1,config%n_g_sw) /)
133    end if
134
135    if (config%i_solver_lw == ISolverSpartacus) then
136      ! SPARTACUS requires g points ordered in approximately
137      ! increasing order of optical depth
138      config%i_g_from_reordered_g_lw = RRTM_GPOINT_REORDERING_LW
139    else
140      ! Implied-do for no reordering
141      config%i_g_from_reordered_g_lw = (/ (irep, irep=1,config%n_g_lw) /)
142    end if
143
144    config%i_band_from_reordered_g_sw &
145         = config%i_band_from_g_sw(config%i_g_from_reordered_g_sw)
146
147    config%i_band_from_reordered_g_lw &
148         = config%i_band_from_g_lw(config%i_g_from_reordered_g_lw)
149
150    ! The i_spec_* variables are used solely for storing spectral
151    ! data, and this can either be by band or by g-point
152    if (config%do_save_spectral_flux) then
153      if (config%do_save_gpoint_flux) then
154        config%n_spec_sw = config%n_g_sw
155        config%n_spec_lw = config%n_g_lw
156        config%i_spec_from_reordered_g_sw => config%i_g_from_reordered_g_sw
157        config%i_spec_from_reordered_g_lw => config%i_g_from_reordered_g_lw
158      else
159        config%n_spec_sw = config%n_bands_sw
160        config%n_spec_lw = config%n_bands_lw
161        config%i_spec_from_reordered_g_sw => config%i_band_from_reordered_g_sw
162        config%i_spec_from_reordered_g_lw => config%i_band_from_reordered_g_lw
163      end if
164    else
165      config%n_spec_sw = 0
166      config%n_spec_lw = 0
167      nullify(config%i_spec_from_reordered_g_sw)
168      nullify(config%i_spec_from_reordered_g_lw)
169    end if
170
171    if (lhook) call dr_hook('radiation_ifs_rrtm:setup_gas_optics',1,hook_handle)
172
173  end subroutine setup_gas_optics
174
175
176  !---------------------------------------------------------------------
177  ! Scale gas mixing ratios according to required units
178  subroutine set_gas_units(gas)
179
180    use radiation_gas,           only : gas_type, IMassMixingRatio
181    type(gas_type),    intent(inout) :: gas
182
183    call gas%set_units(IMassMixingRatio)
184
185  end subroutine set_gas_units
186
187
188  !---------------------------------------------------------------------
189  ! Compute gas optical depths, shortwave scattering, Planck function
190  ! and incoming shortwave radiation at top-of-atmosphere
191  subroutine gas_optics(ncol,nlev,istartcol,iendcol, &
192       &  config, single_level, thermodynamics, gas, &
193       &  od_lw, od_sw, ssa_sw, lw_albedo, planck_hl, lw_emission, &
194       &  incoming_sw)
195
196    use parkind1,                 only : jprb, jpim
197
198    USE PARRRTM  , ONLY : JPBAND, JPXSEC, JPINPX
199    USE YOERRTM  , ONLY : JPGPT_LW => JPGPT
200    USE YOESRTM  , ONLY : JPGPT_SW => JPGPT 
201    !USE YOMDIMV  , ONLY : YRDIMV
202    use yomhook  , only : lhook, dr_hook
203
204    use radiation_config,         only : config_type, ISolverSpartacus
205    use radiation_thermodynamics, only : thermodynamics_type
206    use radiation_single_level,   only : single_level_type
207    use radiation_gas
208
209    integer, intent(in) :: ncol               ! number of columns
210    integer, intent(in) :: nlev               ! number of model levels
211    integer, intent(in) :: istartcol, iendcol ! range of columns to process
212    type(config_type), intent(in) :: config
213    type(single_level_type),  intent(in) :: single_level
214    type(thermodynamics_type),intent(in) :: thermodynamics
215    type(gas_type),           intent(in) :: gas
216
217    ! Longwave albedo of the surface
218    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
219         &  intent(in), optional :: lw_albedo
220
221    ! Gaseous layer optical depth in longwave and shortwave, and
222    ! shortwave single scattering albedo (i.e. fraction of extinction
223    ! due to Rayleigh scattering) at each g-point
224    real(jprb), dimension(config%n_g_lw,nlev,istartcol:iendcol), intent(out) :: &
225         &   od_lw
226    real(jprb), dimension(config%n_g_sw,nlev,istartcol:iendcol), intent(out) :: &
227         &   od_sw, ssa_sw
228
229    ! The Planck function (emitted flux from a black body) at half
230    ! levels at each longwave g-point
231    real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), &
232         &   intent(out), optional :: planck_hl
233    ! Planck function for the surface (W m-2)
234    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
235         &   intent(out), optional :: lw_emission
236
237    ! The incoming shortwave flux into a plane perpendicular to the
238    ! incoming radiation at top-of-atmosphere in each of the shortwave
239    ! g-points
240    real(jprb), dimension(config%n_g_sw,istartcol:iendcol), &
241         &   intent(out), optional :: incoming_sw
242
243    real(jprb) :: incoming_sw_scale(istartcol:iendcol)
244
245    ! The variables in capitals are used in the same way as the
246    ! equivalent routine in the IFS
247
248    real(jprb) :: ZOD_LW(JPGPT_LW,nlev,istartcol:iendcol) ! Note ordering of dimensions
249    real(jprb) :: ZOD_SW(istartcol:iendcol,nlev,JPGPT_SW)
250    real(jprb) :: ZSSA_SW(istartcol:iendcol,nlev,JPGPT_SW)
251    real(jprb) :: ZINCSOL(istartcol:iendcol,JPGPT_SW)
252
253    real(jprb) :: ZCOLMOL(istartcol:iendcol,nlev)
254    real(jprb) :: ZCOLDRY(istartcol:iendcol,nlev)
255    real(jprb) :: ZWBRODL(istartcol:iendcol,nlev) !BROADENING GASES,column density (mol/cm2)
256    real(jprb) :: ZCOLBRD(istartcol:iendcol,nlev) !BROADENING GASES, column amount
257    real(jprb) :: ZWKL(istartcol:iendcol,JPINPX,nlev)
258
259    real(jprb) :: ZWX(istartcol:iendcol,JPXSEC,nlev) ! Amount of trace gases
260   
261    real(jprb) :: ZFLUXFAC, ZPI
262
263    ! - from AER
264    real(jprb) :: ZTAUAERL(istartcol:iendcol,nlev,JPBAND)
265
266    !- from INTFAC     
267    real(jprb) :: ZFAC00(istartcol:iendcol,nlev)
268    real(jprb) :: ZFAC01(istartcol:iendcol,nlev)
269    real(jprb) :: ZFAC10(istartcol:iendcol,nlev)
270    real(jprb) :: ZFAC11(istartcol:iendcol,nlev)
271   
272    !- from FOR
273    real(jprb) :: ZFORFAC(istartcol:iendcol,nlev)
274    real(jprb) :: ZFORFRAC(istartcol:iendcol,nlev)
275    integer    :: INDFOR(istartcol:iendcol,nlev)
276
277    !- from MINOR
278    integer    :: INDMINOR(istartcol:iendcol,nlev)
279    real(jprb) :: ZSCALEMINOR(istartcol:iendcol,nlev)
280    real(jprb) :: ZSCALEMINORN2(istartcol:iendcol,nlev)
281    real(jprb) :: ZMINORFRAC(istartcol:iendcol,nlev)
282   
283    real(jprb)     :: &                 
284         &  ZRAT_H2OCO2(istartcol:iendcol,nlev),ZRAT_H2OCO2_1(istartcol:iendcol,nlev), &
285         &  ZRAT_H2OO3(istartcol:iendcol,nlev) ,ZRAT_H2OO3_1(istartcol:iendcol,nlev), &
286         &  ZRAT_H2ON2O(istartcol:iendcol,nlev),ZRAT_H2ON2O_1(istartcol:iendcol,nlev), &
287         &  ZRAT_H2OCH4(istartcol:iendcol,nlev),ZRAT_H2OCH4_1(istartcol:iendcol,nlev), &
288         &  ZRAT_N2OCO2(istartcol:iendcol,nlev),ZRAT_N2OCO2_1(istartcol:iendcol,nlev), &
289         &  ZRAT_O3CO2(istartcol:iendcol,nlev) ,ZRAT_O3CO2_1(istartcol:iendcol,nlev)
290   
291    !- from INTIND
292    integer :: JP(istartcol:iendcol,nlev)
293    integer :: JT(istartcol:iendcol,nlev)
294    integer :: JT1(istartcol:iendcol,nlev)
295
296    !- from PRECISE             
297    real(jprb) :: ZONEMINUS, ZONEMINUS_ARRAY(istartcol:iendcol)
298
299    !- from PROFDATA             
300    real(jprb) :: ZCOLH2O(istartcol:iendcol,nlev)
301    real(jprb) :: ZCOLCO2(istartcol:iendcol,nlev)
302    real(jprb) :: ZCOLO3(istartcol:iendcol,nlev)
303    real(jprb) :: ZCOLN2O(istartcol:iendcol,nlev)
304    real(jprb) :: ZCOLCH4(istartcol:iendcol,nlev)
305    real(jprb) :: ZCOLO2(istartcol:iendcol,nlev)
306    real(jprb) :: ZCO2MULT(istartcol:iendcol,nlev)
307    integer    :: ILAYTROP(istartcol:iendcol)
308    integer    :: ILAYSWTCH(istartcol:iendcol)
309    integer    :: ILAYLOW(istartcol:iendcol)
310
311    !- from PROFILE             
312    real(jprb) :: ZPAVEL(istartcol:iendcol,nlev)
313    real(jprb) :: ZTAVEL(istartcol:iendcol,nlev)
314    real(jprb) :: ZPZ(istartcol:iendcol,0:nlev)
315    real(jprb) :: ZTZ(istartcol:iendcol,0:nlev)
316   
317    !- from SELF             
318    real(jprb) :: ZSELFFAC(istartcol:iendcol,nlev)
319    real(jprb) :: ZSELFFRAC(istartcol:iendcol,nlev)
320    integer :: INDSELF(istartcol:iendcol,nlev)
321
322    !- from SP             
323    real(jprb) :: ZPFRAC(istartcol:iendcol,JPGPT_LW,nlev)
324   
325    !- from SURFACE             
326    integer :: IREFLECT(istartcol:iendcol)
327
328    real(jprb) :: pressure_fl(ncol, nlev), temperature_fl(ncol, nlev)
329
330    ! If nlev is less than the number of heights at which gas mixing
331    ! ratios are stored, then we assume that the lower part of the
332    ! atmosphere is required. This enables nlev=1 to be passed in to
333    ! the routine, in which case the gas properties of the lowest
334    ! layer are provided, useful for canopy radiative transfer.
335    integer :: istartlev, iendlev
336
337    integer :: jlev, jgreorder, jg, ig, iband, jcol
338
339    real(jprb) :: hook_handle
340
341#include "rrtm_prepare_gases.intfb.h"
342#include "rrtm_setcoef_140gp.intfb.h"
343#include "rrtm_gas_optical_depth.intfb.h"
344#include "srtm_setcoef.intfb.h"
345#include "srtm_gas_optical_depth.intfb.h"
346
347    if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',0,hook_handle)
348
349    ! Compute start and end levels for indexing the gas mixing ratio
350    ! and thermodynamics arrays
351    iendlev   = ubound(gas%mixing_ratio,2)
352    istartlev = iendlev - nlev + 1
353
354    ZPI = 2.0_jprb*ASIN(1.0_jprb)
355    ZFLUXFAC = ZPI * 1.E+4
356    ZONEMINUS = 1.0_jprb - 1.0e-6_jprb
357    ZONEMINUS_ARRAY = ZONEMINUS
358
359!    if (.not. associated(YRDIMV)) then
360!      allocate(YRDIMV)
361!      YRDIMV%NFLEVG = nlev
362!    end if
363
364    do jlev=1,nlev
365      do jcol= istartcol,iendcol
366        pressure_fl(jcol,jlev) &
367            &  = 0.5_jprb * (thermodynamics%pressure_hl(jcol,jlev+istartlev-1) &
368            &               +thermodynamics%pressure_hl(jcol,jlev+istartlev))
369        temperature_fl(jcol,jlev) &
370            &  = 0.5_jprb * (thermodynamics%temperature_hl(jcol,jlev+istartlev-1) &
371            &               +thermodynamics%temperature_hl(jcol,jlev+istartlev))
372      end do
373    end do
374   
375    ! Check we have gas mixing ratios in the right units
376    call gas%assert_units(IMassMixingRatio)
377
378    ! Warning: O2 is hard-coded within the following function so the
379    ! user-provided concentrations of this gas are ignored for both
380    ! the longwave and shortwave
381    CALL RRTM_PREPARE_GASES &
382         & ( istartcol, iendcol, ncol, nlev, &
383         &   thermodynamics%pressure_hl(:,istartlev:iendlev+1), &
384         &   pressure_fl, &
385         &   thermodynamics%temperature_hl(:,istartlev:iendlev+1), &
386         &   temperature_fl, &
387         &   gas%mixing_ratio(:,istartlev:iendlev,IH2O), &
388         &   gas%mixing_ratio(:,istartlev:iendlev,ICO2), &
389         &   gas%mixing_ratio(:,istartlev:iendlev,ICH4), &
390         &   gas%mixing_ratio(:,istartlev:iendlev,IN2O), &
391         &   gas%mixing_ratio(:,istartlev:iendlev,INO2), &
392         &   gas%mixing_ratio(:,istartlev:iendlev,ICFC11), &
393         &   gas%mixing_ratio(:,istartlev:iendlev,ICFC12), &
394         &   gas%mixing_ratio(:,istartlev:iendlev,IHCFC22), &
395         &   gas%mixing_ratio(:,istartlev:iendlev,ICCl4), &
396         &   gas%mixing_ratio(:,istartlev:iendlev,IO3), &
397         &  ZCOLDRY, ZWBRODL,ZWKL, ZWX, &
398         &  ZPAVEL , ZTAVEL , ZPZ , ZTZ, IREFLECT) 
399
400    CALL RRTM_SETCOEF_140GP &
401         &( istartcol, iendcol, nlev , ZCOLDRY  , ZWBRODL , ZWKL , &
402         &  ZFAC00 , ZFAC01   , ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, JP, JT, JT1 , &
403         &  ZCOLH2O, ZCOLCO2  , ZCOLO3 , ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT , ZCOLBRD, &
404         &  ILAYTROP,ILAYSWTCH, ILAYLOW, ZPAVEL , ZTAVEL , ZSELFFAC, ZSELFFRAC, INDSELF, &
405         &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
406         &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
407         &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
408         &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)   
409
410    ZTAUAERL(istartcol:iendcol,:,:) = 0.0_jprb
411
412    CALL RRTM_GAS_OPTICAL_DEPTH &
413         &( istartcol, iendcol, nlev, ZOD_LW, ZPAVEL, ZCOLDRY, ZCOLBRD, ZWX ,&
414         &  ZTAUAERL, ZFAC00 , ZFAC01, ZFAC10 , ZFAC11 , ZFORFAC,ZFORFRAC,INDFOR, &
415         &  JP, JT, JT1, ZONEMINUS ,&
416         &  ZCOLH2O , ZCOLCO2, ZCOLO3, ZCOLN2O, ZCOLCH4, ZCOLO2,ZCO2MULT ,&
417         &  ILAYTROP, ILAYSWTCH,ILAYLOW, ZSELFFAC, ZSELFFRAC, INDSELF, ZPFRAC, &
418         &  INDMINOR,ZSCALEMINOR,ZSCALEMINORN2,ZMINORFRAC,&
419         &  ZRAT_H2OCO2, ZRAT_H2OCO2_1, ZRAT_H2OO3, ZRAT_H2OO3_1, &
420         &  ZRAT_H2ON2O, ZRAT_H2ON2O_1, ZRAT_H2OCH4, ZRAT_H2OCH4_1, &
421         &  ZRAT_N2OCO2, ZRAT_N2OCO2_1, ZRAT_O3CO2, ZRAT_O3CO2_1)     
422
423    if (present(lw_albedo)) then
424   
425      call planck_function_atmos(nlev, istartcol, iendcol, config, &
426           &                     thermodynamics, ZPFRAC, planck_hl)
427
428      if (single_level%is_simple_surface) then
429        call planck_function_surf(istartcol, iendcol, config, &
430             &                    single_level%skin_temperature, ZPFRAC(:,:,1), &
431             &                    lw_emission)
432       
433        ! The following can be used to extract the parameters defined at
434        ! the top of the planck_function routine below:
435        !write(*,'(a,140(e12.5,","),a)') 'ZPFRAC_surf=[', &
436        !&  sum(ZPFRAC(istartcol:iendcol,:,1),1) / (iendcol+1-istartcol), ']'
437       
438        ! lw_emission at this point is actually the planck function of
439        ! the surface
440        lw_emission = lw_emission * (1.0_jprb - lw_albedo)
441      else
442        ! Longwave emission has already been computed
443        if (config%use_canopy_full_spectrum_lw) then
444          lw_emission = transpose(single_level%lw_emission(istartcol:iendcol,:))
445        else
446          lw_emission = transpose(single_level%lw_emission(istartcol:iendcol, &
447               & config%i_emiss_from_band_lw(config%i_band_from_reordered_g_lw)))
448        end if
449      end if
450
451    end if
452
453    if (config%i_solver_lw == ISolverSpartacus) then
454      !    if (.true.) then
455      ! We need to rearrange the gas optics info in memory: reordering
456      ! the g points in order of approximately increasing optical
457      ! depth (for efficient 3D processing on only the regions of the
458      ! spectrum that are optically thin for gases) and reorder in
459      ! pressure since the the functions above treat pressure
460      ! decreasing with increasing index.  Note that the output gas
461      ! arrays have dimensions in a different order to the inputs,
462      ! so there is some inefficiency here.
463      do jgreorder = 1,config%n_g_lw
464        iband = config%i_band_from_reordered_g_lw(jgreorder)
465        ig = config%i_g_from_reordered_g_lw(jgreorder)
466       
467        ! Top-of-atmosphere half level
468        do jlev = 1,nlev
469          do jcol = istartcol,iendcol
470            ! Some g points can return negative optical depths;
471            ! specifically original g points 54-56 which causes
472            ! unphysical single-scattering albedo when combined with
473            ! aerosol
474            od_lw(jgreorder,jlev,jcol) &
475                 &   = max(config%min_gas_od_lw, ZOD_LW(ig,nlev+1-jlev,jcol))
476          end do
477        end do
478      end do
479    else
480      ! G points have not been reordered
481      do jcol = istartcol,iendcol
482        do jlev = 1,nlev
483          ! Check for negative optical depth
484          od_lw(:,jlev,jcol) = max(config%min_gas_od_lw, ZOD_LW(:,nlev+1-jlev,jcol))
485        end do
486      end do
487    end if
488   
489    CALL SRTM_SETCOEF &
490         & ( istartcol, iendcol, nlev,&
491         & ZPAVEL  , ZTAVEL,&
492         & ZCOLDRY , ZWKL,&
493         & ILAYTROP,&
494         & ZCOLCH4  , ZCOLCO2 , ZCOLH2O , ZCOLMOL  , ZCOLO2 , ZCOLO3,&
495         & ZFORFAC , ZFORFRAC , INDFOR  , ZSELFFAC, ZSELFFRAC, INDSELF, &
496         & ZFAC00  , ZFAC01   , ZFAC10  , ZFAC11,&
497         & JP      , JT       , JT1     , single_level%cos_sza(istartcol:iendcol)  &
498         & ) 
499   
500    ! SRTM_GAS_OPTICAL_DEPTH will not initialize profiles when the sun
501    ! is below the horizon, so we do it here
502    ZOD_SW(istartcol:iendcol,:,:)  = 0.0_jprb
503    ZSSA_SW(istartcol:iendcol,:,:) = 0.0_jprb
504    ZINCSOL(istartcol:iendcol,:)   = 0.0_jprb
505
506    CALL SRTM_GAS_OPTICAL_DEPTH &
507         &( istartcol, iendcol , nlev  , ZONEMINUS_ARRAY,&
508         & single_level%cos_sza(istartcol:iendcol), ILAYTROP,&
509         & ZCOLCH4 , ZCOLCO2  , ZCOLH2O, ZCOLMOL , ZCOLO2   , ZCOLO3,&
510         & ZFORFAC , ZFORFRAC , INDFOR , ZSELFFAC, ZSELFFRAC, INDSELF,&
511         & ZFAC00  , ZFAC01   , ZFAC10 , ZFAC11  ,&
512         & JP      , JT       , JT1    ,&
513         & ZOD_SW  , ZSSA_SW  , ZINCSOL )
514   
515    ! Scale the incoming solar per band, if requested
516    if (config%use_spectral_solar_scaling) then
517      do jg = 1,JPGPT_SW
518        do jcol = istartcol,iendcol
519          ZINCSOL(jcol,jg) = ZINCSOL(jcol,jg) * &
520            &   single_level%spectral_solar_scaling(config%i_band_from_reordered_g_sw(jg))
521        end do
522      end do
523    end if
524
525    ! Scaling factor to ensure that the total solar irradiance is as
526    ! requested.  Note that if the sun is below the horizon then
527    ! ZINCSOL will be zero.
528    if (present(incoming_sw)) then
529      do jcol = istartcol,iendcol
530        if (single_level%cos_sza(jcol) > 0.0_jprb) then
531! Added for DWD (2020)
532!NEC$ nounroll
533          incoming_sw_scale(jcol) = single_level%solar_irradiance / sum(ZINCSOL(jcol,:))
534        else
535          incoming_sw_scale(jcol) = 1.0_jprb
536        end if
537      end do
538    end if
539
540    if (config%i_solver_sw == ISolverSpartacus) then
541!    if (.true.) then
542      ! Account for reordered g points
543      do jgreorder = 1,config%n_g_sw
544        ig = config%i_g_from_reordered_g_sw(jgreorder)
545        do jlev = 1,nlev
546          do jcol = istartcol,iendcol
547            ! Check for negative optical depth
548            od_sw (jgreorder,nlev+1-jlev,jcol) &
549                 &  = max(config%min_gas_od_sw, ZOD_SW (jcol,jlev,ig))
550            ssa_sw(jgreorder,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,ig)
551           end do
552        end do
553        if (present(incoming_sw)) then
554          incoming_sw(jgreorder,:) &
555               &  = incoming_sw_scale(:) * ZINCSOL(:,ig)
556        end if
557      end do
558    else
559      ! G points have not been reordered
560      do jcol = istartcol,iendcol
561        do jlev = 1,nlev
562          do jg = 1,config%n_g_sw
563            ! Check for negative optical depth
564            od_sw (jg,nlev+1-jlev,jcol) = max(config%min_gas_od_sw, ZOD_SW(jcol,jlev,jg))
565            ssa_sw(jg,nlev+1-jlev,jcol) = ZSSA_SW(jcol,jlev,jg)
566          end do
567        end do
568        if (present(incoming_sw)) then
569          do jg = 1,config%n_g_sw
570            incoming_sw(jg,jcol) = incoming_sw_scale(jcol) * ZINCSOL(jcol,jg)
571          end do
572        end if
573      end do
574    end if
575   
576    if (lhook) call dr_hook('radiation_ifs_rrtm:gas_optics',1,hook_handle)
577   
578  end subroutine gas_optics
579 
580
581  !---------------------------------------------------------------------
582  ! Compute Planck function of the atmosphere
583  subroutine planck_function_atmos(nlev,istartcol,iendcol, &
584       config, thermodynamics, PFRAC, &
585       planck_hl)
586
587    use parkind1,                 only : jprb, jpim
588
589    USE YOERRTM  , ONLY : JPGPT_LW => JPGPT
590    use yoerrtwn, only : totplnk, delwave
591
592    use yomhook, only : lhook, dr_hook
593
594    use radiation_config,         only : config_type, ISolverSpartacus
595    use radiation_thermodynamics, only : thermodynamics_type
596    !use radiation_gas
597
598    integer, intent(in) :: nlev               ! number of model levels
599    integer, intent(in) :: istartcol, iendcol ! range of columns to process
600    type(config_type), intent(in) :: config
601    type(thermodynamics_type),intent(in) :: thermodynamics
602    real(jprb), intent(in) :: PFRAC(istartcol:iendcol,JPGPT_LW,nlev)
603
604    ! The Planck function (emitted flux from a black body) at half
605    ! levels at each longwave g-point
606    real(jprb), dimension(config%n_g_lw,nlev+1,istartcol:iendcol), intent(out) :: &
607         &   planck_hl
608
609    ! Planck function values per band
610    real(jprb), dimension(istartcol:iendcol, config%n_bands_lw) :: planck_store
611
612    ! Look-up table variables for Planck function
613    real(jprb), dimension(istartcol:iendcol) :: frac
614    integer,    dimension(istartcol:iendcol) :: ind
615
616    ! Temperature (K) of a half-level
617    real(jprb) :: temperature
618
619    real(jprb) :: factor, planck_tmp(istartcol:iendcol,config%n_g_lw)
620    real(jprb) :: ZFLUXFAC
621
622    integer :: jlev, jgreorder, jg, ig, iband, jband, jcol, ilevoffset
623
624    real(jprb) :: hook_handle
625
626    if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_atmos',0,hook_handle)
627
628    ZFLUXFAC = 2.0_jprb*ASIN(1.0_jprb) * 1.0e4_jprb
629   
630    ! nlev may be less than the number of original levels, in which
631    ! case we assume that the user wants the lower part of the
632    ! atmosphere
633    ilevoffset = ubound(thermodynamics%temperature_hl,2)-nlev-1
634
635    ! Work out interpolations: for each half level, the index of the
636    ! lowest interpolation bound, and the fraction into interpolation
637    ! interval
638    do jlev = 1,nlev+1
639      do jcol = istartcol,iendcol
640        temperature = thermodynamics%temperature_hl(jcol,jlev+ilevoffset)
641        if (temperature < 339.0_jprb .and. temperature >= 160.0_jprb) then
642          ! Linear interpolation between -113 and 66 degC
643          ind(jcol)  = int(temperature - 159.0_jprb)
644          frac(jcol) = temperature - int(temperature)
645        else if(temperature >= 339.0_jprb) then
646          ! Extrapolation above 66 degC
647          ind(jcol)  = 180
648          frac(jcol) = temperature - 339.0_jprb
649        else
650          ! Cap below -113 degC (to avoid possible negative Planck
651          ! function values)
652          ind(jcol)  = 1
653          frac(jcol) = 0.0_jprb
654        end if
655      end do
656
657      ! Calculate Planck functions per band
658      do jband = 1,config%n_bands_lw
659        factor = zfluxfac * delwave(jband)
660        do jcol = istartcol,iendcol
661          planck_store(jcol,jband) = factor &
662               &  * (totplnk(ind(jcol),jband) &
663               &  + frac(jcol)*(totplnk(ind(jcol)+1,jband)-totplnk(ind(jcol),jband)))
664        end do
665      end do
666
667      if (config%i_solver_lw == ISolverSpartacus) then
668        ! We need to rearrange the gas optics info in memory:
669        ! reordering the g points in order of approximately increasing
670        ! optical depth (for efficient 3D processing on only the
671        ! regions of the spectrum that are optically thin for gases)
672        ! and reorder in pressure since the the functions above treat
673        ! pressure decreasing with increasing index.
674        if (jlev == 1) then
675          ! Top-of-atmosphere half level - note that PFRAC is on model
676          ! levels not half levels
677          do jgreorder = 1,config%n_g_lw
678            iband = config%i_band_from_reordered_g_lw(jgreorder)
679            ig = config%i_g_from_reordered_g_lw(jgreorder)
680            planck_hl(jgreorder,1,:) = planck_store(:,iband) &
681                 &   * PFRAC(:,ig,nlev)
682          end do
683        else
684          do jgreorder = 1,config%n_g_lw
685            iband = config%i_band_from_reordered_g_lw(jgreorder)
686            ig = config%i_g_from_reordered_g_lw(jgreorder)
687            planck_hl(jgreorder,jlev,:) &
688                   &   = planck_store(:,iband) &
689                   &   * PFRAC(:,ig,nlev+2-jlev)
690          end do
691        end if
692      else
693        ! G points have not been reordered
694        if (jlev == 1) then
695          ! Top-of-atmosphere half level - note that PFRAC is on model
696          ! levels not half levels
697          do jg = 1,config%n_g_lw
698            iband = config%i_band_from_g_lw(jg)
699            planck_hl(jg,1,:) = planck_store(:,iband) * PFRAC(:,jg,nlev)
700          end do
701        else
702          do jg = 1,config%n_g_lw
703            iband = config%i_band_from_g_lw(jg)
704            planck_tmp(:,jg) = planck_store(:,iband) * PFRAC(:,jg,nlev+2-jlev)
705          end do
706          do jcol = istartcol,iendcol
707            planck_hl(:,jlev,jcol) = planck_tmp(jcol,:)
708          end do
709        end if
710      end if
711    end do
712
713    if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_atmos',1,hook_handle)
714
715  end subroutine planck_function_atmos
716
717
718  !---------------------------------------------------------------------
719  ! Compute Planck function of the surface
720  subroutine planck_function_surf(istartcol, iendcol, config, temperature, PFRAC, &
721       &  planck_surf)
722
723    use parkind1,                 only : jprb, jpim
724
725    USE YOERRTM  , ONLY : JPGPT_LW => JPGPT
726    use yoerrtwn, only : totplnk, delwave
727
728    use yomhook, only : lhook, dr_hook
729
730    use radiation_config,         only : config_type, ISolverSpartacus
731    !    use radiation_gas
732
733    integer, intent(in) :: istartcol, iendcol ! range of columns to process
734    type(config_type), intent(in) :: config
735    real(jprb), intent(in) :: temperature(:)
736
737    real(jprb), intent(in) :: PFRAC(istartcol:iendcol,JPGPT_LW)
738
739    ! Planck function of the surface (W m-2)
740    real(jprb), dimension(config%n_g_lw,istartcol:iendcol), &
741         &  intent(out) :: planck_surf
742
743    ! Planck function values per band
744    real(jprb), dimension(istartcol:iendcol, config%n_bands_lw) :: planck_store
745
746    ! Look-up table variables for Planck function
747    real(jprb), dimension(istartcol:iendcol) :: frac
748    integer,    dimension(istartcol:iendcol) :: ind
749
750    ! Temperature (K)
751    real(jprb) :: Tsurf
752
753    real(jprb) :: factor
754    real(jprb) :: ZFLUXFAC
755
756    integer :: jgreorder, jg, ig, iband, jband, jcol
757
758    real(jprb) :: hook_handle
759
760    if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_surf',0,hook_handle)
761
762    ZFLUXFAC = 2.0_jprb*ASIN(1.0_jprb) * 1.0e4_jprb
763
764    ! Work out surface interpolations
765    do jcol = istartcol,iendcol
766      Tsurf = temperature(jcol)
767      if (Tsurf < 339.0_jprb .and. Tsurf >= 160.0_jprb) then
768        ! Linear interpolation between -113 and 66 degC
769        ind(jcol)  = int(Tsurf - 159.0_jprb)
770        frac(jcol) = Tsurf - int(Tsurf)
771      else if(Tsurf >= 339.0_jprb) then
772        ! Extrapolation above 66 degC
773        ind(jcol)  = 180
774        frac(jcol) = Tsurf - 339.0_jprb
775      else
776        ! Cap below -113 degC (to avoid possible negative Planck
777        ! function values)
778        ind(jcol)  = 1
779        frac(jcol) = 0.0_jprb
780      end if
781    end do
782
783    ! Calculate Planck functions per band
784    do jband = 1,config%n_bands_lw
785      factor = zfluxfac * delwave(jband)
786      do jcol = istartcol,iendcol
787        planck_store(jcol,jband) = factor &
788             &  * (totplnk(ind(jcol),jband) &
789             &  + frac(jcol)*(totplnk(ind(jcol)+1,jband)-totplnk(ind(jcol),jband)))
790      end do
791    end do
792
793    if (config%i_solver_lw == ISolverSpartacus) then
794      ! We need to rearrange the gas optics info in memory: reordering
795      ! the g points in order of approximately increasing optical
796      ! depth (for efficient 3D processing on only the regions of
797      ! the spectrum that are optically thin for gases) and reorder
798      ! in pressure since the the functions above treat pressure
799      ! decreasing with increasing index.
800      do jgreorder = 1,config%n_g_lw
801        iband = config%i_band_from_reordered_g_lw(jgreorder)
802        ig = config%i_g_from_reordered_g_lw(jgreorder)
803        planck_surf(jgreorder,:) = planck_store(:,iband) * PFRAC(:,ig)
804      end do
805    else
806      ! G points have not been reordered
807      do jg = 1,config%n_g_lw
808        iband = config%i_band_from_g_lw(jg)
809        planck_surf(jg,:) = planck_store(:,iband) * PFRAC(:,jg)
810      end do
811    end if
812
813    if (lhook) call dr_hook('radiation_ifs_rrtm:planck_function_surf',1,hook_handle)
814   
815  end subroutine planck_function_surf
816
817
818  !---------------------------------------------------------------------
819  ! Externally facing function for computing the Planck function
820  ! without reference to any gas profile; typically this would be used
821  ! for computing the emission by facets of a complex surface.  Note
822  ! that this uses fixed "PFRAC" values, obtained by averaging over
823  ! those derived from RRTM-G for near-surface conditions over a line
824  ! of meridian from the ECMWF model.
825  subroutine planck_function(config, temperature, planck_surf)
826
827    use parkind1,                 only : jprb, jpim
828
829    use radiation_config,         only : config_type
830
831    type(config_type), intent(in) :: config
832    real(jprb), intent(in) :: temperature
833
834    ! Planck function of the surface (W m-2)
835    real(jprb), dimension(config%n_g_lw), &
836         &  intent(out) :: planck_surf
837
838    ! Fraction of each band contributed by each g-point within
839    ! it. Since there are 16 bands, this array sums to 16
840    real(jprb), parameter, dimension(1,140) :: frac &
841         = reshape( (/ 0.21227E+00, 0.18897E+00, 0.25491E+00, 0.17864E+00, 0.11735E+00, 0.38298E-01, 0.57871E-02, &
842         &    0.31753E-02, 0.53169E-03, 0.76476E-04, 0.16388E+00, 0.15241E+00, 0.14290E+00, 0.12864E+00, &
843         &    0.11615E+00, 0.10047E+00, 0.80013E-01, 0.60445E-01, 0.44918E-01, 0.63395E-02, 0.32942E-02, &
844         &    0.54541E-03, 0.15380E+00, 0.15194E+00, 0.14339E+00, 0.13138E+00, 0.11701E+00, 0.10081E+00, &
845         &    0.82296E-01, 0.61735E-01, 0.41918E-01, 0.45918E-02, 0.37743E-02, 0.30121E-02, 0.22500E-02, &
846         &    0.14490E-02, 0.55410E-03, 0.78364E-04, 0.15938E+00, 0.15146E+00, 0.14213E+00, 0.13079E+00, &
847         &    0.11672E+00, 0.10053E+00, 0.81566E-01, 0.61126E-01, 0.41150E-01, 0.44488E-02, 0.36950E-02, &
848         &    0.29101E-02, 0.21357E-02, 0.19609E-02, 0.14134E+00, 0.14390E+00, 0.13913E+00, 0.13246E+00, &
849         &    0.12185E+00, 0.10596E+00, 0.87518E-01, 0.66164E-01, 0.44862E-01, 0.49402E-02, 0.40857E-02, &
850         &    0.32288E-02, 0.23613E-02, 0.15406E-02, 0.58258E-03, 0.82171E-04, 0.29127E+00, 0.28252E+00, &
851         &    0.22590E+00, 0.14314E+00, 0.45494E-01, 0.71792E-02, 0.38483E-02, 0.65712E-03, 0.29810E+00, &
852         &    0.27559E+00, 0.11997E+00, 0.10351E+00, 0.84515E-01, 0.62253E-01, 0.41050E-01, 0.44217E-02, &
853         &    0.36946E-02, 0.29113E-02, 0.34290E-02, 0.55993E-03, 0.31441E+00, 0.27586E+00, 0.21297E+00, &
854         &    0.14064E+00, 0.45588E-01, 0.65665E-02, 0.34232E-02, 0.53199E-03, 0.19811E+00, 0.16833E+00, &
855         &    0.13536E+00, 0.11549E+00, 0.10649E+00, 0.93264E-01, 0.75720E-01, 0.56405E-01, 0.41865E-01, &
856         &    0.59331E-02, 0.26510E-02, 0.40040E-03, 0.32328E+00, 0.26636E+00, 0.21397E+00, 0.14038E+00, &
857         &    0.52142E-01, 0.38852E-02, 0.14601E+00, 0.13824E+00, 0.27703E+00, 0.22388E+00, 0.15446E+00, &
858         &    0.48687E-01, 0.98054E-02, 0.18870E-02, 0.11961E+00, 0.12106E+00, 0.13215E+00, 0.13516E+00, &
859         &    0.25249E+00, 0.16542E+00, 0.68157E-01, 0.59725E-02, 0.49258E+00, 0.33651E+00, 0.16182E+00, &
860         &    0.90984E-02, 0.95202E+00, 0.47978E-01, 0.91716E+00, 0.82857E-01, 0.77464E+00, 0.22536E+00 /), (/ 1,140 /) )
861
862    call planck_function_surf(1, 1, config, spread(temperature,1,1), &
863         &                    frac, planck_surf)
864
865  end subroutine planck_function
866
867end module radiation_ifs_rrtm
Note: See TracBrowser for help on using the repository browser.