SCALE-RM
scale_atmos_phy_rd_mstrnx.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 #include "scalelib.h"
15 module scale_atmos_phy_rd_mstrnx
16  !-----------------------------------------------------------------------------
17  !
18  !++ used modules
19  !
20  use scale_precision
21  use scale_io
22  use scale_prof
23  use scale_cpl_sfc_index
24 
25  use scale_atmos_phy_rd_common, only: &
26  i_sw, &
27  i_lw, &
28  i_dn, &
29  i_up
30 
31  use scale_atmos_hydrometeor, only: &
32  n_hyd
33  use scale_atmos_aerosol, only: &
34  n_ae
35  !-----------------------------------------------------------------------------
36  implicit none
37  private
38  !-----------------------------------------------------------------------------
39  !
40  !++ Public procedure
41  !
42  public :: atmos_phy_rd_mstrnx_setup
43  public :: atmos_phy_rd_mstrnx_flux
44 
45  !-----------------------------------------------------------------------------
46  !
47  !++ Public parameters & variables
48  !
49  !-----------------------------------------------------------------------------
50  !
51  !++ Private procedure
52  !
53  private :: rd_mstrn_setup
54  private :: rd_mstrn_dtrn3
55  private :: rd_mstrn_two_stream
56 
57  !-----------------------------------------------------------------------------
58  !
59  !++ Private parameters & variables
60  !
61  real(RP), private, parameter :: RD_cosSZA_min = 0.017_rp ! minimum SZA (>89.0)
62  real(RP), private, parameter :: RD_EPS = 1.e-4_rp
63 
64  real(RP), private :: RD_TOA = 100.0_rp
65  integer, private :: RD_KADD = 10
66  integer, private, parameter :: RD_naero = n_hyd + n_ae ! # of cloud/aerosol species
67  integer, private, parameter :: RD_hydro_str = 1 ! start index for cloud
68  integer, private, parameter :: RD_hydro_end = n_hyd ! end index for cloud
69  integer, private, parameter :: RD_aero_str = n_hyd + 1 ! start index for aerosol
70  integer, private, parameter :: RD_aero_end = n_hyd + n_ae ! end index for aerosol
71 
72  integer, private :: RD_KMAX ! # of computational cells: z for radiation scheme
73 
74  real(RP), private, allocatable :: RD_zh (:) ! altitude at the interface [km]
75  real(RP), private, allocatable :: RD_z (:) ! altitude at the center [km]
76  real(RP), private, allocatable :: RD_rhodz (:) ! density * delta z [kg/m2]
77  real(RP), private, allocatable :: RD_pres (:) ! pressure at the center [hPa]
78  real(RP), private, allocatable :: RD_presh (:) ! pressure at the interface [hPa]
79  real(RP), private, allocatable :: RD_temp (:) ! temperature at the center [K]
80  real(RP), private, allocatable :: RD_temph (:) ! temperature at the interface [K]
81  real(RP), private, allocatable :: RD_gas (:,:) ! gas species volume mixing ratio [ppmv]
82  real(RP), private, allocatable :: RD_cfc (:,:) ! CFCs volume mixing ratio [ppmv]
83  real(RP), private, allocatable :: RD_aerosol_conc(:,:) ! cloud/aerosol volume mixing ratio [ppmv]
84  real(RP), private, allocatable :: RD_aerosol_radi(:,:) ! cloud/aerosol effective radius [cm]
85  real(RP), private, allocatable :: RD_cldfrac (:) ! cloud fraction (0-1)
86 
87  integer, private :: I_MPAE2RD(RD_naero) ! look-up table between input aerosol category and MSTRN particle type
88 
89  character(len=H_LONG), private :: MSTRN_GASPARA_INPUTFILE = 'PARAG.29'
90  character(len=H_LONG), private :: MSTRN_AEROPARA_INPUTFILE = 'PARAPC.29'
91  character(len=H_LONG), private :: MSTRN_HYGROPARA_INPUTFILE = 'VARDATA.RM29'
92 
93  integer, private :: MSTRN_nband = 29
94 
95  integer, private, parameter :: MSTRN_nstream = 1
96  integer, private, parameter :: MSTRN_ch_limit = 10
97  integer, private, parameter :: MSTRN_nflag = 7 ! # of optical properties flag
98  integer, private, parameter :: MSTRN_nfitP = 26 ! # of fitting point for log(pressure)
99  integer, private, parameter :: MSTRN_nfitT = 3 ! # of fitting point for temperature
100 
101  integer, private, parameter :: MSTRN_ngas = 7
102  ! 1: H2O
103  ! 2: CO2
104  ! 3: O3
105  ! 4: N2O
106  ! 5: CO
107  ! 6: CH4
108  ! 7: O2
109  integer, private, parameter :: MSTRN_ncfc = 28
110  ! 1: CFC-11
111  ! 2: CFC-12
112  ! 3: CFC-13
113  ! 4: CFC-14
114  ! 5: CFC-113
115  ! 6: CFC-114
116  ! 7: CFC-115
117  ! 8: HCFC-21
118  ! 9: HCFC-22
119  ! 10: HCFC-123
120  ! 11: HCFC-124
121  ! 12: HCFC-141b
122  ! 13: HCFC-142b
123  ! 14: HCFC-225ca
124  ! 15: HCFC-225cb
125  ! 16: HFC-32
126  ! 17: HFC-125
127  ! 18: HFC-134
128  ! 19: HFC-134a
129  ! 20: HFC-143a
130  ! 21: HFC-152a
131  ! 22: SF6
132  ! 23: ClONO2
133  ! 24: CCl4
134  ! 25: N2O5
135  ! 26: C2F6
136  ! 27: HNO4
137  integer, private :: MSTRN_nptype = 9
138  ! 1: water cloud
139  ! 2: ice cloud
140  ! 3: Soil dust
141  ! 4: Carbonacerous (BC/OC=0.3)
142  ! 5: Carbonacerous (BC/OC=0.15)
143  ! 6: Carbonacerous (BC/OC=0.)
144  ! 7: Black carbon
145  ! 8: Sulfate
146  ! 9: Sea salt
147 
148  integer, private, parameter :: MSTRN_nsfc = 7
149  ! 1: ocean
150  ! 2: wet land
151  ! 3: dry land
152  ! 4: low plants
153  ! 5: forest
154  ! 6: snow
155  ! 7: ice
156  integer, private, parameter :: MSTRN_nfitPLK = 5
157  integer, private, parameter :: MSTRN_nplkord = 3
158  integer, private, parameter :: MSTRN_nmoment = 6
159  integer, private :: MSTRN_nradius = 8
160  integer, private, parameter :: MSTRN_ncloud = 2
161 
162  logical, private :: ATMOS_PHY_RD_MSTRN_ONLY_QCI = .false.
163  logical, private :: ATMOS_PHY_RD_MSTRN_ONLY_TROPOCLOUD = .false.
164  logical, private :: ATMOS_PHY_RD_MSTRN_USE_AERO = .false.
165 
166 
167 
168  real(RP), private, allocatable :: waveh (:) ! wavenumbers at band boundary [1/cm]
169 
170  real(RP), private, allocatable :: logfitP (:) ! fitting point for log10(pressure)
171  real(RP), private, allocatable :: fitT (:) ! fitting point for temperature
172  real(RP), private, allocatable :: logfitT (:) ! fitting point for log10(temperature)
173  integer, private, allocatable :: iflgb (:,:) ! optical properties flag in each band
174  integer, private, allocatable :: nch (:) ! number of subintervals in each band
175  real(RP), private, allocatable :: wgtch (:,:) ! weights of subintervals in each band
176  integer, private, allocatable :: ngasabs (:) ! number of absorbers(gas) in each band
177  integer, private, allocatable :: igasabs (:,:) ! index of absorbers(gas) in each band
178 
179  real(RP), private, allocatable :: akd (:,:,:,:,:) ! absorption coefficient table
180  real(RP), private, allocatable :: skd (:,:,:,:) ! absorption coefficient table for H2O self broadening
181  real(RP), private, allocatable :: acfc_pow(:,:) ! 10 to the power of absorption coefficient table for CFC
182 
183  real(RP), private, allocatable :: fitPLK (:,:) ! fitting point for planck function
184  real(RP), private, allocatable :: fsol (:) ! solar insolation in each band
185  real(RP), private :: fsol_tot ! total solar insolation
186  real(RP), private, allocatable :: sfc (:,:) ! surface condition in each band
187  real(RP), private, allocatable :: rayleigh(:) ! rayleigh scattering in each band
188  real(RP), private, allocatable :: qmol (:,:) ! moments for rayleigh scattering phase function
189  real(RP), private, allocatable :: q (:,:,:,:) ! moments for aerosol scattering phase function
190 
191  integer, private, allocatable :: hygro_flag(:) ! flag for hygroscopic enlargement
192  real(RP), private, allocatable :: radmode (:,:) ! radius mode for hygroscopic parameter
193 
194  integer, private, allocatable :: ptype_nradius(:) ! number limit of radius mode for cloud / aerosol particle type
195 
196 
197  ! index for optical flag iflgb
198  integer, private, parameter :: I_SWLW = 4
199  integer, private, parameter :: I_H2O_continuum = 5
200  integer, private, parameter :: I_CFC_continuum = 7
201  ! for cloud type
202  integer, private, parameter :: I_ClearSky = 1
203  integer, private, parameter :: I_Cloud = 2
204 
205  ! pre-calc
206  real(RP), private :: RHO_std ! rho(0C,1atm) [kg/m3]
207 
208  real(RP), private :: M(2) ! discrete quadrature mu for two-stream approximation
209  real(RP), private :: W(2) ! discrete quadrature w for two-stream approximation
210  real(RP), private :: Wmns(2), Wpls(2) ! W-, W+
211  real(RP), private :: Wbar(2), Wscale(2)
212  !$acc declare create(m, w, wmns, wpls, wscale)
213 
214  !-----------------------------------------------------------------------------
215 contains
216  !-----------------------------------------------------------------------------
218  subroutine atmos_phy_rd_mstrnx_setup( &
219  KA, KS, KE, &
220  CZ, FZ )
221  use scale_prc, only: &
222  prc_abort
223  use scale_time, only: &
224  time_nowdate
225  use scale_atmos_grid_cartesc_real, only: &
226  atmos_grid_cartesc_real_basepoint_lat
227  use scale_atmos_phy_rd_profile, only: &
228  rd_profile_setup => atmos_phy_rd_profile_setup, &
229  rd_profile_setup_zgrid => atmos_phy_rd_profile_setup_zgrid, &
230  rd_profile_read => atmos_phy_rd_profile_read
231  use scale_atmos_hydrometeor, only: &
232  n_hyd
233  use scale_atmos_aerosol, only: &
234  n_ae
235  implicit none
236  integer, intent(in) :: ka, ks, ke
237 
238  real(rp), intent(in) :: cz(ka)
239  real(rp), intent(in) :: fz(0:ka)
240 
241  real(rp) :: atmos_phy_rd_mstrn_toa
242  integer :: atmos_phy_rd_mstrn_kadd
243  character(len=H_LONG) :: atmos_phy_rd_mstrn_gaspara_in_filename
244  character(len=H_LONG) :: atmos_phy_rd_mstrn_aeropara_in_filename
245  character(len=H_LONG) :: atmos_phy_rd_mstrn_hygropara_in_filename
246  integer :: atmos_phy_rd_mstrn_nband
247  integer :: atmos_phy_rd_mstrn_nptype
248  integer :: atmos_phy_rd_mstrn_nradius
249  integer :: atmos_phy_rd_mstrn_nradius_cloud = -1
250  integer :: atmos_phy_rd_mstrn_nradius_aero = -1
251 
252  namelist / param_atmos_phy_rd_mstrn / &
253  atmos_phy_rd_mstrn_toa, &
254  atmos_phy_rd_mstrn_kadd, &
255  atmos_phy_rd_mstrn_gaspara_in_filename, &
256  atmos_phy_rd_mstrn_aeropara_in_filename, &
257  atmos_phy_rd_mstrn_hygropara_in_filename, &
258  atmos_phy_rd_mstrn_nband, &
259  atmos_phy_rd_mstrn_nptype, &
260  atmos_phy_rd_mstrn_nradius, &
261  atmos_phy_rd_mstrn_nradius_cloud, &
262  atmos_phy_rd_mstrn_nradius_aero, &
263  atmos_phy_rd_mstrn_only_qci, &
264  atmos_phy_rd_mstrn_only_tropocloud, &
265  atmos_phy_rd_mstrn_use_aero
266 
267  integer :: kmax
268  integer :: ngas, ncfc
269  integer :: ierr
270  !---------------------------------------------------------------------------
271 
272  log_newline
273  log_info("ATMOS_PHY_RD_mstrnx_setup",*) 'Setup'
274  log_info("ATMOS_PHY_RD_mstrnx_setup",*) 'Sekiguchi and Nakajima (2008) mstrnX radiation process'
275 
276  !--- read namelist
277  atmos_phy_rd_mstrn_toa = rd_toa
278  atmos_phy_rd_mstrn_kadd = rd_kadd
279  atmos_phy_rd_mstrn_gaspara_in_filename = mstrn_gaspara_inputfile
280  atmos_phy_rd_mstrn_aeropara_in_filename = mstrn_aeropara_inputfile
281  atmos_phy_rd_mstrn_hygropara_in_filename = mstrn_hygropara_inputfile
282  atmos_phy_rd_mstrn_nband = mstrn_nband
283  atmos_phy_rd_mstrn_nptype = mstrn_nptype
284  atmos_phy_rd_mstrn_nradius = mstrn_nradius
285 
286  rewind(io_fid_conf)
287  read(io_fid_conf,nml=param_atmos_phy_rd_mstrn,iostat=ierr)
288  if( ierr < 0 ) then !--- missing
289  log_info("ATMOS_PHY_RD_mstrnx_setup",*) 'Not found namelist. Default used.'
290  elseif( ierr > 0 ) then !--- fatal error
291  log_error("ATMOS_PHY_RD_mstrnx_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_RD_MSTRN. Check!'
292  call prc_abort
293  endif
294  log_nml(param_atmos_phy_rd_mstrn)
295 
296  rd_toa = atmos_phy_rd_mstrn_toa
297  rd_kadd = atmos_phy_rd_mstrn_kadd
298  mstrn_gaspara_inputfile = atmos_phy_rd_mstrn_gaspara_in_filename
299  mstrn_aeropara_inputfile = atmos_phy_rd_mstrn_aeropara_in_filename
300  mstrn_hygropara_inputfile = atmos_phy_rd_mstrn_hygropara_in_filename
301  mstrn_nband = atmos_phy_rd_mstrn_nband
302  mstrn_nptype = atmos_phy_rd_mstrn_nptype
303 
304  if ( atmos_phy_rd_mstrn_nradius_cloud < 0 ) then ! set default
305  atmos_phy_rd_mstrn_nradius_cloud = atmos_phy_rd_mstrn_nradius
306  endif
307  if ( atmos_phy_rd_mstrn_nradius_aero < 0 ) then ! set default
308  atmos_phy_rd_mstrn_nradius_aero = atmos_phy_rd_mstrn_nradius
309  endif
310 
311  mstrn_nradius = max( atmos_phy_rd_mstrn_nradius, &
312  atmos_phy_rd_mstrn_nradius_cloud, &
313  atmos_phy_rd_mstrn_nradius_aero )
314 
315  !--- setup MP/AE - RD(particle) coupling parameter
316  allocate( ptype_nradius(mstrn_nptype) )
317 
318  if ( mstrn_nptype == 9 ) then
319 
320  i_mpae2rd( 1) = 1 ! cloud water Qc -> Water
321  i_mpae2rd( 2) = 1 ! rain Qr -> Water
322  i_mpae2rd( 3) = 2 ! cloud ice Qi -> Ice
323  i_mpae2rd( 4) = 2 ! snow Qs -> Ice
324  i_mpae2rd( 5) = 2 ! graupel Qg -> Ice
325  i_mpae2rd( 6) = 2 ! hail Qh -> Ice
326  i_mpae2rd( 7) = 3 ! aerosol type1 -> Soil dust
327  i_mpae2rd( 8) = 4 ! aerosol type2 -> Carbonacerous (BC/OC=0.3)
328  i_mpae2rd( 9) = 5 ! aerosol type3 -> Carbonacerous (BC/OC=0.15)
329  i_mpae2rd(10) = 6 ! aerosol type4 -> Carbonacerous (BC/OC=0.)
330  i_mpae2rd(11) = 7 ! aerosol type5 -> Black carbon
331  i_mpae2rd(12) = 8 ! aerosol type6 -> Sulfate
332  i_mpae2rd(13) = 9 ! aerosol type7 -> Sea salt
333 
334  ptype_nradius( 1) = atmos_phy_rd_mstrn_nradius_cloud
335  ptype_nradius( 2) = atmos_phy_rd_mstrn_nradius_cloud
336  ptype_nradius( 3) = atmos_phy_rd_mstrn_nradius_aero
337  ptype_nradius( 4) = atmos_phy_rd_mstrn_nradius_aero
338  ptype_nradius( 5) = atmos_phy_rd_mstrn_nradius_aero
339  ptype_nradius( 6) = atmos_phy_rd_mstrn_nradius_aero
340  ptype_nradius( 7) = atmos_phy_rd_mstrn_nradius_aero
341  ptype_nradius( 8) = atmos_phy_rd_mstrn_nradius_aero
342  ptype_nradius( 9) = atmos_phy_rd_mstrn_nradius_aero
343 
344  elseif( mstrn_nptype == 12 ) then
345 
346  i_mpae2rd( 1) = 1 ! cloud water Qc -> CLOUD
347  i_mpae2rd( 2) = 2 ! rain Qr -> RAIN
348  i_mpae2rd( 3) = 3 ! cloud ice Qi -> ICE
349  i_mpae2rd( 4) = 4 ! snow Qs -> SNOW
350  i_mpae2rd( 5) = 5 ! graupel Qg -> GRAUPEL
351  i_mpae2rd( 6) = 5 ! hail Qh -> GRAUPEL
352  i_mpae2rd( 7) = 6 ! aerosol type1 -> Soil dust
353  i_mpae2rd( 8) = 7 ! aerosol type2 -> Carbonacerous (BC/OC=0.3)
354  i_mpae2rd( 9) = 8 ! aerosol type3 -> Carbonacerous (BC/OC=0.15)
355  i_mpae2rd(10) = 9 ! aerosol type4 -> Carbonacerous (BC/OC=0.)
356  i_mpae2rd(11) = 10 ! aerosol type5 -> Black carbon
357  i_mpae2rd(12) = 11 ! aerosol type6 -> Sulfate
358  i_mpae2rd(13) = 12 ! aerosol type7 -> Sea salt
359 
360  ptype_nradius( 1) = atmos_phy_rd_mstrn_nradius_cloud
361  ptype_nradius( 2) = atmos_phy_rd_mstrn_nradius_cloud
362  ptype_nradius( 3) = atmos_phy_rd_mstrn_nradius_cloud
363  ptype_nradius( 4) = atmos_phy_rd_mstrn_nradius_cloud
364  ptype_nradius( 5) = atmos_phy_rd_mstrn_nradius_cloud
365  ptype_nradius( 6) = atmos_phy_rd_mstrn_nradius_aero
366  ptype_nradius( 7) = atmos_phy_rd_mstrn_nradius_aero
367  ptype_nradius( 8) = atmos_phy_rd_mstrn_nradius_aero
368  ptype_nradius( 9) = atmos_phy_rd_mstrn_nradius_aero
369  ptype_nradius(10) = atmos_phy_rd_mstrn_nradius_aero
370  ptype_nradius(11) = atmos_phy_rd_mstrn_nradius_aero
371  ptype_nradius(12) = atmos_phy_rd_mstrn_nradius_aero
372 
373  endif
374 
375  !--- setup MSTRN parameter
376  call rd_mstrn_setup( ngas, & ! [OUT]
377  ncfc ) ! [OUT]
378 
379  !--- setup climatological profile
380  call rd_profile_setup
381 
382  kmax = ke - ks + 1
383  rd_kmax = kmax + rd_kadd
384 
385  !--- allocate arrays
386  ! input
387  allocate( rd_zh(rd_kmax+1) )
388  allocate( rd_z(rd_kmax ) )
389 
390  allocate( rd_rhodz(rd_kmax ) )
391  allocate( rd_pres(rd_kmax ) )
392  allocate( rd_presh(rd_kmax+1) )
393  allocate( rd_temp(rd_kmax ) )
394  allocate( rd_temph(rd_kmax+1) )
395 
396  allocate( rd_gas(rd_kmax,ngas ) )
397  allocate( rd_cfc(rd_kmax,ncfc ) )
398  allocate( rd_aerosol_conc(rd_kmax,rd_naero) )
399  allocate( rd_aerosol_radi(rd_kmax,rd_naero) )
400  allocate( rd_cldfrac(rd_kmax ) )
401 
402  !--- setup vartical grid for radiation (larger TOA than Model domain)
403  call rd_profile_setup_zgrid( &
404  ka, ks, ke, &
405  rd_kmax, rd_kadd, & ! [IN]
406  rd_toa, cz(:), fz(:), & ! [IN]
407  rd_zh(:), rd_z(:) ) ! [OUT]
408 
409  !--- read climatological profile
410  call rd_profile_read( rd_kmax, & ! [IN]
411  ngas, & ! [IN]
412  ncfc, & ! [IN]
413  rd_naero, & ! [IN]
414  atmos_grid_cartesc_real_basepoint_lat, & ! [IN]
415  time_nowdate(:), & ! [IN]
416  rd_zh(:), & ! [IN]
417  rd_z(:), & ! [IN]
418  rd_rhodz(:), & ! [OUT]
419  rd_pres(:), & ! [OUT]
420  rd_presh(:), & ! [OUT]
421  rd_temp(:), & ! [OUT]
422  rd_temph(:), & ! [OUT]
423  rd_gas(:,:), & ! [OUT]
424  rd_cfc(:,:), & ! [OUT]
425  rd_aerosol_conc(:,:), & ! [OUT]
426  rd_aerosol_radi(:,:), & ! [OUT]
427  rd_cldfrac(:) ) ! [OUT]
428 
429  return
430  end subroutine atmos_phy_rd_mstrnx_setup
431 
432  !-----------------------------------------------------------------------------
434  subroutine atmos_phy_rd_mstrnx_flux( &
435  KA, KS, KE, IA, IS, IE, JA, JS, JE, &
436  DENS, TEMP, PRES, QV, &
437  CZ, FZ, &
438  fact_ocean, &
439  fact_land, &
440  fact_urban, &
441  temp_sfc, albedo_sfc, &
442  solins, cosSZA, &
443  CLDFRAC, MP_Re, MP_Qe, &
444  AE_Re, AE_Qe, &
445  flux_rad, &
446  flux_rad_top, &
447  flux_rad_sfc_dn, &
448  dtau_s, dem_s )
449  !Jval )
450  use scale_const, only: &
451  eps => const_eps, &
452  mdry => const_mdry, &
453  mvap => const_mvap, &
454  ppm => const_ppm
455  use scale_time, only: &
456  time_nowdate
457  use scale_atmos_grid_cartesc_real, only: &
458  atmos_grid_cartesc_real_basepoint_lat
459  use scale_atmos_hydrometeor, only: &
460  n_hyd, &
461  i_hc, &
462  i_hi, &
463  hyd_dens
464  use scale_atmos_aerosol, only: &
465  n_ae, &
466  ae_dens
467  use scale_atmos_phy_rd_profile, only: &
468  rd_profile_read => atmos_phy_rd_profile_read, &
469  rd_profile_use_climatology => atmos_phy_rd_profile_use_climatology
470  implicit none
471  integer, intent(in) :: ka, ks, ke
472  integer, intent(in) :: ia, is, ie
473  integer, intent(in) :: ja, js, je
474 
475  real(rp), intent(in) :: dens (ka,ia,ja)
476  real(rp), intent(in) :: temp (ka,ia,ja)
477  real(rp), intent(in) :: pres (ka,ia,ja)
478  real(rp), intent(in) :: qv (ka,ia,ja)
479  real(rp), intent(in) :: cz ( ka,ia,ja) ! UNUSED
480  real(rp), intent(in) :: fz (0:ka,ia,ja)
481  real(rp), intent(in) :: fact_ocean (ia,ja)
482  real(rp), intent(in) :: fact_land (ia,ja)
483  real(rp), intent(in) :: fact_urban (ia,ja)
484  real(rp), intent(in) :: temp_sfc (ia,ja)
485  real(rp), intent(in) :: albedo_sfc (ia,ja,n_rad_dir,n_rad_rgn)
486  real(rp), intent(in) :: solins (ia,ja)
487  real(rp), intent(in) :: cossza (ia,ja)
488  real(rp), intent(in) :: cldfrac (ka,ia,ja)
489  real(rp), intent(in) :: mp_re (ka,ia,ja,n_hyd)
490  real(rp), intent(in) :: mp_qe (ka,ia,ja,n_hyd)
491  real(rp), intent(in) :: ae_re (ka,ia,ja,n_ae)
492  real(rp), intent(in) :: ae_qe (ka,ia,ja,n_ae)
493  real(rp), intent(out) :: flux_rad (ka,ia,ja,2,2,2)
494  real(rp), intent(out) :: flux_rad_top (ia,ja,2,2,2)
495  real(rp), intent(out) :: flux_rad_sfc_dn(ia,ja,n_rad_dir,n_rad_rgn)
496 
497  real(rp), intent(out), optional :: dtau_s(ka,ia,ja) ! 0.67 micron cloud optical depth
498  real(rp), intent(out), optional :: dem_s (ka,ia,ja) ! 10.5 micron cloud emissivity
499 ! real(RP), intent(out), optional :: Jval (KA,IA,JA,CH_QA_photo)
500 
501  integer :: tropopause(ia,ja)
502  real(rp) :: gamma
503 
504  real(rp), parameter :: min_cldfrac = 1.e-8_rp
505 
506  real(rp) :: rhodz_merge (rd_kmax,ia,ja)
507  real(rp) :: pres_merge (rd_kmax,ia,ja)
508  real(rp) :: temp_merge (rd_kmax,ia,ja)
509  real(rp) :: temph_merge (rd_kmax+1,ia,ja)
510 
511  real(rp) :: gas_merge (rd_kmax,ia,ja,mstrn_ngas)
512  real(rp) :: cfc_merge (rd_kmax,ia,ja,mstrn_ncfc)
513  real(rp) :: aerosol_conc_merge(rd_kmax,ia,ja,rd_naero )
514  real(rp) :: aerosol_radi_merge(rd_kmax,ia,ja,rd_naero )
515  real(rp) :: cldfrac_merge (rd_kmax,ia,ja)
516 
517  ! output
518  real(rp) :: flux_rad_merge(rd_kmax+1,ia,ja,2,2,mstrn_ncloud)
519  real(rp) :: taucld_067u (rd_kmax,ia,ja) ! 0.67 micron cloud optical depth
520  real(rp) :: emiscld_105u (rd_kmax,ia,ja) ! 10.5 micron cloud emissivity
521 
522  real(rp) :: zerosw
523 
524  integer :: ihydro, iaero
525  integer :: rd_k, k, i, j, v, ic
526  !---------------------------------------------------------------------------
527 
528  log_progress(*) 'atmosphere / physics / radiation / mstrnX'
529 
530  call prof_rapstart('RD_Profile', 3)
531 
532  if ( atmos_phy_rd_mstrn_only_tropocloud ) then
533  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
534  !$omp private(k,i,j, &
535  !$omp gamma) &
536  !$omp shared (tropopause,pres,temp,CZ, &
537  !$omp KS,KE,IS,IE,JS,JE)
538  do j = js, je
539  do i = is, ie
540  tropopause(i,j) = ke+1
541  do k = ke, ks, -1
542  if ( pres(k,i,j) >= 300.e+2_rp ) then
543  exit
544  elseif( pres(k,i,j) < 50.e+2_rp ) then
545  tropopause(i,j) = k
546  else
547  gamma = ( temp(k+1,i,j) - temp(k,i,j) ) / ( cz(k+1,i,j) - cz(k,i,j) )
548  if( gamma > 1.e-4_rp ) tropopause(i,j) = k
549  endif
550  enddo
551  enddo
552  enddo
553  else
554  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
555  !$omp private(i,j) &
556  !$omp shared (tropopause, &
557  !$omp KE,IS,IE,JS,JE)
558 !OCL XFILL
559  do j = js, je
560  do i = is, ie
561  tropopause(i,j) = ke+1
562  end do
563  end do
564  endif
565 
566  ! marge basic profile and value in model domain
567 
568  if ( rd_profile_use_climatology ) then
569  call rd_profile_read( rd_kmax, & ! [IN]
570  mstrn_ngas, & ! [IN]
571  mstrn_ncfc, & ! [IN]
572  rd_naero, & ! [IN]
573  atmos_grid_cartesc_real_basepoint_lat, & ! [IN]
574  time_nowdate(:), & ! [IN]
575  rd_zh(:), & ! [IN]
576  rd_z(:), & ! [IN]
577  rd_rhodz(:), & ! [OUT]
578  rd_pres(:), & ! [OUT]
579  rd_presh(:), & ! [OUT]
580  rd_temp(:), & ! [OUT]
581  rd_temph(:), & ! [OUT]
582  rd_gas(:,:), & ! [OUT]
583  rd_cfc(:,:), & ! [OUT]
584  rd_aerosol_conc(:,:), & ! [OUT]
585  rd_aerosol_radi(:,:), & ! [OUT]
586  rd_cldfrac(:) ) ! [OUT]
587  endif
588 
589 !OCL XFILL
590  !$omp parallel do default(none) &
591  !$omp shared(JS,JE,IS,IE,RD_KADD,temph_merge,RD_temph,KE,RD_KMAX,KS,temp,CZ,FZ) &
592  !$omp private(i,j,k,RD_k) OMP_SCHEDULE_ collapse(2)
593  do j = js, je
594  do i = is, ie
595  do rd_k = 1, rd_kadd
596  temph_merge(rd_k,i,j) = rd_temph(rd_k)
597  enddo
598 
599  temph_merge(rd_kadd+1,i,j) = temp(ke,i,j)
600  do rd_k = rd_kadd+2, rd_kmax
601  k = ks + rd_kmax - rd_k ! reverse axis
602 
603  temph_merge(rd_k,i,j) = ( temp(k ,i,j) * (cz(k+1,i,j)-fz(k,i,j)) / (cz(k+1,i,j)-cz(k,i,j)) &
604  + temp(k+1,i,j) * (fz(k ,i,j)-cz(k,i,j)) / (cz(k+1,i,j)-cz(k,i,j)) )
605  enddo
606  temph_merge(rd_kmax+1,i,j) = temp(ks,i,j)
607  enddo
608  enddo
609 
610 !OCL XFILL
611  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
612  !$omp private(k,i,j,RD_k) &
613  !$omp shared (RD_temp,dens,FZ,pres,temp, &
614  !$omp rhodz_merge,RD_rhodz,pres_merge,RD_pres,temp_merge, &
615  !$omp KS,JS,JE,IS,IE,RD_KMAX,RD_KADD)
616  do j = js, je
617  do i = is, ie
618  do rd_k = 1, rd_kadd
619  rhodz_merge(rd_k,i,j) = rd_rhodz(rd_k)
620  pres_merge(rd_k,i,j) = rd_pres(rd_k)
621  temp_merge(rd_k,i,j) = rd_temp(rd_k)
622  enddo
623 
624  do rd_k = rd_kadd+1, rd_kmax
625  k = ks + rd_kmax - rd_k ! reverse axis
626 
627  rhodz_merge(rd_k,i,j) = dens(k,i,j) * ( fz(k,i,j)-fz(k-1,i,j) ) ! [kg/m2]
628  pres_merge(rd_k,i,j) = pres(k,i,j) * 1.e-2_rp ! [hPa]
629  temp_merge(rd_k,i,j) = temp(k,i,j)
630  enddo
631  enddo
632  enddo
633 
634 !OCL XFILL
635 !OCL SERIAL
636  do v = 1, mstrn_ngas
637 !OCL PARALLEL
638  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
639  !$omp private(i,j,RD_k) &
640  !$omp shared (gas_merge,RD_gas, &
641  !$omp IS,IE,JS,JE,RD_KMAX,v)
642  do j = js, je
643  do i = is, ie
644  do rd_k = 1, rd_kmax
645  gas_merge(rd_k,i,j,v) = rd_gas(rd_k,v)
646  enddo
647  enddo
648  enddo
649  enddo
650 
651  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
652  !$omp private(k,i,j,RD_k, &
653  !$omp zerosw ) &
654  !$omp shared (gas_merge,QV,EPS,Mvap,Mdry, &
655  !$omp KS,IS,IE,JS,JE,RD_KADD,RD_KMAX)
656  do j = js, je
657  do i = is, ie
658  do rd_k = rd_kadd+1, rd_kmax
659  k = ks + rd_kmax - rd_k ! reverse axis
660  zerosw = sign(0.5_rp, qv(k,i,j)-eps) + 0.5_rp
661  gas_merge(rd_k,i,j,1) = qv(k,i,j) / mvap * mdry / ppm * zerosw ! [PPM]
662  enddo
663  enddo
664  enddo
665 
666 !OCL XFILL
667 !OCL SERIAL
668  do v = 1, mstrn_ncfc
669  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
670  !$omp private(i,j,RD_k) &
671  !$omp shared (cfc_merge,RD_cfc, &
672  !$omp IS,IE,JS,JE,RD_KMAX,v)
673 !OCL PARALLEL
674  do j = js, je
675  do i = is, ie
676  do rd_k = 1, rd_kmax
677  cfc_merge(rd_k,i,j,v) = rd_cfc(rd_k,v)
678  enddo
679  enddo
680  enddo
681  enddo
682 
683  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
684  !$omp private(k,i,j,RD_k) &
685  !$omp shared (cldfrac_merge,RD_cldfrac,cldfrac, &
686  !$omp KS,IS,IE,JS,JE,RD_KMAX,RD_KADD)
687  do j = js, je
688  do i = is, ie
689  do rd_k = 1, rd_kadd
690  cldfrac_merge(rd_k,i,j) = rd_cldfrac(rd_k)
691  enddo
692 
693  do rd_k = rd_kadd+1, rd_kmax
694  k = ks + rd_kmax - rd_k ! reverse axis
695 
696  cldfrac_merge(rd_k,i,j) = min( max( cldfrac(k,i,j), 0.0_rp ), 1.0_rp )
697 ! cldfrac_merge(RD_k,i,j) = 0.5_RP + sign( 0.5_RP, cldfrac(k,i,j)-min_cldfrac )
698  enddo
699  enddo
700  enddo
701 
702 !OCL XFILL
703 !OCL SERIAL
704  do v = 1, rd_naero
705  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
706  !$omp private(i,j,RD_k) &
707  !$omp shared (aerosol_conc_merge,RD_aerosol_conc,aerosol_radi_merge,RD_aerosol_radi, &
708  !$omp IS,IE,JS,JE,RD_KADD,v)
709 !OCL PARALLEL
710  do j = js, je
711  do i = is, ie
712  do rd_k = 1, rd_kadd
713  aerosol_conc_merge(rd_k,i,j,v) = rd_aerosol_conc(rd_k,v)
714  aerosol_radi_merge(rd_k,i,j,v) = rd_aerosol_radi(rd_k,v)
715  enddo
716  enddo
717  enddo
718  enddo
719 
720 !OCL SERIAL
721  do ihydro = 1, n_hyd
722  if ( atmos_phy_rd_mstrn_only_qci .and. &
723  ( ihydro /= i_hc .and. ihydro /= i_hi ) ) then
724 !OCL XFILL
725  aerosol_conc_merge(:,:,:,ihydro) = 0.0_rp
726  cycle
727  end if
728  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
729  !$omp private(k,i,j,RD_k) &
730  !$omp shared (aerosol_conc_merge,tropopause,MP_Qe,HYD_DENS,RHO_std, &
731  !$omp KS,KE,IS,IE,JS,JE,RD_KMAX,RD_KADD,ihydro)
732 !OCL PARALLEL
733  do j = js, je
734  do i = is, ie
735  do rd_k = rd_kadd+1, rd_kadd+1 + ke - tropopause(i,j)
736  aerosol_conc_merge(rd_k,i,j,ihydro) = 0.0_rp
737  end do
738  do rd_k = rd_kadd+1 + ke - tropopause(i,j) + 1, rd_kmax
739  k = ks + rd_kmax - rd_k ! reverse axis
740  aerosol_conc_merge(rd_k,i,j,ihydro) = max( mp_qe(k,i,j,ihydro), 0.0_rp ) &
741  / hyd_dens(ihydro) * rho_std / ppm ! [PPM to standard air]
742  enddo
743  enddo
744  enddo
745  enddo
746 
747 !OCL SERIAL
748  do ihydro = 1, n_hyd
749  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
750  !$omp private(k,i,j,RD_k) &
751  !$omp shared (aerosol_radi_merge,MP_Re, &
752  !$omp KS,IS,IE,JS,JE,RD_KMAX,RD_KADD,ihydro)
753 !OCL PARALLEL
754  do j = js, je
755  do i = is, ie
756  do rd_k = rd_kadd+1, rd_kmax
757  k = ks + rd_kmax - rd_k ! reverse axis
758  aerosol_radi_merge(rd_k,i,j,ihydro) = mp_re(k,i,j,ihydro)
759  enddo
760  enddo
761  enddo
762  enddo
763 
764 !OCL SERIAL
765  do iaero = 1, n_ae
766 
767  if ( atmos_phy_rd_mstrn_use_aero ) then
768  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
769  !$omp private(k,i,j,RD_k) &
770  !$omp shared (aerosol_conc_merge,aerosol_radi_merge,AE_Qe,RHO_std,AE_Re, &
771  !$omp KS,IS,IE,JS,JE,RD_KMAX,RD_KADD,iaero)
772 !OCL PARALLEL
773  do j = js, je
774  do i = is, ie
775  do rd_k = rd_kadd+1, rd_kmax
776  k = ks + rd_kmax - rd_k ! reverse axis
777  aerosol_conc_merge(rd_k,i,j,n_hyd+iaero) = max( ae_qe(k,i,j,iaero), 0.0_rp ) &
778  / ae_dens(iaero) * rho_std / ppm ! [PPM to standard air]
779  aerosol_radi_merge(rd_k,i,j,n_hyd+iaero) = ae_re(k,i,j,iaero)
780  enddo
781  enddo
782  enddo
783  else
784  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
785  !$omp private(i,j,RD_k) &
786  !$omp shared (aerosol_conc_merge,RD_aerosol_conc,aerosol_radi_merge,RD_aerosol_radi, &
787  !$omp IS,IE,JS,JE,RD_KMAX,RD_KADD,iaero)
788 !OCL PARALLEL
789  do j = js, je
790  do i = is, ie
791  do rd_k = rd_kadd+1, rd_kmax
792  aerosol_conc_merge(rd_k,i,j,n_hyd+iaero) = rd_aerosol_conc(rd_k,n_hyd+iaero)
793  aerosol_radi_merge(rd_k,i,j,n_hyd+iaero) = rd_aerosol_radi(rd_k,n_hyd+iaero)
794  enddo
795  enddo
796  enddo
797  endif
798 
799  enddo
800 
801  call prof_rapend ('RD_Profile', 3)
802  call prof_rapstart('RD_MSTRN_DTRN3', 3)
803 
804  ! calc radiative transfer
805  call rd_mstrn_dtrn3( rd_kmax, ia, is, ie, ja, js, je, & ! [IN]
806  mstrn_ngas, mstrn_ncfc, rd_naero, & ! [IN]
807  rd_hydro_str, rd_hydro_end, rd_aero_str, rd_aero_end, & ! [IN]
808  solins(:,:), cossza(:,:), & ! [IN]
809  rhodz_merge(:,:,:), pres_merge(:,:,:), & ! [IN]
810  temp_merge(:,:,:), temph_merge(:,:,:), temp_sfc(:,:), & ! [IN]
811  gas_merge(:,:,:,:), cfc_merge(:,:,:,:), & ! [IN]
812  aerosol_conc_merge(:,:,:,:), aerosol_radi_merge(:,:,:,:), & ! [IN]
813  i_mpae2rd(:), & ! [IN]
814  cldfrac_merge(:,:,:), & ! [IN]
815  albedo_sfc(:,:,:,:), & ! [IN]
816  fact_ocean(:,:), fact_land(:,:), fact_urban(:,:), & ! [IN]
817  flux_rad_merge(:,:,:,:,:,:), flux_rad_sfc_dn(:,:,:,:), & ! [OUT]
818  taucld_067u(:,:,:), emiscld_105u(:,:,:) ) ! [OUT]
819 
820  call prof_rapend ('RD_MSTRN_DTRN3', 3)
821 
822  ! return to grid coordinate of model domain
823 !OCL SERIAL
824  do ic = 1, 2
825  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
826  !$omp private(k,i,j,RD_k) &
827  !$omp shared (flux_rad,flux_rad_merge, &
828  !$omp KS,IS,IE,JS,JE,RD_KMAX,RD_KADD,ic)
829 !OCL PARALLEL
830  do j = js, je
831  do i = is, ie
832  do rd_k = rd_kadd+1, rd_kmax+1
833  k = ks + rd_kmax - rd_k ! reverse axis
834 
835  flux_rad(k,i,j,i_lw,i_up,ic) = flux_rad_merge(rd_k,i,j,i_lw,i_up,ic)
836  flux_rad(k,i,j,i_lw,i_dn,ic) = flux_rad_merge(rd_k,i,j,i_lw,i_dn,ic)
837  flux_rad(k,i,j,i_sw,i_up,ic) = flux_rad_merge(rd_k,i,j,i_sw,i_up,ic)
838  flux_rad(k,i,j,i_sw,i_dn,ic) = flux_rad_merge(rd_k,i,j,i_sw,i_dn,ic)
839  enddo
840  enddo
841  enddo
842  enddo
843 
844 !OCL XFILL
845 !OCL SERIAL
846  do ic = 1, 2
847  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
848  !$omp private(i,j) &
849  !$omp shared (flux_rad_top,flux_rad_merge, &
850  !$omp IS,IE,JS,JE,ic)
851 !OCL PARALLEL
852  do j = js, je
853  do i = is, ie
854  flux_rad_top(i,j,i_lw,i_up,ic) = flux_rad_merge(1,i,j,i_lw,i_up,ic)
855  flux_rad_top(i,j,i_lw,i_dn,ic) = flux_rad_merge(1,i,j,i_lw,i_dn,ic)
856  flux_rad_top(i,j,i_sw,i_up,ic) = flux_rad_merge(1,i,j,i_sw,i_up,ic)
857  flux_rad_top(i,j,i_sw,i_dn,ic) = flux_rad_merge(1,i,j,i_sw,i_dn,ic)
858  enddo
859  enddo
860  enddo
861 
862  if ( present( dtau_s ) ) then
863  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
864  !$omp private(k,i,j,RD_k) &
865  !$omp shared (dtau_s,tauCLD_067u, &
866  !$omp KS,IS,IE,JS,JE,RD_KMAX,RD_KADD)
867  do j = js, je
868  do i = is, ie
869  do rd_k = rd_kadd+1, rd_kmax
870  k = ks + rd_kmax - rd_k ! reverse axis
871  dtau_s(k,i,j) = taucld_067u(rd_k,i,j)
872  enddo
873  enddo
874  enddo
875  end if
876 
877  if ( present( dem_s ) ) then
878  !$omp parallel do default(none) OMP_SCHEDULE_ collapse(2) &
879  !$omp private(k,i,j,RD_k) &
880  !$omp shared (dem_s,emisCLD_105u, &
881  !$omp KS,IS,IE,JS,JE,RD_KMAX,RD_KADD)
882  do j = js, je
883  do i = is, ie
884  do rd_k = rd_kadd+1, rd_kmax
885  k = ks + rd_kmax - rd_k ! reverse axis
886  dem_s(k,i,j) = emiscld_105u(rd_k,i,j)
887  enddo
888  enddo
889  enddo
890  end if
891 
892  return
893  end subroutine atmos_phy_rd_mstrnx_flux
894 
895  !-----------------------------------------------------------------------------
897  subroutine rd_mstrn_setup( &
898  ngas, &
899  ncfc )
900  use scale_prc, only: &
901  prc_abort
902  use scale_const, only: &
903  grav => const_grav, &
904  rdry => const_rdry, &
905  pstd => const_pstd, &
906  tem00 => const_tem00
907  implicit none
908 
909  integer, intent(out) :: ngas
910  integer, intent(out) :: ncfc
911 
912  real(rp) :: acfc(mstrn_ncfc)
913 
914  integer :: nband, nstream, nfitp, nfitt, nflag
915  integer :: nsfc, nptype, nplkord, nfitplk
916  integer :: nradius
917 
918  character(len=H_LONG) :: dummy
919 
920  integer :: fid, ierr
921  integer :: iw, ich, ip, it, igas, icfc, iptype, im
922  !---------------------------------------------------------------------------
923 
924  !---< gas absorption parameter input >---
925  ngas = mstrn_ngas
926 
927  ! allocate arrays
928  allocate( waveh(mstrn_nband+1) )
929  allocate( logfitp(mstrn_nfitp) )
930  allocate( fitt(mstrn_nfitt) )
931  allocate( logfitt(mstrn_nfitt) )
932 
933  allocate( iflgb(mstrn_nflag, mstrn_nband) )
934  allocate( nch( mstrn_nband) )
935  allocate( wgtch(mstrn_ch_limit,mstrn_nband) )
936  allocate( ngasabs( mstrn_nband) )
937  allocate( igasabs(mstrn_ngas, mstrn_nband) )
938 
939  allocate( akd(mstrn_ch_limit,mstrn_nfitp,mstrn_nfitt,mstrn_ngas,mstrn_nband) )
940  allocate( skd(mstrn_ch_limit,mstrn_nfitp,mstrn_nfitt, mstrn_nband) )
941  allocate( acfc_pow(mstrn_ncfc,mstrn_nband) )
942 
943  fid = io_get_available_fid()
944 
945  open( fid, &
946  file = trim(mstrn_gaspara_inputfile), &
947  form = 'formatted', &
948  status = 'old', &
949  iostat = ierr )
950 
951  if ( ierr /= 0 ) then
952  log_error("RD_MSTRN_setup",*) 'Input data file does not found! ', trim(mstrn_gaspara_inputfile)
953  call prc_abort
954  endif
955 
956  ! read gas parameters for check
957  read(fid,*) nband, nstream, nfitp, nfitt, nflag, ncfc
958 
959  if ( nband /= mstrn_nband ) then
960  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nband(given,file)=', mstrn_nband, nband
961  call prc_abort
962  endif
963  if ( nstream /= mstrn_nstream ) then
964  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nstream(given,file)=', mstrn_nstream, nstream
965  call prc_abort
966  endif
967  if ( nfitp /= mstrn_nfitp ) then
968  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nfitP(given,file)=', mstrn_nfitp, nfitp
969  call prc_abort
970  endif
971  if ( nfitt /= mstrn_nfitt ) then
972  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nfitT(given,file)=', mstrn_nfitt, nfitt
973  call prc_abort
974  endif
975  if ( nflag /= mstrn_nflag ) then
976  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nflag(given,file)=', mstrn_nflag, nflag
977  call prc_abort
978  endif
979  if ( ncfc /= mstrn_ncfc ) then
980  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! ncfc(given,file)=', mstrn_ncfc, ncfc
981  call prc_abort
982  endif
983 
984  ! wave band boundaries
985  read(fid,*) dummy
986  read(fid,*) waveh(:)
987  ! fitting point for log(pressure)
988  read(fid,*) dummy
989  read(fid,*) logfitp(:)
990  ! fitting point for temperature
991  read(fid,*) dummy
992  read(fid,*) fitt(:)
993 
994  logfitt(:) = log10( fitt(:) )
995 
996  ! for each band
997  do iw = 1, mstrn_nband
998 
999  ! optical properties flag
1000  read(fid,*) dummy
1001  read(fid,*) iflgb(:,iw)
1002  ! number of subintervals
1003  read(fid,*) dummy
1004  read(fid,*) nch(iw)
1005  ! weights for subintervals
1006  read(fid,*) dummy
1007  read(fid,*) (wgtch(ich,iw),ich=1,nch(iw))
1008  ! number of considering gases
1009  read(fid,*) dummy
1010  read(fid,*) ngasabs(iw)
1011 
1012  ! major gas absorption
1013  if ( ngasabs(iw) > 0 ) then
1014  do igas = 1, ngasabs(iw)
1015  read(fid,*) igasabs(igas,iw)
1016  do it = 1, mstrn_nfitt
1017  do ip = 1, mstrn_nfitp
1018  read(fid,*) (akd(ich,ip,it,igasabs(igas,iw),iw),ich=1,nch(iw))
1019  enddo
1020  enddo
1021  enddo
1022  endif
1023 
1024  ! H2O continuum
1025  if ( iflgb(i_h2o_continuum,iw) > 0 ) then
1026  read(fid,*) dummy
1027  do it = 1, mstrn_nfitt
1028  do ip = 1, mstrn_nfitp
1029  read(fid,*) (skd(ich,ip,it,iw),ich=1,nch(iw))
1030  enddo
1031  enddo
1032  endif
1033 
1034  ! CFC absorption
1035  if ( iflgb(i_cfc_continuum,iw) > 0 ) then
1036  read(fid,*) dummy
1037  read(fid,*) (acfc(icfc),icfc=1,mstrn_ncfc)
1038 
1039  do icfc = 1, mstrn_ncfc
1040  acfc_pow(icfc,iw) = 10.0_rp**acfc(icfc)
1041  end do
1042  endif
1043 
1044  enddo ! band loop
1045 
1046  close(fid)
1047 
1048  !---< aerosol(particle) parameter input >---
1049 
1050  ! allocate arrays
1051  allocate( fitplk(mstrn_nfitplk,mstrn_nband) )
1052  allocate( fsol( mstrn_nband) )
1053  allocate( sfc(mstrn_nsfc, mstrn_nband) )
1054  allocate( rayleigh( mstrn_nband) )
1055 
1056  allocate( qmol( mstrn_nmoment,mstrn_nband) )
1057  allocate( q(-1:mstrn_nradius+1,mstrn_nptype,mstrn_nmoment,mstrn_nband) )
1058  do iptype = 1, mstrn_nptype
1059  nradius = ptype_nradius(iptype)
1060 
1061  q( -1:0 ,iptype,:,:) = 0.d0 ! dummy for NaN
1062  q(nradius+1:mstrn_nradius+1,iptype,:,:) = 0.d0 ! dummy for extrapolation
1063  enddo
1064 
1065 
1066  open( fid, &
1067  file = trim(mstrn_aeropara_inputfile), &
1068  form = 'formatted', &
1069  status = 'old', &
1070  iostat = ierr )
1071 
1072  if ( ierr /= 0 ) then
1073  log_error("RD_MSTRN_setup",*) 'Input data file does not found! ', trim(mstrn_aeropara_inputfile)
1074  call prc_abort
1075  endif
1076 
1077  ! read aerosol/surface parameters for check
1078  read(fid,*) nband, nsfc, nptype, nstream, nplkord, nfitplk
1079 
1080  if ( nband /= mstrn_nband ) then
1081  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nband(given,file)=', mstrn_nband, nband
1082  call prc_abort
1083  endif
1084  if ( nsfc /= mstrn_nsfc ) then
1085  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nsfc(given,file)=', mstrn_nsfc, nsfc
1086  call prc_abort
1087  endif
1088  if ( nptype /= mstrn_nptype ) then
1089  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nptype(given,file)=', mstrn_nptype, nptype
1090  call prc_abort
1091  endif
1092  if ( nstream /= mstrn_nstream ) then
1093  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nstream(given,file)=', mstrn_nstream, nstream
1094  call prc_abort
1095  endif
1096  if ( nplkord /= mstrn_nplkord ) then
1097  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nplkord(given,file)=', mstrn_nplkord, nplkord
1098  call prc_abort
1099  endif
1100  if ( nfitplk /= mstrn_nfitplk ) then
1101  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nfitPLK(given,file)=', mstrn_nfitplk, nfitplk
1102  call prc_abort
1103  endif
1104 
1105  ! wave band boundaries
1106  read(fid,*) dummy
1107  read(fid,*) waveh(:)
1108 
1109  ! for each band
1110  do iw = 1, mstrn_nband
1111 
1112  ! fitting point for planck functions
1113  read(fid,*) dummy
1114  read(fid,*) fitplk(:,iw)
1115  ! solar insolation
1116  read(fid,*) dummy
1117  read(fid,*) fsol(iw)
1118  ! surface properties
1119  read(fid,*) dummy
1120  read(fid,*) sfc(:,iw)
1121  ! rayleigh scattering
1122  read(fid,*) dummy
1123  read(fid,*) rayleigh(iw)
1124 
1125  ! moments
1126  read(fid,*) dummy
1127  do im = 1, mstrn_nmoment
1128  ! for rayleigh scattering phase function
1129  read(fid,*) qmol(im,iw)
1130  ! for aerosol scattering phase function
1131  do iptype = 1, mstrn_nptype
1132  nradius = ptype_nradius(iptype)
1133  read(fid,*) q(1:nradius,iptype,im,iw)
1134  enddo
1135  enddo
1136 
1137  enddo
1138 
1139  close(fid)
1140 
1141  fsol_tot = 0.0_rp
1142  do iw = 1, mstrn_nband
1143  fsol_tot = fsol_tot + fsol(iw)
1144  enddo
1145 
1146  !---< radius mode & hygroscopic parameter input >---
1147 
1148  ! allocate arrays
1149  allocate( hygro_flag(mstrn_nptype) )
1150  allocate( radmode(mstrn_nptype,mstrn_nradius) )
1151 
1152  open( fid, &
1153  file = trim(mstrn_hygropara_inputfile), &
1154  form = 'formatted', &
1155  status = 'old', &
1156  iostat = ierr )
1157 
1158  if ( ierr /= 0 ) then
1159  log_error("RD_MSTRN_setup",*) 'Input data file does not found! ', trim(mstrn_hygropara_inputfile)
1160  call prc_abort
1161  endif
1162 
1163  read(fid,*) nptype
1164 
1165  if ( nptype /= mstrn_nptype ) then
1166  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nptype(given,file)=', mstrn_nptype, nptype
1167  call prc_abort
1168  endif
1169 
1170  do iptype = 1, mstrn_nptype
1171  read(fid,*) dummy
1172  read(fid,*) hygro_flag(iptype), nradius
1173 
1174  if ( nradius /= ptype_nradius(iptype) ) then
1175  log_error("RD_MSTRN_setup",*) 'Inconsistent parameter value! nradius(given,file)=', ptype_nradius(iptype), nradius
1176  call prc_abort
1177  endif
1178 
1179  read(fid,*) radmode(iptype,1:nradius)
1180  enddo
1181 
1182  close(fid)
1183 
1184  !----- report data -----
1185  log_newline
1186  log_info("RD_MSTRN_setup",'(1x,A,F12.7)') 'Baseline of total solar insolation : ', fsol_tot
1187 
1188  !---< constant parameter for main scheme >---
1189  rho_std = pstd / ( rdry * tem00 ) ! [kg/m3]
1190 
1191  m(i_sw) = 1.0_rp / sqrt(3.0_rp)
1192  w(i_sw) = 1.0_rp
1193  wbar(i_sw) = 1.0_rp
1194 
1195  m(i_lw) = 1.0_rp / 1.66_rp
1196  w(i_lw) = 1.0_rp
1197  wbar(i_lw) = 2.0_rp * m(i_lw)
1198 
1199  do im = 1, 2
1200  wmns(im) = sqrt( w(im) / m(im) )
1201  wpls(im) = sqrt( w(im) * m(im) )
1202  wscale(im) = wpls(im) / wbar(im)
1203  end do
1204 
1205  !$acc enter data &
1206  !$acc& pcopyin(wgtch, fitPLK, logfitP, logfitT, fitT) &
1207  !$acc& pcopyin(radmode, ngasabs, igasabs, radmode, ptype_nradius) &
1208  !$acc& pcopyin(fsol, q, qmol, rayleigh, acfc_pow, nch, AKD, SKD) &
1209  !$acc& pcopyin(Wmns, Wpls, Wscale, W, M)
1210 
1211  return
1212  end subroutine rd_mstrn_setup
1213 
1214  !-----------------------------------------------------------------------------
1216  subroutine rd_mstrn_dtrn3( &
1217  rd_kmax, IA, IS, IE, JA, JS, JE, &
1218  ngas, &
1219  ncfc, &
1220  naero, &
1221  hydro_str, &
1222  hydro_end, &
1223  aero_str, &
1224  aero_end, &
1225  solins, &
1226  cosSZA, &
1227  rhodz, &
1228  pres, &
1229  temp, &
1230  temph, &
1231  temp_sfc, &
1232  gas, &
1233  cfc, &
1234  aerosol_conc, &
1235  aerosol_radi, &
1236  aero2ptype, &
1237  cldfrac, &
1238  albedo_sfc, &
1239  fact_ocean, &
1240  fact_land, &
1241  fact_urban, &
1242  rflux, &
1243  rflux_sfc_dn, &
1244  tauCLD_067u, &
1245  emisCLD_105u )
1246  use scale_prc, only: &
1247  prc_abort
1248  use scale_const, only: &
1249  grav => const_grav, &
1250  pstd => const_pstd, &
1251  ppm => const_ppm
1252  implicit none
1253  integer, intent(in) :: rd_kmax
1254  integer, intent(in) :: ia, is, ie
1255  integer, intent(in) :: ja, js, je
1256 
1257  integer, intent(in) :: ngas
1258  integer, intent(in) :: ncfc
1259  integer, intent(in) :: naero
1260  integer, intent(in) :: hydro_str
1261  integer, intent(in) :: hydro_end
1262  integer, intent(in) :: aero_str
1263  integer, intent(in) :: aero_end
1264  real(rp), intent(in) :: solins (ia,ja)
1265  real(rp), intent(in) :: cossza (ia,ja)
1266  real(rp), intent(in) :: rhodz (rd_kmax ,ia,ja)
1267  real(rp), intent(in) :: pres (rd_kmax ,ia,ja)
1268  real(rp), intent(in) :: temp (rd_kmax ,ia,ja)
1269  real(rp), intent(in) :: temph (rd_kmax+1,ia,ja)
1270  real(rp), intent(in) :: temp_sfc (ia,ja)
1271  real(rp), intent(in) :: gas (rd_kmax,ia,ja,ngas )
1272  real(rp), intent(in) :: cfc (rd_kmax,ia,ja,ncfc )
1273  real(rp), intent(in) :: aerosol_conc(rd_kmax,ia,ja,naero)
1274  real(rp), intent(in) :: aerosol_radi(rd_kmax,ia,ja,naero)
1275  integer, intent(in) :: aero2ptype (naero)
1276  real(rp), intent(in) :: cldfrac (rd_kmax,ia,ja)
1277  real(rp), intent(in) :: albedo_sfc (ia,ja,n_rad_dir,n_rad_rgn)
1278  real(rp), intent(in) :: fact_ocean (ia,ja)
1279  real(rp), intent(in) :: fact_land (ia,ja)
1280  real(rp), intent(in) :: fact_urban (ia,ja)
1281  real(rp), intent(out) :: rflux (rd_kmax+1,ia,ja,2,2,mstrn_ncloud)
1282  real(rp), intent(out) :: rflux_sfc_dn(ia,ja,n_rad_dir,n_rad_rgn) ! surface downward radiation flux (direct/diffuse,IR/NIR/VIS)
1283  real(rp), intent(out) :: taucld_067u (rd_kmax,ia,ja) ! 0.67 micron cloud optical depth
1284  real(rp), intent(out) :: emiscld_105u(rd_kmax,ia,ja) ! 10.5 micron cloud emissivity
1285 
1286  ! for P-T fitting
1287  real(rp) :: dz_std (rd_kmax) ! layer thickness at 0C, 1atm [cm]
1288  real(rp) :: logp ! log10(pres)
1289  real(rp) :: logt ! log10(temp)
1290  integer :: indexp (rd_kmax) ! index for interpolation in P-fitting
1291  real(rp) :: factp (rd_kmax) ! interpolation factor in P-fitting
1292  real(rp) :: factt32(rd_kmax) ! interpolation factor in T-fitting
1293  real(rp) :: factt21(rd_kmax) ! interpolation factor in T-fitting
1294  integer :: indexr (rd_kmax,naero) ! index for interpolation in R-fitting
1295  real(rp) :: factr (rd_kmax,naero) ! interpolation factor in R-fitting
1296 
1297  ! for optical thickness by gas
1298  real(rp) :: taugas(rd_kmax,mstrn_ch_limit) ! optical thickness by gas absorption (total)
1299  real(rp) :: a1, a2, a3, factpt
1300  real(rp) :: qv, length
1301  integer :: gasno
1302 
1303  ! for optical thickness by particles
1304  real(rp) :: taupr (rd_kmax,mstrn_ncloud) ! optical thickness by Rayleigh/cloud/aerosol
1305  real(rp) :: omgpr (rd_kmax,mstrn_ncloud) ! single scattering albedo by Rayleigh/cloud/aerosol
1306  real(dp) :: optparam(rd_kmax,mstrn_nmoment,mstrn_ncloud) ! optical parameters
1307  real(rp) :: q_fit, dp_p
1308 
1309  ! for planck functions
1310  real(rp) :: bbar (rd_kmax ) ! planck functions for thermal source at the interface
1311  real(rp) :: bbarh(rd_kmax+1) ! planck functions for thermal source at the center
1312  real(rp) :: b_sfc ! planck functions for thermal source at the surface
1313  real(rp) :: wl, beta
1314 
1315  ! for two-stream
1316  real(rp) :: tau(rd_kmax, mstrn_ncloud) ! total optical thickness
1317  real(rp) :: omg(rd_kmax, mstrn_ncloud) ! single scattering albedo
1318  real(rp) :: g (rd_kmax,0:2,mstrn_ncloud) ! two-stream approximation factors
1319  ! 0: always 1
1320  ! 1: asymmetry factor
1321  ! 2: truncation factor
1322  real(rp) :: b (rd_kmax,0:2,mstrn_ncloud) ! planck expansion coefficients (zero if SW)
1323  real(rp) :: fsol_rgn ! solar insolation (zero if LW)
1324 
1325  real(rp) :: flux (rd_kmax+1,2,mstrn_ncloud) ! upward/downward flux
1326  real(rp) :: flux_direct(rd_kmax+1, mstrn_ncloud) ! downward flux (direct solar)
1327 
1328  ! for satellite simulator
1329  real(rp) :: taucld (rd_kmax) ! cloud optical thickness
1330  real(rp) :: emiscld(rd_kmax) ! cloud emissivity
1331 
1332  real(rp) :: zerosw
1333  real(rp) :: valsum
1334  integer :: chmax
1335  integer :: ip, ir, irgn, irgn_alb
1336  integer :: igas, icfc, iaero, iptype, nradius
1337  integer :: iw, ich, iplk, icloud, im
1338  integer :: k, i, j
1339  !---------------------------------------------------------------------------
1340 
1341  !$acc data &
1342  !$acc& pcopyin(solins, cosSZA, rhodz, pres, temp, temph, temp_sfc, gas, cfc) &
1343  !$acc& pcopyin(aerosol_conc, aerosol_radi, aero2ptype, cldfrac, albedo_sfc) &
1344  !$acc& pcopyin(fact_ocean, fact_land, fact_urban) &
1345  !$acc& pcopyout(rflux, rflux_sfc_dn, tauCLD_067u, emisCLD_105u)
1346 
1347  !$acc kernels
1348 
1349  !$acc loop gang
1350  !$omp parallel do default(none) &
1351  !$omp private(ip, ir, irgn, irgn_alb, iptype, nradius, valsum, chmax, A1, A2, A3, factPT, qv, length, gasno, zerosw, &
1352  !$omp q_fit, dp_P, wl, beta, &
1353  !$omp dz_std, logP, logT, indexP, factP, factT32, factT21, indexR, factR, tauGAS, &
1354  !$omp tauPR, omgPR, optparam, bbar, bbarh, b_sfc, &
1355  !$omp tau, omg, b, g, fsol_rgn, flux, flux_direct, tauCLD, emisCLD) &
1356  !$omp shared(solins, cosSZA, rhodz, pres, temp, temph, temp_sfc, gas, cfc, &
1357  !$omp aerosol_conc, aerosol_radi, aero2ptype, cldfrac, albedo_sfc, &
1358  !$omp rflux, rflux_sfc_dn, tauCLD_067u, emisCLD_105u, &
1359  !$omp GRAV, Pstd, &
1360  !$omp rd_kmax, IA, IS, IE, JA, JS, JE, ncfc, naero, hydro_str, hydro_end, aero_str, aero_end, &
1361  !$omp MSTRN_nband, RHO_std, logfitP, fitT, logfitT, fitPLK, ptype_nradius, &
1362  !$omp nch, ngasabs, igasabs, radmode, waveh, iflgb, AKD, SKD, acfc_pow, rayleigh, qmol, q, fsol, fsol_tot, wgtch)
1363  do j = js, je
1364  !$acc loop gang vector &
1365  !$acc& private(dz_std, logP, logT, indexP, factP, factT32, factT21, indexR, factR, tauGAS, &
1366  !$acc& tauPR, omgPR, optparam, bbar, bbarh, b_sfc, &
1367  !$acc& tau, omg, b, g, fsol_rgn, flux, flux_direct, tauCLD, emisCLD)
1368  do i = is, ie
1369 
1370  !$acc loop gang vector
1371  do k = 1, rd_kmax
1372  dz_std(k) = rhodz(k,i,j) / rho_std * 100.0_rp ! [cm]
1373  enddo
1374 
1375  !$acc loop gang vector
1376  do k = 1, rd_kmax
1377  logp = log10( pres(k,i,j) )
1378  logt = log10( temp(k,i,j) )
1379 
1380  indexp(k) = mstrn_nfitp
1381  !$acc loop seq
1382  do ip = mstrn_nfitp, 2, -1
1383  if( logp >= logfitp(ip) ) indexp(k) = ip
1384  enddo
1385 
1386  ip = indexp(k)
1387 
1388  factp(k) = ( logp - logfitp(ip-1) ) / ( logfitp(ip) - logfitp(ip-1) )
1389 
1390  factt32(k) = ( logt - logfitt(2) ) / ( logfitt(3) - logfitt(2) ) &
1391  * ( temp(k,i,j) - fitt(1) ) / ( fitt(3) - fitt(1) )
1392  factt21(k) = ( logt - logfitt(2) ) / ( logfitt(2) - logfitt(1) ) &
1393  * ( fitt(3) - temp(k,i,j) ) / ( fitt(3) - fitt(1) )
1394  enddo
1395 
1396  !---< interpolation of mode radius & hygroscopic parameter (R-fitting) >---
1397  do iaero = 1, naero
1398  iptype = aero2ptype(iaero)
1399  nradius = ptype_nradius(iptype)
1400 
1401  !$acc loop gang vector
1402  do k = 1, rd_kmax
1403  ! [Note] Extrapolation sometimes makes unexpected value
1404  if ( aerosol_radi(k,i,j,iaero) <= radmode(iptype,1) ) then
1405 
1406  indexr(k,iaero) = 1
1407  factr(k,iaero) = 0.0_rp
1408 
1409  elseif( aerosol_radi(k,i,j,iaero) > radmode(iptype,nradius) ) then
1410 
1411  indexr(k,iaero) = nradius
1412  factr(k,iaero) = 1.0_rp
1413 
1414  else
1415  indexr(k,iaero) = -1
1416  !$acc loop seq
1417  do ir = 1, nradius-1
1418  if ( aerosol_radi(k,i,j,iaero) <= radmode(iptype,ir+1) &
1419  .AND. aerosol_radi(k,i,j,iaero) > radmode(iptype,ir ) ) then ! interpolation
1420 
1421  indexr(k,iaero) = ir
1422  factr(k,iaero) = ( aerosol_radi(k,i,j,iaero) - radmode(iptype,ir) ) &
1423  / ( radmode(iptype,ir+1) - radmode(iptype,ir) )
1424 
1425  exit
1426  endif
1427  enddo
1428  endif
1429  enddo
1430  enddo
1431 
1432  ! initialize
1433  rflux(:,i,j,:,:,:) = 0.0_rp
1434  rflux_sfc_dn( i,j,:,: ) = 0.0_rp
1435 
1436  do iw = 1, mstrn_nband
1437  if ( iflgb(i_swlw,iw) == 0 ) then ! Long wave region
1438 
1439  irgn = i_lw
1440  irgn_alb = i_r_ir ! IR
1441 
1442  else ! Short wave region
1443 
1444  irgn = i_sw
1445  if ( waveh(iw) >= 10000.0_rp / 0.7_rp ) then
1446  irgn_alb = i_r_vis ! Visible
1447  else
1448  irgn_alb = i_r_nir ! Near-IR
1449  endif
1450 
1451  endif
1452 
1453  chmax = nch(iw)
1454 
1455  !---< interpolation of gas parameters (P-T fitting) >---
1456  do ich = 1, chmax
1457  !$acc loop gang vector
1458  do k = 1, rd_kmax
1459  taugas(k,ich) = 0.0_rp
1460  enddo
1461  enddo
1462 
1463  do k = 1, rd_kmax
1464  taucld(k) = 0.0_rp
1465  enddo
1466 
1467  !--- Gas line absorption
1468  do igas = 1, ngasabs(iw)
1469  gasno = igasabs(igas,iw)
1470 
1471  !$acc loop gang vector
1472  do k = 1, rd_kmax
1473  ip = indexp(k)
1474 
1475  length = gas(k,i,j,igasabs(igas,iw)) * ppm * dz_std(k)
1476 
1477  !$acc loop seq
1478  do ich = 1, chmax
1479  a1 = akd(ich,ip-1,1,gasno,iw) * ( 1.0_rp - factp(k) )&
1480  + akd(ich,ip ,1,gasno,iw) * ( factp(k) )
1481  a2 = akd(ich,ip-1,2,gasno,iw) * ( 1.0_rp - factp(k) )&
1482  + akd(ich,ip ,2,gasno,iw) * ( factp(k) )
1483  a3 = akd(ich,ip-1,3,gasno,iw) * ( 1.0_rp - factp(k) )&
1484  + akd(ich,ip ,3,gasno,iw) * ( factp(k) )
1485 
1486  factpt = factt32(k)*(a3-a2) + a2 + factt21(k)*(a2-a1)
1487 
1488  taugas(k,ich) = taugas(k,ich) + 10.0_rp**factpt * length
1489  enddo ! channel loop
1490  enddo
1491  enddo ! gas loop
1492 
1493  !--- Gas broad absorption
1494  if ( iflgb(i_h2o_continuum,iw) == 1 ) then
1495  !$acc loop gang vector
1496  do k = 1, rd_kmax
1497  ip = indexp(k)
1498 
1499  qv = gas(k,i,j,1) * ppm * dz_std(k)
1500  length = qv*qv / ( qv + dz_std(k) )
1501 
1502  !$acc loop seq
1503  do ich = 1, chmax
1504  a1 = skd(ich,ip-1,1,iw) * ( 1.0_rp-factp(k) )&
1505  + skd(ich,ip ,1,iw) * ( factp(k) )
1506  a2 = skd(ich,ip-1,2,iw) * ( 1.0_rp-factp(k) )&
1507  + skd(ich,ip ,2,iw) * ( factp(k) )
1508  a3 = skd(ich,ip-1,3,iw) * ( 1.0_rp-factp(k) )&
1509  + skd(ich,ip ,3,iw) * ( factp(k) )
1510 
1511  factpt = factt32(k)*(a3-a2) + a2 + factt21(k)*(a2-a1)
1512 
1513  taugas(k,ich) = taugas(k,ich) + 10.0_rp**factpt * length
1514  enddo ! channel loop
1515  enddo
1516  endif
1517 
1518  if ( iflgb(i_cfc_continuum,iw) == 1 ) then
1519  !$acc loop gang vector
1520  do k = 1, rd_kmax
1521  valsum = 0.0_rp
1522  !$acc loop seq
1523  do icfc = 1, ncfc
1524  valsum = valsum + acfc_pow(icfc,iw) * cfc(k,i,j,icfc)
1525  enddo
1526  valsum = valsum * ppm * dz_std(k)
1527 
1528  !$acc loop seq
1529  do ich = 1, chmax
1530  taugas(k,ich) = taugas(k,ich) + valsum
1531  enddo
1532  enddo
1533  endif
1534 
1535  !---< particle (Rayleigh/Cloud/Aerosol) >---
1536 
1537  ! optical thickness, phase function
1538  ! im=1: extinction coefficient
1539  ! im=2,3,4: moments of the volume scattering phase function
1540 
1541  !--- Rayleigh scattering
1542  !$acc loop gang vector
1543  do k = 1, rd_kmax
1544  dp_p = rhodz(k,i,j) * grav / pstd
1545  length = rayleigh(iw) * dp_p
1546 
1547  !$acc loop seq
1548  do im = 1, mstrn_nstream*2+2
1549  optparam(k,im,i_cloud ) = qmol(im,iw) * length
1550  optparam(k,im,i_clearsky) = qmol(im,iw) * length
1551  enddo
1552  enddo
1553 
1554  !--- Cloud scattering
1555  do iaero = hydro_str, hydro_end
1556  iptype = aero2ptype(iaero)
1557 
1558  !$acc loop gang vector
1559  do k = 1, rd_kmax
1560  ir = indexr(k,iaero)
1561 
1562  length = aerosol_conc(k,i,j,iaero) * ppm * dz_std(k)
1563 
1564  !$acc loop seq
1565  do im = 1, mstrn_nstream*2+2
1566  q_fit = q(ir ,iptype,im,iw) * ( 1.0_rp-factr(k,iaero) ) &
1567  + q(ir+1,iptype,im,iw) * ( factr(k,iaero) )
1568  q_fit = max( q_fit, 0.0_rp )
1569 
1570  optparam(k,im,i_cloud) = optparam(k,im,i_cloud) + q_fit * length
1571  enddo
1572 
1573  q_fit = q(ir ,iptype,1,iw) * ( 1.0_rp-factr(k,iaero) ) &
1574  + q(ir+1,iptype,1,iw) * ( factr(k,iaero) )
1575  q_fit = max( q_fit, 0.0_rp )
1576 
1577  taucld(k) = taucld(k) + q_fit * length
1578  enddo
1579  enddo
1580 
1581  !--- Aerosol scattering
1582  do iaero = aero_str, aero_end
1583  iptype = aero2ptype(iaero)
1584 
1585  !$acc loop gang vector
1586  do k = 1, rd_kmax
1587  ir = indexr(k,iaero)
1588 
1589  length = aerosol_conc(k,i,j,iaero) * ppm * dz_std(k)
1590 
1591  !$acc loop seq
1592  do im = 1, mstrn_nstream*2+2
1593  q_fit = q(ir ,iptype,im,iw) * ( 1.0_rp-factr(k,iaero) ) &
1594  + q(ir+1,iptype,im,iw) * ( factr(k,iaero) )
1595  q_fit = max( q_fit, 0.0_rp )
1596 
1597  optparam(k,im,i_cloud ) = optparam(k,im,i_cloud ) + q_fit * length
1598  optparam(k,im,i_clearsky) = optparam(k,im,i_clearsky) + q_fit * length
1599  enddo
1600  enddo
1601  enddo
1602 
1603  do icloud = 1, mstrn_ncloud
1604  !$acc loop gang vector
1605  do k = 1, rd_kmax
1606  taupr(k,icloud) = optparam(k,1,icloud)
1607  omgpr(k,icloud) = optparam(k,1,icloud) - optparam(k,2,icloud)
1608 
1609  !--- g
1610  zerosw = 0.5_rp - sign(0.5_rp,omgpr(k,icloud)-rd_eps)
1611 
1612  g(k,0,icloud) = 1.0_rp
1613  g(k,1,icloud) = optparam(k,3,icloud) * ( 1.0_rp-zerosw ) / ( omgpr(k,icloud)+zerosw )
1614  g(k,2,icloud) = optparam(k,4,icloud) * ( 1.0_rp-zerosw ) / ( omgpr(k,icloud)+zerosw )
1615  !g(k,1,icloud) = max( optparam(k,3,icloud) * ( 1.0_RP-zerosw ) / ( omgPR(k,icloud)+zerosw ), 0.0_RP )
1616  !g(k,2,icloud) = max( optparam(k,4,icloud) * ( 1.0_RP-zerosw ) / ( omgPR(k,icloud)+zerosw ), 0.0_RP )
1617  enddo
1618  enddo
1619 
1620  if ( waveh(iw) <= 1.493e+4_rp .AND. 1.493e+4_rp < waveh(iw+1) ) then ! 0.67 micron
1621  !$acc loop gang vector
1622  do k = 1, rd_kmax
1623  taucld_067u(k,i,j) = taucld(k) ! 0.67 micron tau for resolved clouds
1624  enddo
1625  endif
1626 
1627  if ( irgn == i_sw ) then ! solar
1628  do icloud = 1, 2
1629  !$acc loop gang vector
1630  do k = 1, rd_kmax
1631  b(k,0,icloud) = 0.0_rp
1632  b(k,1,icloud) = 0.0_rp
1633  b(k,2,icloud) = 0.0_rp
1634  enddo
1635  enddo
1636 
1637  b_sfc = 0.0_rp
1638  fsol_rgn = fsol(iw) / fsol_tot * solins(i,j)
1639 
1640  elseif( irgn == i_lw ) then ! IR
1641  !--- set planck functions
1642  wl = 10000.0_rp / sqrt( waveh(iw) * waveh(iw+1) )
1643 
1644  ! from temp at cell center
1645  !$acc loop gang vector
1646  do k = 1, rd_kmax
1647  beta = 0.0_rp
1648  !$acc loop seq
1649  do iplk = mstrn_nfitplk, 1, -1
1650  beta = beta / ( wl*temp(k,i,j) ) + fitplk(iplk,iw)
1651  enddo
1652  bbar(k) = exp(-beta) * temp(k,i,j) / (wl*wl)
1653  enddo
1654 
1655  ! from temp at cell wall
1656  !$acc loop gang vector
1657  do k = 1, rd_kmax+1
1658  beta = 0.0_rp
1659  !$acc loop seq
1660  do iplk = mstrn_nfitplk, 1, -1
1661  beta = beta / ( wl*temph(k,i,j) ) + fitplk(iplk,iw)
1662  enddo
1663  bbarh(k) = exp(-beta) * temph(k,i,j) / (wl*wl)
1664  enddo
1665 
1666  ! from temp_sfc
1667  beta = 0.0_rp
1668  !$acc loop seq
1669  do iplk = mstrn_nfitplk, 1, -1
1670  beta = beta / ( wl*temp_sfc(i,j) ) + fitplk(iplk,iw)
1671  enddo
1672 
1673  b_sfc = exp(-beta) * temp_sfc(i,j) / (wl*wl)
1674 
1675  fsol_rgn = 0.0_rp
1676  endif
1677 
1678  ! sub-channel loop
1679  do ich = 1, chmax
1680 
1681  !--- total tau & omega
1682  do icloud = 1, 2
1683  !$acc loop gang vector
1684  do k = 1, rd_kmax
1685  tau(k,icloud) = taugas(k,ich) + taupr(k,icloud)
1686  !tau(k,icloud) = max( tauGAS(k,ich) + tauPR(k,icloud), 0.0_RP )
1687 
1688  zerosw = 0.5_rp - sign( 0.5_rp, tau(k,icloud)-rd_eps ) ! if tau < EPS, zerosw = 1
1689 
1690  omg(k,icloud) = ( 1.0_rp-zerosw ) * omgpr(k,icloud) / ( tau(k,icloud)-zerosw ) &
1691  + ( zerosw ) * 1.0_rp
1692 
1693  !omg(k,icloud) = min( max( omg(k,icloud), 0.0_RP ), 1.0_RP )
1694  enddo
1695  enddo
1696 
1697  if( irgn == i_lw ) then ! IR
1698 
1699  do icloud = 1, 2
1700  !$acc loop gang vector
1701  do k = 1, rd_kmax
1702  zerosw = 0.5_rp - sign( 0.5_rp, tau(k,icloud)-rd_eps ) ! if tau < EPS, zerosw = 1
1703 
1704  b(k,0,icloud) = bbarh(k)
1705  b(k,1,icloud) = ( 1.0_rp-zerosw ) &
1706  * ( - bbarh(k+1) &
1707  + 4.0_rp * bbar(k ) &
1708  - 3.0_rp * bbarh(k ) &
1709  ) / ( tau(k,icloud)-zerosw )
1710  b(k,2,icloud) = ( 1.0_rp-zerosw ) &
1711  * ( + bbarh(k+1) &
1712  - 2.0_rp * bbar(k ) &
1713  + bbarh(k ) &
1714  ) / ( tau(k,icloud)*tau(k,icloud)-zerosw ) * 2.0_rp
1715  enddo
1716  enddo
1717 
1718  endif ! solar/IR switch
1719 
1720  ! two-stream transfer
1721  !call PROF_rapstart('RD_MSTRN_twst', 3)
1722  call rd_mstrn_two_stream( rd_kmax, & ! [IN]
1723  cossza(i,j), & ! [IN]
1724  fsol_rgn, & ! [IN]
1725  irgn, & ! [IN]
1726  tau(:,:), & ! [IN]
1727  omg(:,:), & ! [IN]
1728  g(:,:,:), & ! [IN]
1729  b(:,:,:), & ! [IN]
1730  b_sfc, & ! [IN]
1731  albedo_sfc(i,j,:,irgn_alb), & ! [IN]
1732  cldfrac(:,i,j), & ! [IN]
1733  flux(:,:,:), & ! [OUT]
1734  flux_direct(:,:), & ! [OUT]
1735  emiscld(:) ) ! [OUT]
1736  !call PROF_rapend ('RD_MSTRN_twst', 3)
1737 
1738  do icloud = 1, 2
1739  !$acc loop gang vector
1740  do k = 1, rd_kmax+1
1741  rflux(k,i,j,irgn,i_up,icloud) = rflux(k,i,j,irgn,i_up,icloud) + flux(k,i_up,icloud) * wgtch(ich,iw)
1742  rflux(k,i,j,irgn,i_dn,icloud) = rflux(k,i,j,irgn,i_dn,icloud) + flux(k,i_dn,icloud) * wgtch(ich,iw)
1743  enddo
1744  enddo
1745 
1746  rflux_sfc_dn(i,j,i_r_direct ,irgn_alb) = rflux_sfc_dn(i,j,i_r_direct ,irgn_alb) &
1747  + ( flux_direct(rd_kmax+1,i_cloud) ) * wgtch(ich,iw)
1748  rflux_sfc_dn(i,j,i_r_diffuse,irgn_alb) = rflux_sfc_dn(i,j,i_r_diffuse,irgn_alb) &
1749  + ( flux(rd_kmax+1,i_dn,i_cloud) - flux_direct(rd_kmax+1,i_cloud) ) * wgtch(ich,iw)
1750 
1751 
1752  if ( waveh(iw) <= 952.0_rp .AND. 952.0_rp < waveh(iw+1) ) then ! 10.5 micron
1753  !$acc loop gang vector
1754  do k = 1, rd_kmax
1755  emiscld_105u(k,i,j) = emiscld(k) ! 10.5 micron emissivity for resolved clouds
1756  enddo
1757  endif
1758 
1759  enddo ! ICH loop
1760  enddo ! IW loop
1761 
1762  end do
1763  end do
1764  !$acc end kernels
1765 
1766  !$acc wait
1767 
1768  !$acc end data
1769 
1770  return
1771  end subroutine rd_mstrn_dtrn3
1772 
1773  !-----------------------------------------------------------------------------
1775 !OCL SERIAL
1776  subroutine rd_mstrn_two_stream( &
1777  RD_KMAX, &
1778  cosSZA0, &
1779  fsol, &
1780  irgn, &
1781  tau, &
1782  omg, &
1783  g, &
1784  b, &
1785  b_sfc, &
1786  albedo_sfc, &
1787  cldfrac, &
1788  flux, &
1789  flux_direct, &
1790  emisCLD )
1791  !$acc routine vector
1792  use scale_const, only: &
1793  pi => const_pi, &
1794  eps => const_eps, &
1795  eps1 => const_eps1
1796  implicit none
1797 
1798  integer, intent(in) :: rd_kmax
1799  real(rp), intent(in) :: cossza0 ! cos(SZA) = mu0
1800  real(rp), intent(in) :: fsol ! solar radiation intensity
1801  integer, intent(in) :: irgn ! 1:LW 2:SW
1802  real(rp), intent(in) :: tau (rd_kmax, mstrn_ncloud) ! total optical thickness (clear-sky/cloud)
1803  real(rp), intent(in) :: omg (rd_kmax, mstrn_ncloud) ! single scattering albedo (clear-sky/cloud)
1804  real(rp), intent(in) :: g (rd_kmax,0:2,mstrn_ncloud) ! two-stream approximation factors (clear-sky/cloud)
1805  real(rp), intent(in) :: b (rd_kmax,0:2,mstrn_ncloud) ! planck expansion coefficients (clear-sky/cloud)
1806  real(rp), intent(in) :: b_sfc ! planck function at surface
1807  real(rp), intent(in) :: albedo_sfc (n_rad_dir) ! surface albedo (DIRECT/DIFFUSE)
1808  real(rp), intent(in) :: cldfrac (rd_kmax) ! cloud fraction
1809  real(rp), intent(out) :: flux (rd_kmax+1,2,mstrn_ncloud) ! upward(sfc->TOA)/downward(TOA->sfc) flux (clear-sky/cloud)
1810  real(rp), intent(out) :: flux_direct(rd_kmax+1, mstrn_ncloud) ! downward(TOA->sfc) flux, solar direct (clear-sky/cloud)
1811  real(rp), intent(out) :: emiscld (rd_kmax) ! cloud emissivity factor (cloud)
1812 
1813  ! parameters with two-stream truncation
1814  real(rp) :: tau_new ! optical thickness : two-stream truncation
1815  real(rp) :: omg_new ! single scattering albedo : two-stream truncation
1816  real(rp) :: g_new ! asymmetric factor : two-stream truncation
1817  real(rp) :: b_new0 ! planck function : two-stream truncation
1818  real(rp) :: b_new1
1819  real(rp) :: b_new2
1820  real(rp) :: c0
1821  real(rp) :: c1
1822  real(rp) :: c2
1823  real(rp) :: pmns, ppls ! Phase function : two-stream truncation
1824  real(rp) :: smns, spls ! Source function : two-stream truncation
1825 
1826  ! working
1827  real(rp) :: cossza
1828  real(rp) :: x, y ! X-, X+
1829  real(rp) :: lamda ! eigenvalue of X-, X+
1830  real(rp) :: e
1831  real(rp) :: apls_mns, bpls_mns ! A+/A-, B+/B-
1832  real(dp) :: v0mns, v0pls ! V0-, V0+
1833  real(dp) :: v1mns, v1pls ! V1-, V1+
1834  real(dp) :: dmns0, dmns1, dmns2 ! D0-, D1-, D2-
1835  real(dp) :: dpls0, dpls1, dpls2 ! D0+, D1+, D2+
1836  real(rp) :: sigmns, sigpls ! sigma-, sigma+
1837  real(rp) :: qgamma ! Q * gamma
1838  real(rp) :: zerosw, tmp
1839 
1840  ! main factors
1841  real(rp) :: tdir0(rd_kmax,mstrn_ncloud) ! transmission factor for solar direct (clear-sky/cloud)
1842  real(rp) :: r0 (rd_kmax,mstrn_ncloud) ! reflection factor (clear-sky/cloud)
1843  real(rp) :: t0 (rd_kmax,mstrn_ncloud) ! transmission factor (clear-sky/cloud)
1844  real(rp) :: em_lw(rd_kmax,mstrn_ncloud) ! thermal source (sfc->TOA) (clear-sky/cloud)
1845  real(rp) :: em_sw(rd_kmax,mstrn_ncloud) ! solar source (sfc->TOA) (clear-sky/cloud)
1846  real(rp) :: ep_lw(rd_kmax,mstrn_ncloud) ! thermal source (TOA->sfc) (clear-sky/cloud)
1847  real(rp) :: ep_sw(rd_kmax,mstrn_ncloud) ! solar source (TOA->sfc) (clear-sky/cloud)
1848 
1849  ! Averaged factors, considering cloud overwrap
1850  real(rp) :: cf (rd_kmax ) ! cloud fraction
1851  real(rp) :: tau_bar_sol(rd_kmax+1) ! solar insolation through accumulated optical thickness at each layer
1852  real(rp) :: tdir (rd_kmax+1) ! transmission factor for solar direct
1853  real(rp) :: r (rd_kmax+1) ! reflection factor
1854  real(rp) :: t (rd_kmax+1) ! transmission factor
1855  real(rp) :: em (rd_kmax+1) ! source (sfc->TOA)
1856  real(rp) :: ep (rd_kmax+1) ! source (TOA->sfc)
1857 
1858  ! Doubling-Adding
1859  real(rp) :: r12mns(rd_kmax+1) ! reflection factor in doubling method
1860  real(rp) :: r12pls(rd_kmax+1) ! reflection factor in doubling method
1861  real(rp) :: e12mns(rd_kmax+1) ! source function in doubling method
1862  real(rp) :: e12pls(rd_kmax+1) ! source function in doubling method
1863  real(rp) :: umns, upls ! flux intensity
1864 
1865  real(rp) :: em0(mstrn_ncloud)
1866  real(rp) :: factor
1867  real(rp) :: wmns_irgn, m_irgn, w_irgn, wpls_irgn, wscale_irgn
1868 
1869  integer, parameter :: i_sfc2toa = 1
1870  integer, parameter :: i_toa2sfc = 2
1871 
1872  real(rp) :: sw
1873  integer :: k, icloud
1874  !---------------------------------------------------------------------------
1875 
1876  m_irgn = m(irgn)
1877  w_irgn = w(irgn)
1878  wmns_irgn = wmns(irgn)
1879  wpls_irgn = wpls(irgn)
1880  wscale_irgn = wscale(irgn)
1881 
1882  cossza = max( cossza0, rd_cossza_min )
1883 
1884  do icloud = 1, 2
1885  !$acc loop vector
1886  do k = 1, rd_kmax
1887 
1888  !---< two-stream truncation >---
1889  tau_new = ( 1.0_rp - omg(k,icloud)*g(k,2,icloud) ) * tau(k,icloud)
1890 
1891  omg_new = ( 1.0_rp - g(k,2,icloud) ) / ( 1.0_rp - omg(k,icloud)*g(k,2,icloud) ) * omg(k,icloud)
1892  omg_new = min( omg_new, eps1 )
1893 
1894  g_new = ( g(k,1,icloud) - g(k,2,icloud) ) / ( 1.0_rp - g(k,2,icloud) )
1895 
1896 #if defined(PGI) || defined(SX)
1897  tdir0(k,icloud) = exp( -min( tau_new/cossza, 1.e+3_rp ) ) ! apply exp limiter
1898 #else
1899  tdir0(k,icloud) = exp(-tau_new/cossza)
1900 #endif
1901 
1902  factor = ( 1.0_rp - omg(k,icloud)*g(k,2,icloud) )
1903  b_new0 = b(k,0,icloud)
1904  b_new1 = b(k,1,icloud) / factor
1905  b_new2 = b(k,2,icloud) / (factor*factor)
1906  c0 = wmns_irgn * 2.0_rp * pi * ( 1.0_rp - omg_new ) * b_new0
1907  c1 = wmns_irgn * 2.0_rp * pi * ( 1.0_rp - omg_new ) * b_new1
1908  c2 = wmns_irgn * 2.0_rp * pi * ( 1.0_rp - omg_new ) * b_new2
1909 
1910  !--- P+, P-
1911  pmns = omg_new * 0.5_rp * ( 1.0_rp - 3.0_rp * g_new * m_irgn*m_irgn )
1912  ppls = omg_new * 0.5_rp * ( 1.0_rp + 3.0_rp * g_new * m_irgn*m_irgn )
1913 
1914  !--- S+, S-
1915  smns = omg_new * 0.5_rp * ( 1.0_rp - 3.0_rp * g_new * m_irgn*cossza )
1916  spls = omg_new * 0.5_rp * ( 1.0_rp + 3.0_rp * g_new * m_irgn*cossza )
1917 
1918  !---< calculate R, T, e+, e- >---
1919  sw = 0.5_rp + sign(0.5_rp,tau_new-rd_eps)
1920  tau_new = max( tau_new, rd_eps )
1921 
1922  !--- X, Y
1923  x = ( 1.0_rp - w_irgn * ( ppls - pmns ) ) / m_irgn
1924  y = ( 1.0_rp - w_irgn * ( ppls + pmns ) ) / m_irgn
1925  !X = max( ( 1.0_RP - W_irgn * ( Ppls - Pmns ) ) / M_irgn, 1.E-30 )
1926  !Y = max( ( 1.0_RP - W_irgn * ( Ppls + Pmns ) ) / M_irgn, 1.E-30 )
1927  lamda = sqrt(x*y)
1928 #if defined(PGI) || defined(SX)
1929  e = exp( -min( lamda*tau_new, 1.e+3_rp ) ) ! apply exp limiter
1930 #else
1931  e = exp(-lamda*tau_new)
1932 #endif
1933 
1934  !--- A+/A-, B+/B-
1935  apls_mns = ( x * ( 1.0_dp+e ) - lamda * ( 1.0_dp-e ) ) &
1936  / ( x * ( 1.0_dp+e ) + lamda * ( 1.0_dp-e ) )
1937  bpls_mns = ( x * ( 1.0_dp-e ) - lamda * ( 1.0_dp+e ) ) &
1938  / ( x * ( 1.0_dp-e ) + lamda * ( 1.0_dp+e ) )
1939 
1940  !--- R, T
1941  r0(k,icloud) = ( sw ) * 0.5_rp * ( apls_mns + bpls_mns ) &
1942  + ( 1.0_rp-sw ) * ( tau_new * ( pmns ) / m_irgn )
1943  t0(k,icloud) = ( sw ) * 0.5_rp * ( apls_mns - bpls_mns ) &
1944  + ( 1.0_rp-sw ) * ( 1.0_rp - tau_new * ( 1.0_rp - ppls ) / m_irgn )
1945 
1946  !--- thermal source
1947  dmns0 = c0 / y + 2.0_rp * c2 / (x*y*y) + c1 / (x*y)
1948  dpls0 = c0 / y + 2.0_rp * c2 / (x*y*y) - c1 / (x*y)
1949  dmns1 = c1 / y + 2.0_rp * c2 / (x*y)
1950  dpls1 = c1 / y - 2.0_rp * c2 / (x*y)
1951  dmns2 = c2 / y
1952  dpls2 = c2 / y
1953 
1954  v0mns = dmns0
1955  v0pls = dpls0
1956  v1mns = dmns0 + dmns1*tau_new + dmns2*tau_new*tau_new
1957  v1pls = dpls0 + dpls1*tau_new + dpls2*tau_new*tau_new
1958 
1959  em_lw(k,icloud) = ( sw ) * ( v0mns - r0(k,icloud) * v0pls - t0(k,icloud) * v1mns ) &
1960  + ( 1.0_rp-sw ) * 0.5_rp * tau_new * ( 2.0_rp*c0 + c1*tau_new + c2*tau_new*tau_new )
1961  ep_lw(k,icloud) = ( sw ) * ( v1pls - t0(k,icloud) * v0pls - r0(k,icloud) * v1mns ) &
1962  + ( 1.0_rp-sw ) * 0.5_rp * tau_new * ( 2.0_rp*c0 + c1*tau_new + c2*tau_new*tau_new )
1963 
1964  !--- solar source
1965  sigmns = wmns_irgn * ( spls - smns )
1966  sigpls = wmns_irgn * ( spls + smns )
1967 
1968  tmp = x*y*cossza-1.0/cossza
1969  zerosw = 1.0_rp - sign(1.0_rp,abs(tmp)-eps) ! if abs(tmp)<EPS then 2, otherwise 0
1970  qgamma = ( sigpls*x*cossza + sigmns ) / ( tmp + zerosw*eps )
1971 
1972  v0pls = 0.5_rp * ( ( 1.0_rp + 1.0_rp/(x*cossza) ) * qgamma + sigmns / x )
1973  v0mns = 0.5_rp * ( ( 1.0_rp - 1.0_rp/(x*cossza) ) * qgamma - sigmns / x )
1974 
1975  v1pls = v0pls * tdir0(k,icloud)
1976  v1mns = v0mns * tdir0(k,icloud)
1977 
1978  em_sw(k,icloud) = ( sw ) * ( v0mns - r0(k,icloud) * v0pls - t0(k,icloud) * v1mns ) &
1979  + ( 1.0_rp-sw ) * wmns_irgn * smns * tau_new * sqrt( tdir0(k,icloud) )
1980  ep_sw(k,icloud) = ( sw ) * ( v1pls - t0(k,icloud) * v0pls - r0(k,icloud) * v1mns ) &
1981  + ( 1.0_rp-sw ) * wmns_irgn * spls * tau_new * sqrt( tdir0(k,icloud) )
1982 
1983  enddo
1984  enddo ! cloud loop
1985 
1986  !---< consider partial cloud layer: semi-random over-wrapping >---
1987 
1988  do icloud = 1, 2
1989  if ( icloud == 1 ) then
1990  !$acc loop vector
1991  do k = 1, rd_kmax
1992  cf(k) = 0.0_rp
1993  enddo
1994  else
1995  !$acc loop vector
1996  do k = 1, rd_kmax
1997  cf(k) = cldfrac(k)
1998  enddo
1999  endif
2000 
2001  !$acc loop vector
2002  do k = 1, rd_kmax
2003  tdir(k) = ( cf(k) ) * tdir0(k,i_cloud ) &
2004  + ( 1.0_rp-cf(k) ) * tdir0(k,i_clearsky)
2005  enddo
2006 
2007  tau_bar_sol(1) = fsol ! k-recurrence
2008  !$acc loop seq
2009  do k = 2, rd_kmax+1
2010  tau_bar_sol(k) = tau_bar_sol(k-1) * tdir(k-1)
2011  enddo
2012 
2013  !$acc loop vector
2014  do k = 1, rd_kmax
2015  em(k) = ( cf(k) ) * ( em_lw(k,i_cloud ) &
2016  + em_sw(k,i_cloud ) * tau_bar_sol(k) ) &
2017  + ( 1.0_rp-cf(k) ) * ( em_lw(k,i_clearsky) &
2018  + em_sw(k,i_clearsky) * tau_bar_sol(k) )
2019 
2020  ep(k) = ( cf(k) ) * ( ep_lw(k,i_cloud ) &
2021  + ep_sw(k,i_cloud ) * tau_bar_sol(k) ) &
2022  + ( 1.0_rp-cf(k) ) * ( ep_lw(k,i_clearsky) &
2023  + ep_sw(k,i_clearsky) * tau_bar_sol(k) )
2024 
2025  flux_direct(k,icloud) = cossza * tau_bar_sol(k)
2026  enddo
2027 
2028  ! at lambert surface
2029  r(rd_kmax+1) = ( cf(rd_kmax) ) * albedo_sfc(i_r_diffuse) &
2030  + ( 1.0_rp-cf(rd_kmax) ) * albedo_sfc(i_r_diffuse)
2031  t(rd_kmax+1) = 0.0_rp
2032 
2033  flux_direct(rd_kmax+1,icloud) = cossza * tau_bar_sol(rd_kmax+1)
2034 
2035  em0(i_cloud ) = wpls_irgn * ( flux_direct(rd_kmax+1,icloud) * albedo_sfc(i_r_direct) / (w_irgn*m_irgn) &
2036  + 2.0_rp * pi * ( 1.0_rp-albedo_sfc(i_r_diffuse) ) * b_sfc )
2037  em0(i_clearsky) = wpls_irgn * ( flux_direct(rd_kmax+1,icloud) * albedo_sfc(i_r_direct) / (w_irgn*m_irgn) &
2038  + 2.0_rp * pi * ( 1.0_rp-albedo_sfc(i_r_diffuse) ) * b_sfc )
2039 
2040  em(rd_kmax+1) = ( cf(rd_kmax) ) * em0(i_cloud ) &
2041  + ( 1.0_rp-cf(rd_kmax) ) * em0(i_clearsky)
2042  ep(rd_kmax+1) = 0.0_rp
2043 
2044  !---< Adding-Doubling method >---
2045  ! [note] TOA->Surface is positive direction. "pls" means upper to lower altitude.
2046 
2047  ! adding: surface to TOA
2048 
2049  r12pls(rd_kmax+1) = r(rd_kmax+1)
2050  e12mns(rd_kmax+1) = em(rd_kmax+1)
2051 
2052  !$acc loop seq
2053  do k = rd_kmax, 1, -1
2054  r(k) = ( cf(k) ) * r0(k,i_cloud ) &
2055  + ( 1.0_rp-cf(k) ) * r0(k,i_clearsky)
2056  t(k) = ( cf(k) ) * t0(k,i_cloud ) &
2057  + ( 1.0_rp-cf(k) ) * t0(k,i_clearsky)
2058 
2059  r12pls(k) = r(k) + t(k) / ( 1.0_rp - r12pls(k+1) * r(k) ) &
2060  * ( r12pls(k+1) * t(k) )
2061  e12mns(k) = em(k) + t(k) / ( 1.0_rp - r12pls(k+1) * r(k) ) &
2062  * ( r12pls(k+1) * ep(k) + e12mns(k+1) )
2063  enddo
2064 
2065  ! adding: TOA to surface
2066 
2067  r12mns(1) = r(1)
2068  e12pls(1) = ep(1)
2069 
2070  !$acc loop seq
2071  do k = 2, rd_kmax+1
2072  r12mns(k) = r(k) + t(k) / ( 1.0_rp - r12mns(k-1) * r(k) ) &
2073  * ( r12mns(k-1) *t(k) )
2074  e12pls(k) = ep(k) + t(k) / ( 1.0_rp - r12mns(k-1) * r(k) ) &
2075  * ( r12mns(k-1)*em(k) + e12pls(k-1) )
2076  enddo
2077 
2078 
2079  !--- radiative flux at cell wall:
2080  ! [note] "d" means upper to lower altitude.
2081 
2082 
2083  ! TOA boundary
2084  upls = 0.0_rp
2085  umns = e12mns(1) + r12pls(1) * upls
2086 
2087  flux(1,i_up,icloud) = wscale_irgn * umns
2088  flux(1,i_dn,icloud) = wscale_irgn * upls + flux_direct(1,icloud)
2089 
2090  !$acc loop vector
2091  do k = 2, rd_kmax+1
2092  upls = ( e12pls(k-1) + r12mns(k-1)*e12mns(k) ) / ( 1.0_rp - r12mns(k-1)*r12pls(k) )
2093  umns = e12mns(k) + r12pls(k) * upls
2094 
2095  flux(k,i_up,icloud) = wscale_irgn * umns
2096  flux(k,i_dn,icloud) = wscale_irgn * upls + flux_direct(k,icloud)
2097  enddo
2098 
2099  if ( icloud == 2 ) then ! cloud emissivity
2100  !$acc loop vector
2101  do k = 1, rd_kmax
2102  emiscld(k) = 1.0_rp - r(k) - t(k)
2103  enddo
2104  endif
2105 
2106  enddo ! cloud loop
2107 
2108 
2109  return
2110  end subroutine rd_mstrn_two_stream
2111 
2112 end module scale_atmos_phy_rd_mstrnx
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:46
scale_cpl_sfc_index::n_rad_dir
integer, parameter, public n_rad_dir
Definition: scale_cpl_sfc_index.F90:36
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
scale_cpl_sfc_index::i_r_direct
integer, parameter, public i_r_direct
Definition: scale_cpl_sfc_index.F90:37
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_use_climatology
logical, public atmos_phy_rd_profile_use_climatology
use climatology?
Definition: scale_atmos_phy_rd_profile.F90:38
scale_cpl_sfc_index::i_r_diffuse
integer, parameter, public i_r_diffuse
Definition: scale_cpl_sfc_index.F90:38
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_aerosol::n_ae
integer, parameter, public n_ae
Definition: scale_atmos_aerosol.F90:33
scale_atmos_phy_rd_common::i_dn
integer, parameter, public i_dn
Definition: scale_atmos_phy_rd_common.F90:35
scale_atmos_phy_rd_common::i_lw
integer, parameter, public i_lw
Definition: scale_atmos_phy_rd_common.F90:37
scale_cpl_sfc_index::i_r_ir
integer, parameter, public i_r_ir
Definition: scale_cpl_sfc_index.F90:29
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_prof::prof_rapstart
subroutine, public prof_rapstart(rapname_base, level, disable_barrier)
Start raptime.
Definition: scale_prof.F90:159
scale_atmos_phy_rd_common::i_up
integer, parameter, public i_up
Definition: scale_atmos_phy_rd_common.F90:34
scale_atmos_hydrometeor
module atmosphere / hydrometeor
Definition: scale_atmos_hydrometeor.F90:12
scale_io::io_get_available_fid
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_io.F90:321
scale_atmos_grid_cartesc_real
module Atmosphere GRID CartesC Real(real space)
Definition: scale_atmos_grid_cartesC_real.F90:11
scale_atmos_grid_cartesc_real::atmos_grid_cartesc_real_basepoint_lat
real(rp), public atmos_grid_cartesc_real_basepoint_lat
position of base point in real world [rad,-pi,pi]
Definition: scale_atmos_grid_cartesC_real.F90:36
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:31
scale_atmos_hydrometeor::i_hi
integer, parameter, public i_hi
ice water cloud
Definition: scale_atmos_hydrometeor.F90:83
scale_io
module STDIO
Definition: scale_io.F90:10
scale_const::const_pstd
real(rp), public const_pstd
standard pressure [Pa]
Definition: scale_const.F90:87
scale_cpl_sfc_index::i_r_nir
integer, parameter, public i_r_nir
Definition: scale_cpl_sfc_index.F90:30
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_mdry
real(rp), public const_mdry
mass weight (dry air) [g/mol]
Definition: scale_const.F90:54
scale_const::const_eps1
real(rp), public const_eps1
small number
Definition: scale_const.F90:34
scale_const::const_ppm
real(rp), parameter, public const_ppm
parts par million
Definition: scale_const.F90:91
scale_atmos_phy_rd_mstrnx::atmos_phy_rd_mstrnx_setup
subroutine, public atmos_phy_rd_mstrnx_setup(KA, KS, KE, CZ, FZ)
Setup.
Definition: scale_atmos_phy_rd_mstrnx.F90:221
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_setup_zgrid
subroutine, public atmos_phy_rd_profile_setup_zgrid(KA, KS, KE, KMAX, KADD, toa, CZ, FZ, zh, z)
Setup vertical grid for radiation.
Definition: scale_atmos_phy_rd_profile.F90:1086
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_time
module TIME
Definition: scale_time.F90:11
scale_atmos_hydrometeor::i_hc
integer, parameter, public i_hc
liquid water cloud
Definition: scale_atmos_hydrometeor.F90:81
scale_atmos_aerosol::ae_dens
real(rp), dimension(n_ae), parameter, public ae_dens
Definition: scale_atmos_aerosol.F90:53
scale_atmos_aerosol
module atmosphere / aerosol
Definition: scale_atmos_aerosol.F90:12
scale_atmos_phy_rd_common::i_sw
integer, parameter, public i_sw
Definition: scale_atmos_phy_rd_common.F90:38
scale_const::const_tem00
real(rp), parameter, public const_tem00
temperature reference (0C) [K]
Definition: scale_const.F90:90
scale_cpl_sfc_index
module coupler / surface-atmospehre
Definition: scale_cpl_sfc_index.F90:11
scale_time::time_nowdate
integer, dimension(6), public time_nowdate
current time [YYYY MM DD HH MM SS]
Definition: scale_time.F90:66
scale_cpl_sfc_index::i_r_vis
integer, parameter, public i_r_vis
Definition: scale_cpl_sfc_index.F90:31
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_read
subroutine, public atmos_phy_rd_profile_read(kmax, ngas, ncfc, naero, real_lat, now_date, zh, z, rhodz, pres, presh, temp, temph, gas, cfc, aerosol_conc, aerosol_radi, cldfrac)
Read profile for radiation.
Definition: scale_atmos_phy_rd_profile.F90:536
scale_const::const_rdry
real(rp), public const_rdry
specific gas constant (dry air) [J/kg/K]
Definition: scale_const.F90:55
scale_atmos_phy_rd_common
module atmosphere / physics / radiation / common
Definition: scale_atmos_phy_rd_common.F90:12
scale_prof::prof_rapend
subroutine, public prof_rapend(rapname_base, level, disable_barrier)
Save raptime.
Definition: scale_prof.F90:217
scale_const::const_mvap
real(rp), public const_mvap
mass weight (water vapor) [g/mol]
Definition: scale_const.F90:62
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_setup
subroutine, public atmos_phy_rd_profile_setup
Setup.
Definition: scale_atmos_phy_rd_profile.F90:134
scale_cpl_sfc_index::n_rad_rgn
integer, parameter, public n_rad_rgn
Definition: scale_cpl_sfc_index.F90:28
scale_atmos_phy_rd_mstrnx::atmos_phy_rd_mstrnx_flux
subroutine, public atmos_phy_rd_mstrnx_flux(KA, KS, KE, IA, IS, IE, JA, JS, JE, DENS, TEMP, PRES, QV, CZ, FZ, fact_ocean, fact_land, fact_urban, temp_sfc, albedo_sfc, solins, cosSZA, CLDFRAC, MP_Re, MP_Qe, AE_Re, AE_Qe, flux_rad, flux_rad_top, flux_rad_sfc_dn, dtau_s, dem_s)
Radiation main.
Definition: scale_atmos_phy_rd_mstrnx.F90:449
scale_atmos_hydrometeor::hyd_dens
real(rp), dimension(n_hyd), public hyd_dens
Definition: scale_atmos_hydrometeor.F90:95
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_atmos_hydrometeor::n_hyd
integer, parameter, public n_hyd
Definition: scale_atmos_hydrometeor.F90:79
scale_atmos_phy_rd_mstrnx
module atmosphere / physics / radiation / mstrnX
Definition: scale_atmos_phy_rd_mstrnx.F90:15
scale_atmos_phy_rd_profile
module atmosphere / physics/ radiation / profile
Definition: scale_atmos_phy_rd_profile.F90:15
Generated by   doxygen 1.8.17