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