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

Last change on this file since 4586 was 3908, checked in by idelkadi, 3 years ago

Online implementation of the radiative transfer code ECRAD in the LMDZ model.

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