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