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