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