SCALE-RM
scale_atmos_phy_rd_profile.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
13 !-------------------------------------------------------------------------------
14 #include "scalelib.h"
15 module scale_atmos_phy_rd_profile
16  !-----------------------------------------------------------------------------
17  !
18  !++ used modules
19  !
20  use scale_precision
21  use scale_io
22  use scale_prof
23  !-----------------------------------------------------------------------------
24  implicit none
25  private
26  !-----------------------------------------------------------------------------
27  !
28  !++ Public procedure
29  !
30  public :: atmos_phy_rd_profile_setup
31  public :: atmos_phy_rd_profile_setup_zgrid
32  public :: atmos_phy_rd_profile_finalize
33  public :: atmos_phy_rd_profile_read
34 
35  !-----------------------------------------------------------------------------
36  !
37  !++ Public parameters & variables
38  !
39  logical, public :: atmos_phy_rd_profile_use_climatology = .true.
40 
41  !-----------------------------------------------------------------------------
42  !
43  !++ Private procedure
44  !
45  private :: profile_setup_cira86
46  private :: profile_setup_mipas2001
47  private :: readfile_mipas2001
48  private :: profile_read_climatology
49  private :: profile_read_cira86
50  private :: profile_read_mipas2001
51  private :: profile_read_user
52  private :: profile_interp
53 
54  !-----------------------------------------------------------------------------
55  !
56  !++ Private parameters & variables
57  !
58  character(len=H_LONG), private :: profile_cira86_fname = "cira.nc"
59  character(len=H_LONG), private :: profile_mipas2001_dir = "."
60  character(len=H_LONG), private :: profile_user_fname = ""
61  logical, private :: atmos_phy_rd_profile_use_h2o = .true.
62  logical, private :: atmos_phy_rd_profile_use_co2 = .true.
63  logical, private :: atmos_phy_rd_profile_use_o3 = .true.
64  logical, private :: atmos_phy_rd_profile_use_n2o = .true.
65  logical, private :: atmos_phy_rd_profile_use_co = .true.
66  logical, private :: atmos_phy_rd_profile_use_ch4 = .true.
67  logical, private :: atmos_phy_rd_profile_use_o2 = .true.
68  logical, private :: atmos_phy_rd_profile_use_cfc = .true.
69  logical, private :: debug = .false.
70 
71  integer, private :: cira_ntime
72  integer, private :: cira_nplev
73  integer, private :: cira_nlat
74  real(rp), private, allocatable :: cira_nd (:) ! [day]
75  real(rp), private, allocatable :: cira_plog(:) ! log([hPa])
76  real(rp), private, allocatable :: cira_lat (:) ! [rad]
77  real(rp), private, allocatable :: cira_temp(:,:,:) ! [K]
78  real(rp), private, allocatable :: cira_z (:,:,:) ! [km]
79 
80  real(rp), private, allocatable :: interp_temp(:) ! [K]
81  real(rp), private, allocatable :: interp_z (:) ! [km]
82 
83  integer, private, parameter :: mipas_kmax = 121
84  integer, private, parameter :: mipas_ntime = 2
85  real(rp), private :: mipas_nd (0:mipas_ntime+1) ! [day]
86  real(rp), private :: mipas_lat (5) ! [rad]
87  real(rp), private :: mipas_z (mipas_kmax,4) ! [km]
88  real(rp), private :: mipas_pres(mipas_kmax,4) ! (not used) [hPa]
89  real(rp), private :: mipas_temp(mipas_kmax,4) ! (not used) [K]
90  real(rp), private :: mipas_gas (mipas_kmax,30,4) ! [ppmv]
91 
92  integer, private, parameter :: i_tropic = 1
93  integer, private, parameter :: i_midlat = 2
94  integer, private, parameter :: i_polarsum = 3
95  integer, private, parameter :: i_polarwin = 4
96 
97  integer, private, parameter :: i_n2 = 1
98  integer, private, parameter :: i_o2 = 2
99  integer, private, parameter :: i_co2 = 3
100  integer, private, parameter :: i_o3 = 4
101  integer, private, parameter :: i_h2o = 5
102  integer, private, parameter :: i_ch4 = 6
103  integer, private, parameter :: i_n2o = 7
104  integer, private, parameter :: i_hno3 = 8
105  integer, private, parameter :: i_co = 9
106  integer, private, parameter :: i_no2 = 10
107  integer, private, parameter :: i_n2o5 = 11
108  integer, private, parameter :: i_clo = 12
109  integer, private, parameter :: i_hocl = 13
110  integer, private, parameter :: i_clono2 = 14
111  integer, private, parameter :: i_no = 15
112  integer, private, parameter :: i_hno4 = 16
113  integer, private, parameter :: i_hcn = 17
114  integer, private, parameter :: i_nh3 = 18
115  integer, private, parameter :: i_f11 = 19
116  integer, private, parameter :: i_f12 = 20
117  integer, private, parameter :: i_f14 = 21
118  integer, private, parameter :: i_f22 = 22
119  integer, private, parameter :: i_ccl4 = 23
120  integer, private, parameter :: i_cof2 = 24
121  integer, private, parameter :: i_h2o2 = 25
122  integer, private, parameter :: i_c2h2 = 26
123  integer, private, parameter :: i_c2h6 = 27
124  integer, private, parameter :: i_ocs = 28
125  integer, private, parameter :: i_so2 = 29
126  integer, private, parameter :: i_sf6 = 30
127 
128  logical, private :: report_firsttime = .true.
129 
130  !-----------------------------------------------------------------------------
131 contains
132  !-----------------------------------------------------------------------------
134  subroutine atmos_phy_rd_profile_setup
135  use scale_prc, only: &
136  prc_abort
137  implicit none
138 
139  character(len=H_LONG) :: atmos_phy_rd_profile_cira86_in_filename
140  character(len=H_LONG) :: atmos_phy_rd_profile_mipas2001_in_basename
141  character(len=H_LONG) :: atmos_phy_rd_profile_user_in_filename
142 
143  namelist / param_atmos_phy_rd_profile / &
144  atmos_phy_rd_profile_use_climatology, &
145  atmos_phy_rd_profile_cira86_in_filename, &
146  atmos_phy_rd_profile_mipas2001_in_basename, &
147  atmos_phy_rd_profile_user_in_filename, &
148  atmos_phy_rd_profile_use_h2o, &
149  atmos_phy_rd_profile_use_co2, &
150  atmos_phy_rd_profile_use_o3, &
151  atmos_phy_rd_profile_use_n2o, &
152  atmos_phy_rd_profile_use_co, &
153  atmos_phy_rd_profile_use_ch4, &
154  atmos_phy_rd_profile_use_o2, &
155  atmos_phy_rd_profile_use_cfc, &
156  debug
157 
158  integer :: ierr
159  !---------------------------------------------------------------------------
160 
161  log_newline
162  log_info("ATMOS_PHY_RD_PROFILE_setup",*) 'Setup'
163 
164  atmos_phy_rd_profile_cira86_in_filename = profile_cira86_fname
165  atmos_phy_rd_profile_mipas2001_in_basename = profile_mipas2001_dir
166  atmos_phy_rd_profile_user_in_filename = profile_user_fname
167 
168  !--- read namelist
169  rewind(io_fid_conf)
170  read(io_fid_conf,nml=param_atmos_phy_rd_profile,iostat=ierr)
171 
172  if( ierr < 0 ) then !--- missing
173  log_info("ATMOS_PHY_RD_PROFILE_setup",*) 'Not found namelist. Default used.'
174  elseif( ierr > 0 ) then !--- fatal error
175  log_error("ATMOS_PHY_RD_PROFILE_setup",*) 'Not appropriate names in namelist PARAM_ATMOS_PHY_RD_PROFILE. Check!'
176  call prc_abort
177  endif
178  log_nml(param_atmos_phy_rd_profile)
179 
180  profile_cira86_fname = atmos_phy_rd_profile_cira86_in_filename
181  profile_mipas2001_dir = atmos_phy_rd_profile_mipas2001_in_basename
182  profile_user_fname = atmos_phy_rd_profile_user_in_filename
183 
184  log_newline
185  log_info("ATMOS_PHY_RD_PROFILE_setup",*) 'Climatological profile for radiation'
186 
187  if ( atmos_phy_rd_profile_use_climatology ) then
188 
189  call profile_setup_cira86
190 
191  call profile_setup_mipas2001
192 
193  endif
194 
195  return
196  end subroutine atmos_phy_rd_profile_setup
197 
198  !-----------------------------------------------------------------------------
200  subroutine atmos_phy_rd_profile_finalize
201 
202  if ( allocated( cira_nd ) ) deallocate( cira_nd )
203  if ( allocated( cira_plog ) ) deallocate( cira_plog )
204  if ( allocated( cira_lat ) ) deallocate( cira_lat )
205  if ( allocated( cira_temp ) ) deallocate( cira_temp )
206  if ( allocated( cira_z ) ) deallocate( cira_z )
207 
208  if ( allocated( interp_temp ) ) deallocate( interp_temp )
209  if ( allocated( interp_z ) ) deallocate( interp_z )
210 
211  return
212  end subroutine atmos_phy_rd_profile_finalize
213 
214  !-----------------------------------------------------------------------------
216  subroutine profile_setup_cira86
217  use scale_prc, only: &
218  prc_abort
219  use scale_const, only: &
220  const_d2r
221  use scale_calendar, only: &
222  calendar_date2daysec
223  use scale_file, only: &
224  file_open, &
225  file_get_shape, &
226  file_read, &
227  file_close
228  implicit none
229 
230  integer :: fid ! for netCDF
231  integer :: dims(3)
232 
233  integer, allocatable :: cira_date(:,:)
234  integer :: nday
235  real(dp) :: nsec
236  real(dp) :: subsec = 0.0_dp
237  integer :: offset_year = 0
238 
239  real(rp), allocatable :: tmp1d(:)
240  real(rp), allocatable :: tmp3d(:,:,:)
241 
242  integer :: n, m, t
243  !---------------------------------------------------------------------------
244 
245  log_newline
246  log_info("PROFILE_setup_CIRA86",*) 'Read CIRA86 climatology, filename : ', trim(profile_cira86_fname)
247 
248  ! open CIRA86 datafile (netCDF)
249  call file_open( profile_cira86_fname, & ! (in)
250  fid, & ! (out)
251  aggregate = .false., & ! (in)
252  postfix = "" ) ! (in)
253 
254  call file_get_shape( fid, "ta", dims(:) )
255  cira_nlat = dims(1)
256  cira_nplev = dims(2)
257  cira_ntime = dims(3)
258 
259  !print *, "CIRA_ntime", CIRA_ntime
260  !print *, "CIRA_nplev", CIRA_nplev
261  !print *, "CIRA_nlat", CIRA_nlat
262 
263  allocate( cira_nd( 0:cira_ntime+1) )
264 
265  allocate( cira_plog(cira_nplev) )
266  allocate( cira_lat(cira_nlat ) )
267 
268  allocate( cira_temp(cira_nplev,cira_nlat,0:cira_ntime+1) )
269  allocate( cira_z(cira_nplev,cira_nlat,0:cira_ntime+1) )
270 
271  ! read time list
272  allocate( cira_date(6,0:cira_ntime+1) )
273 
274  cira_date(:, 0) = (/ 1985, 12, 15, 12, 0, 0 /)
275  cira_date(:, 1) = (/ 1986, 1, 15, 12, 0, 0 /)
276  cira_date(:, 2) = (/ 1986, 2, 15, 12, 0, 0 /)
277  cira_date(:, 3) = (/ 1986, 3, 15, 12, 0, 0 /)
278  cira_date(:, 4) = (/ 1986, 4, 15, 12, 0, 0 /)
279  cira_date(:, 5) = (/ 1986, 5, 15, 12, 0, 0 /)
280  cira_date(:, 6) = (/ 1986, 6, 15, 12, 0, 0 /)
281  cira_date(:, 7) = (/ 1986, 7, 15, 12, 0, 0 /)
282  cira_date(:, 8) = (/ 1986, 8, 15, 12, 0, 0 /)
283  cira_date(:, 9) = (/ 1986, 9, 15, 12, 0, 0 /)
284  cira_date(:,10) = (/ 1986, 10, 15, 12, 0, 0 /)
285  cira_date(:,11) = (/ 1986, 11, 15, 12, 0, 0 /)
286  cira_date(:,12) = (/ 1986, 12, 15, 12, 0, 0 /)
287  cira_date(:,13) = (/ 1987, 1, 15, 12, 0, 0 /)
288 
289  do t = 0, cira_ntime+1
290  call calendar_date2daysec( nday, nsec, & ! [OUT]
291  cira_date(:,t), subsec, & ! [IN]
292  offset_year ) ! [IN]
293 
294  cira_nd(t) = real(nday,kind=rp) + nsec / 86400.0_rp
295  enddo
296  deallocate( cira_date )
297 
298  ! read pressure level [hPa]
299  allocate( tmp1d(cira_nplev) )
300 
301  call file_read( fid, "plev", tmp1d(:) )
302  do n = 1, cira_nplev
303  cira_plog(n) = log( real(tmp1d(n),kind=rp) )
304  enddo
305  deallocate( tmp1d )
306 
307  ! read latitude bin
308  allocate( tmp1d(cira_nlat) )
309 
310  call file_read( fid, "latitude", tmp1d(:) )
311  do n = 1, cira_nlat
312  cira_lat(n) = real(tmp1d(n),kind=rp) * const_d2r ! [deg]->[rad]
313  enddo
314  deallocate( tmp1d )
315 
316  !print *, "CIRA_plog", CIRA_plog
317  !print *, "CIRA_pres", exp(CIRA_plog)
318  !print *, "CIRA_lat", CIRA_lat
319 
320  ! read temperature [K]
321  allocate( tmp3d(cira_nlat,cira_nplev,cira_ntime) )
322 
323  call file_read( fid, "ta", tmp3d(:,:,:) )
324  do m = 1, cira_nlat
325  do n = 1, cira_nplev
326  do t = 1, cira_ntime
327  cira_temp(n,m,t) = real(tmp3d(m,n,t),kind=rp)
328  enddo
329  enddo
330  enddo
331  ! cyclic condition
332  cira_temp(:,:,0 ) = cira_temp(:,:,cira_ntime)
333  cira_temp(:,:,cira_ntime+1) = cira_temp(:,:,1 )
334 
335  ! avoid missing value
336  do m = 1, cira_nlat
337  do t = 1, cira_ntime
338  do n = 2, cira_nplev
339  if( cira_temp(n,m,t) >= 999.9_rp ) cira_temp(n,m,t) = cira_temp(n-1,m,t)
340  enddo
341  enddo
342  enddo
343 
344  !print *, "CIRA_temp", CIRA_temp
345 
346  ! read geopotencial height [m]
347  call file_read( fid, "zg", tmp3d(:,:,:) )
348  do m = 1, cira_nlat
349  do n = 1, cira_nplev
350  do t = 1, cira_ntime
351  cira_z(n,m,t) = real(tmp3d(m,n,t),kind=rp) * 1.e-3_rp ! [m]->[km]
352  enddo
353  enddo
354  enddo
355  ! cyclic condition
356  cira_z(:,:,0 ) = cira_z(:,:,cira_ntime)
357  cira_z(:,:,cira_ntime+1) = cira_z(:,:,1 )
358 
359  ! avoid missing value
360  do m = 1, cira_nlat
361  do t = 1, cira_ntime
362  do n = 2, cira_nplev
363  if( cira_z(n,m,t) == 0.999_rp ) cira_z(n,m,t) = cira_z(n-1,m,t)
364  enddo
365  enddo
366  enddo
367 
368  deallocate( tmp3d )
369 
370  ! close CIRA86 datafile (netCDF)
371  call file_close( fid )
372 
373  allocate( interp_temp(cira_nplev) )
374  allocate( interp_z(cira_nplev) )
375 
376  return
377  end subroutine profile_setup_cira86
378 
379  !-----------------------------------------------------------------------------
381  subroutine profile_setup_mipas2001
382  use scale_prc, only: &
383  prc_abort
384  use scale_const, only: &
385  const_d2r
386  use scale_calendar, only: &
387  calendar_date2daysec
388  implicit none
389 
390  character(len=H_LONG) :: fname
391 
392  character(len=7), parameter :: mipas_fname(4) = (/"equ.atm","day.atm","sum.atm","win.atm"/)
393 
394  integer :: mipas_date(6,0:mipas_ntime+1)
395  integer :: nday
396  real(dp) :: nsec
397  real(dp) :: subsec = 0.0_dp
398  integer :: offset_year = 0
399 
400  character(len=H_LONG) :: dummy
401  integer :: fid, ierr
402  integer :: t, l, rgn
403  !---------------------------------------------------------------------------
404 
405  log_newline
406  log_info("PROFILE_setup_MIPAS2001",*) 'Read MIPAS2001 climatology '
407 
408  mipas_date(:, 0) = (/ 2000, 12, 22, 12, 0, 0 /)
409  mipas_date(:, 1) = (/ 2001, 6, 21, 12, 0, 0 /)
410  mipas_date(:, 2) = (/ 2001, 12, 22, 12, 0, 0 /)
411  mipas_date(:, 3) = (/ 2002, 6, 21, 12, 0, 0 /)
412 
413  do t = 0, mipas_ntime+1
414  call calendar_date2daysec( nday, nsec, & ! [OUT]
415  mipas_date(:,t), subsec, & ! [IN]
416  offset_year ) ! [IN]
417 
418  mipas_nd(t) = real(nday,kind=rp) + nsec / 86400.0_rp
419  enddo
420 
421  mipas_lat(1) = -90.0_rp * const_d2r
422  mipas_lat(2) = -45.0_rp * const_d2r
423  mipas_lat(3) = 0.0_rp * const_d2r
424  mipas_lat(4) = 45.0_rp * const_d2r
425  mipas_lat(5) = 90.0_rp * const_d2r
426 
427  do rgn = i_tropic, i_polarwin
428  call io_get_fname(fname, trim(profile_mipas2001_dir)//'/'//mipas_fname(rgn))
429  log_info("PROFILE_setup_MIPAS2001",*) 'filename : ', trim(fname)
430 
431  fid = io_get_available_fid()
432  open( unit = fid, &
433  file = fname, &
434  form = 'formatted', &
435  status = 'old', &
436  iostat = ierr )
437 
438  if ( ierr /= 0 ) then !--- missing
439  log_error("PROFILE_setup_MIPAS2001",*) 'File not found. check!'
440  call prc_abort
441  endif
442 
443  do l = 1, 24
444  read(fid,*) dummy
445  enddo
446 
447  call readfile_mipas2001( fid, mipas_z(:,rgn) )
448  call readfile_mipas2001( fid, mipas_pres(:,rgn) )
449  call readfile_mipas2001( fid, mipas_temp(:,rgn) )
450 
451  call readfile_mipas2001( fid, mipas_gas(:,i_n2 ,rgn) )
452  call readfile_mipas2001( fid, mipas_gas(:,i_o2 ,rgn) )
453  call readfile_mipas2001( fid, mipas_gas(:,i_co2 ,rgn) )
454  call readfile_mipas2001( fid, mipas_gas(:,i_o3 ,rgn) )
455  call readfile_mipas2001( fid, mipas_gas(:,i_h2o ,rgn) )
456  call readfile_mipas2001( fid, mipas_gas(:,i_ch4 ,rgn) )
457  call readfile_mipas2001( fid, mipas_gas(:,i_n2o ,rgn) )
458  call readfile_mipas2001( fid, mipas_gas(:,i_hno3 ,rgn) )
459  call readfile_mipas2001( fid, mipas_gas(:,i_co ,rgn) )
460  call readfile_mipas2001( fid, mipas_gas(:,i_no2 ,rgn) )
461  call readfile_mipas2001( fid, mipas_gas(:,i_n2o5 ,rgn) )
462  call readfile_mipas2001( fid, mipas_gas(:,i_clo ,rgn) )
463  call readfile_mipas2001( fid, mipas_gas(:,i_hocl ,rgn) )
464  call readfile_mipas2001( fid, mipas_gas(:,i_clono2,rgn) )
465  call readfile_mipas2001( fid, mipas_gas(:,i_no ,rgn) )
466  call readfile_mipas2001( fid, mipas_gas(:,i_hno4 ,rgn) )
467  call readfile_mipas2001( fid, mipas_gas(:,i_hcn ,rgn) )
468  call readfile_mipas2001( fid, mipas_gas(:,i_nh3 ,rgn) )
469  call readfile_mipas2001( fid, mipas_gas(:,i_f11 ,rgn) )
470  call readfile_mipas2001( fid, mipas_gas(:,i_f12 ,rgn) )
471  call readfile_mipas2001( fid, mipas_gas(:,i_f14 ,rgn) )
472  call readfile_mipas2001( fid, mipas_gas(:,i_f22 ,rgn) )
473  call readfile_mipas2001( fid, mipas_gas(:,i_ccl4 ,rgn) )
474  call readfile_mipas2001( fid, mipas_gas(:,i_cof2 ,rgn) )
475  call readfile_mipas2001( fid, mipas_gas(:,i_h2o2 ,rgn) )
476  call readfile_mipas2001( fid, mipas_gas(:,i_c2h2 ,rgn) )
477  call readfile_mipas2001( fid, mipas_gas(:,i_c2h6 ,rgn) )
478  call readfile_mipas2001( fid, mipas_gas(:,i_ocs ,rgn) )
479  call readfile_mipas2001( fid, mipas_gas(:,i_so2 ,rgn) )
480  call readfile_mipas2001( fid, mipas_gas(:,i_sf6 ,rgn) )
481  close(fid)
482  enddo
483 
484  return
485  end subroutine profile_setup_mipas2001
486 
487  !-----------------------------------------------------------------------------
489  subroutine readfile_mipas2001( &
490  fid, &
491  var )
492  implicit none
493 
494  integer, intent(in) :: fid
495  real(rp), intent(out) :: var(121)
496 
497  character(len=H_LONG) :: dummy
498  real(rp) :: tmp5(5), tmp1
499 
500  integer :: nstr, l
501  !---------------------------------------------------------------------------
502 
503  read(fid,*) dummy
504  !LOG_INFO("readfile_MIPAS2001",*) dummy
505 
506  nstr = mipas_kmax
507  do l = 1, 24
508  read(fid,*) tmp5(:)
509  !LOG_INFO("readfile_MIPAS2001",*) l, tmp5(:)
510  var(nstr ) = tmp5(1)
511  var(nstr-1) = tmp5(2)
512  var(nstr-2) = tmp5(3)
513  var(nstr-3) = tmp5(4)
514  var(nstr-4) = tmp5(5)
515 
516  nstr = nstr - 5
517  enddo
518 
519  read(fid,*) tmp1
520  var(nstr) = tmp1
521 
522  return
523  end subroutine readfile_mipas2001
524 
525  !-----------------------------------------------------------------------------
527  subroutine atmos_phy_rd_profile_read( &
528  kmax, &
529  ngas, &
530  ncfc, &
531  naero, &
532  real_lat, &
533  now_date, &
534  zh, &
535  z, &
536  rhodz, &
537  pres, &
538  presh, &
539  temp, &
540  temph, &
541  gas, &
542  cfc, &
543  aerosol_conc, &
544  aerosol_radi, &
545  cldfrac )
546  use scale_const, only: &
547  grav => const_grav
548  use scale_atmos_solarins, only: &
549  solarins_fixedlatlon => atmos_solarins_fixedlatlon, &
550  solarins_fixeddate => atmos_solarins_fixeddate, &
551  solarins_lat => atmos_solarins_lat, &
552  solarins_date => atmos_solarins_date
553  implicit none
554 
555  integer, intent(in) :: kmax
556  integer, intent(in) :: ngas
557  integer, intent(in) :: ncfc
558  integer, intent(in) :: naero
559  real(rp), intent(in) :: real_lat
560  integer, intent(in) :: now_date(6)
561  real(rp), intent(in) :: zh(kmax+1)
562  real(rp), intent(in) :: z (kmax)
563  real(rp), intent(out) :: rhodz (kmax)
564  real(rp), intent(out) :: pres (kmax)
565  real(rp), intent(out) :: presh (kmax+1)
566  real(rp), intent(out) :: temp (kmax)
567  real(rp), intent(out) :: temph (kmax+1)
568  real(rp), intent(out) :: gas (kmax,ngas)
569  real(rp), intent(out) :: cfc (kmax,ncfc)
570  real(rp), intent(out) :: aerosol_conc(kmax,naero)
571  real(rp), intent(out) :: aerosol_radi(kmax,naero)
572  real(rp), intent(out) :: cldfrac (kmax)
573 
574  real(rp) :: lat
575  integer :: date(6)
576 
577  integer :: k
578  !---------------------------------------------------------------------------
579 
580  !$acc data copyout(rhodz,pres,presh,temp,temph,gas,cfc,aerosol_conc,aerosol_radi,cldfrac)
581 
582  if ( atmos_phy_rd_profile_use_climatology ) then
583 
584  if ( solarins_fixedlatlon ) then
585  lat = solarins_lat
586  else
587  lat = real_lat
588  endif
589 
590  date = now_date
591  if ( solarins_fixeddate ) then
592  if( solarins_date(1) >= 0 ) date(1) = solarins_date(1)
593  if( solarins_date(2) >= 1 ) date(2) = solarins_date(2)
594  if( solarins_date(3) >= 1 ) date(3) = solarins_date(3)
595  if( solarins_date(4) >= 0 ) date(4) = solarins_date(4)
596  if( solarins_date(5) >= 0 ) date(5) = solarins_date(5)
597  if( solarins_date(6) >= 0 ) date(6) = solarins_date(6)
598  endif
599 
600  call profile_read_climatology( kmax, & ! [IN]
601  ngas, & ! [IN]
602  ncfc, & ! [IN]
603  naero, & ! [IN]
604  lat, & ! [IN], tentative treatment
605  date, & ! [IN]
606  zh, & ! [IN]
607  z, & ! [IN]
608  pres, & ! [OUT]
609  presh, & ! [OUT]
610  temp, & ! [OUT]
611  temph, & ! [OUT]
612  gas, & ! [OUT]
613  cfc ) ! [OUT]
614 
615  else
616 
617  call profile_read_user( kmax, & ! [IN]
618  ngas, & ! [IN]
619  ncfc, & ! [IN]
620  naero, & ! [IN]
621  zh, & ! [IN]
622  z, & ! [IN]
623  pres, & ! [OUT]
624  presh, & ! [OUT]
625  temp, & ! [OUT]
626  temph, & ! [OUT]
627  gas, & ! [OUT]
628  cfc ) ! [OUT]
629 
630  endif
631 
632 
633  do k = 1, kmax
634  rhodz(k) = ( presh(k+1) - presh(k) ) * 100.0_rp / grav
635  enddo
636 
637  ! no cloud/aerosol
638  aerosol_conc(:,:) = 0.0_rp
639  aerosol_radi(:,:) = 0.0_rp
640  cldfrac(:) = 0.0_rp
641 
642  if ( .NOT. atmos_phy_rd_profile_use_h2o ) gas(:,1) = 0.0_rp
643  if ( .NOT. atmos_phy_rd_profile_use_co2 ) gas(:,2) = 0.0_rp
644  if ( .NOT. atmos_phy_rd_profile_use_o3 ) gas(:,3) = 0.0_rp
645  if ( .NOT. atmos_phy_rd_profile_use_n2o ) gas(:,4) = 0.0_rp
646  if ( .NOT. atmos_phy_rd_profile_use_co ) gas(:,5) = 0.0_rp
647  if ( .NOT. atmos_phy_rd_profile_use_ch4 ) gas(:,6) = 0.0_rp
648  if ( .NOT. atmos_phy_rd_profile_use_o2 ) gas(:,7) = 0.0_rp
649  if ( .NOT. atmos_phy_rd_profile_use_cfc ) cfc(:,:) = 0.0_rp
650 
651  !----- report data -----
652  if ( debug .AND. report_firsttime ) then
653  report_firsttime = .false.
654 
655  log_newline
656  log_info("ATMOS_PHY_RD_PROFILE_read",'(1x,A)') 'Vertical Coordinate'
657  log_info_cont('(1x,A)') '| -GRID CENTER- -GRID INTERFACE- |'
658  log_info_cont('(1x,A)') '| k z pres temp zh pres temp k |'
659  log_info_cont('(1x,A)') '| [km] [hPa] [K] [km] [hPa] [K] |'
660  k = 1
661  log_info_cont('(1x,A,F8.3,F10.4,F8.2,I5,A)') '| ',zh(k),presh(k),temph(k),k, ' | TOA'
662  log_info_cont('(1x,A,I5,F8.3,F10.4,F8.2,A)') '|',k,z(k),pres(k),temp(k), ' | '
663  do k = 2, kmax-1
664  log_info_cont('(1x,A,F8.3,F10.4,F8.2,I5,A)') '| ',zh(k),presh(k),temph(k),k, ' | '
665  log_info_cont('(1x,A,I5,F8.3,F10.4,F8.2,A)') '|',k,z(k),pres(k),temp(k), ' | '
666  enddo
667  k = kmax
668  log_info_cont('(1x,A,F8.3,F10.4,F8.2,I5,A)') '| ',zh(k),presh(k),temph(k),k, ' | '
669  log_info_cont('(1x,A,I5,F8.3,F10.4,F8.2,A)') '|',k,z(k),pres(k),temp(k), ' | '
670  k = kmax+1
671  log_info_cont('(1x,A,F8.3,F10.4,F8.2,I5,A)') '| ',zh(k),presh(k),temph(k),k, ' | Ground'
672  log_info_cont('(1x,A)') '|================================================================|'
673 
674  log_newline
675  log_info_cont('(1x,A)') '|=====================================================================================|'
676  log_info_cont('(1x,A)') '| -Gas concetrations [ppmv]- |'
677  log_info_cont('(1x,A)') '| k z H2O CO2 O3 N2O CO CH4 O2 |'
678  do k = 1, kmax
679  log_info_cont('(1x,A,I5,1F9.3,7ES10.3,A)') '|',k,z(k),gas(k,:),' | '
680  enddo
681  log_info_cont('(1x,A)') '|=====================================================================================|'
682 
683  endif
684 
685  !$acc update device(rhodz,pres,presh,temp,temph,gas,cfc,aerosol_conc,aerosol_radi,cldfrac)
686 
687  !$acc end data
688 
689  return
690  end subroutine atmos_phy_rd_profile_read
691 
692  !-----------------------------------------------------------------------------
694  subroutine profile_read_climatology( &
695  kmax, &
696  ngas, &
697  ncfc, &
698  naero, &
699  lat, &
700  now_date, &
701  zh, &
702  z, &
703  pres, &
704  presh, &
705  temp, &
706  temph, &
707  gas, &
708  cfc )
709  implicit none
710 
711  integer, intent(in) :: kmax
712  integer, intent(in) :: ngas
713  integer, intent(in) :: ncfc
714  integer, intent(in) :: naero
715  real(rp), intent(in) :: lat
716  integer, intent(in) :: now_date(6)
717  real(rp), intent(in) :: zh (kmax+1)
718  real(rp), intent(in) :: z (kmax)
719  real(rp), intent(out) :: pres (kmax)
720  real(rp), intent(out) :: presh(kmax+1)
721  real(rp), intent(out) :: temp (kmax)
722  real(rp), intent(out) :: temph(kmax+1)
723  real(rp), intent(out) :: gas (kmax,ngas)
724  real(rp), intent(out) :: cfc (kmax,ncfc)
725  !---------------------------------------------------------------------------
726 
727  log_info("PROFILE_read_climatology",*) 'Update climatological profile for radiation'
728 
729  call profile_read_cira86( kmax, & ! [IN]
730  lat, & ! [IN]
731  now_date(:), & ! [IN]
732  zh(:), & ! [IN]
733  z(:), & ! [IN]
734  presh(:), & ! [OUT]
735  temph(:), & ! [OUT]
736  pres(:), & ! [OUT]
737  temp(:) ) ! [OUT]
738 
739  call profile_read_mipas2001( kmax, & ! [IN]
740  ngas, & ! [IN]
741  ncfc, & ! [IN]
742  lat, & ! [IN]
743  now_date(:), & ! [IN]
744  z(:), & ! [IN]
745  gas(:,:), & ! [OUT]
746  cfc(:,:) ) ! [OUT]
747 
748  return
749  end subroutine profile_read_climatology
750 
751  !-----------------------------------------------------------------------------
753  subroutine profile_read_cira86( &
754  kmax, &
755  lat, &
756  now_date, &
757  zh, &
758  z, &
759  presh, &
760  temph, &
761  pres, &
762  temp )
763  use scale_calendar, only: &
764  calendar_date2daysec
765  implicit none
766 
767  integer, intent(in) :: kmax
768  real(rp), intent(in) :: lat
769  integer, intent(in) :: now_date(6)
770  real(rp), intent(in) :: zh (kmax+1)
771  real(rp), intent(in) :: z (kmax)
772  real(rp), intent(out) :: presh(kmax+1)
773  real(rp), intent(out) :: temph(kmax+1)
774  real(rp), intent(out) :: pres (kmax)
775  real(rp), intent(out) :: temp (kmax)
776 
777  real(rp) :: plogh(kmax+1)
778  real(rp) :: plog (kmax)
779 
780  integer :: now_date_mod(6), nday
781  real(rp) :: nd
782  real(dp) :: nsec
783  real(dp) :: subsec = 0.0_dp
784  integer :: offset_year = 0
785 
786  integer :: nplev_mod
787  integer :: indexlat, indexd
788  real(rp) :: factlat, factd
789 
790  integer :: n, t
791  !---------------------------------------------------------------------------
792 
793  ! latitude interpolation
794  if ( lat < cira_lat(1) ) then ! extrapolation
795  indexlat = 1
796  factlat = 0.0_rp
797  elseif( lat >= cira_lat(cira_nlat) ) then ! extrapolation
798  indexlat = cira_nlat - 1
799  factlat = 1.0_rp
800  else
801  do n = 1, cira_nlat-1
802  if ( lat >= cira_lat(n ) &
803  .AND. lat < cira_lat(n+1) ) then ! interpolation
804  indexlat = n
805  factlat = ( lat-cira_lat(n) ) / ( cira_lat(n+1)-cira_lat(n) )
806  endif
807  enddo
808  endif
809 
810  ! time interpolation
811  now_date_mod(2:6) = now_date(2:6)
812  now_date_mod(1) = 1986
813 
814  call calendar_date2daysec( nday, nsec, & ! [OUT]
815  now_date_mod(:), subsec, & ! [IN]
816  offset_year ) ! [IN]
817 
818  nd = real(nday,kind=rp) + nsec / 86400.0_rp
819 
820  do t = 0, cira_ntime
821  if ( nd >= cira_nd(t ) &
822  .AND. nd < cira_nd(t+1) ) then ! interpolation
823  indexd = t
824  factd = ( nd-cira_nd(t) ) / ( cira_nd(t+1)-cira_nd(t) )
825  endif
826  enddo
827 
828  interp_z(:) = cira_z(:,indexlat ,indexd ) * ( 1.0_rp-factlat ) * ( 1.0_rp-factd ) &
829  + cira_z(:,indexlat+1,indexd ) * ( factlat ) * ( 1.0_rp-factd ) &
830  + cira_z(:,indexlat ,indexd+1) * ( 1.0_rp-factlat ) * ( factd ) &
831  + cira_z(:,indexlat+1,indexd+1) * ( factlat ) * ( factd )
832 
833  interp_temp(:) = cira_temp(:,indexlat ,indexd ) * ( 1.0_rp-factlat ) * ( 1.0_rp-factd ) &
834  + cira_temp(:,indexlat+1,indexd ) * ( factlat ) * ( 1.0_rp-factd ) &
835  + cira_temp(:,indexlat ,indexd+1) * ( 1.0_rp-factlat ) * ( factd ) &
836  + cira_temp(:,indexlat+1,indexd+1) * ( factlat ) * ( factd )
837 
838  ! avoid missing value
839  nplev_mod = cira_nplev
840  do n = cira_nplev, 1, -1
841  if ( interp_temp(n) == interp_temp(n-1) ) then
842  nplev_mod = nplev_mod-1
843  else
844  exit
845  endif
846  enddo
847 
848  !do n = 1, nplev_mod
849  ! print *, n, interp_z(n), interp_temp(n), indexLAT, indexD, factLAT, factD
850  !enddo
851 
852  ! pressure interpolation
853  call profile_interp( nplev_mod, & ! [IN]
854  interp_z(1:nplev_mod), & ! [IN]
855  cira_plog(1:nplev_mod), & ! [IN]
856  kmax+1, & ! [IN]
857  zh(:), & ! [IN]
858  plogh(:) ) ! [OUT]
859 
860  presh(:) = exp( plogh(:) )
861 
862  call profile_interp( kmax+1, zh(:), plogh(:), kmax, z(:), plog(:) )
863  pres(:) = exp( plog(:) )
864 
865  call profile_interp( nplev_mod, & ! [IN]
866  interp_z(1:nplev_mod), & ! [IN]
867  interp_temp(1:nplev_mod), & ! [IN]
868  kmax, & ! [IN]
869  z(:), & ! [IN]
870  temp(:) ) ! [OUT]
871 
872  call profile_interp( nplev_mod, & ! [IN]
873  interp_z(1:nplev_mod), & ! [IN]
874  interp_temp(1:nplev_mod), & ! [IN]
875  kmax+1, & ! [IN]
876  zh(:), & ! [IN]
877  temph(:) ) ! [OUT]
878 
879  return
880  end subroutine profile_read_cira86
881 
882  !-----------------------------------------------------------------------------
884  subroutine profile_read_mipas2001( &
885  kmax, &
886  ngas, &
887  ncfc, &
888  lat, &
889  now_date, &
890  z, &
891  gas, &
892  cfc )
893  use scale_calendar, only: &
894  calendar_date2daysec
895  implicit none
896 
897  integer, intent(in) :: kmax
898  integer, intent(in) :: ngas
899  integer, intent(in) :: ncfc
900  real(rp), intent(in) :: lat
901  integer, intent(in) :: now_date(6)
902  real(rp), intent(in) :: z (kmax)
903  real(rp), intent(inout) :: gas(kmax,ngas)
904  real(rp), intent(inout) :: cfc(kmax,ncfc)
905 
906  real(rp) :: interp_gas(mipas_kmax,30)
907  real(rp) :: interp_z (mipas_kmax)
908 
909  integer :: now_date_mod(6), nday
910  real(rp) :: nd
911  real(dp) :: nsec
912  real(dp) :: subsec = 0.0_dp
913  integer :: offset_year = 0
914 
915  integer :: indexd1, indexd2
916  real(rp) :: factlat, factd
917  !---------------------------------------------------------------------------
918 
919  ! time interpolation
920  now_date_mod(2:6) = now_date(2:6)
921  now_date_mod(1) = 2001
922 
923  call calendar_date2daysec( nday, nsec, & ! [OUT]
924  now_date_mod(:), subsec, & ! [IN]
925  offset_year ) ! [IN]
926 
927  nd = real(nday,kind=rp) + nsec / 86400.0_rp
928 
929  if ( nd >= mipas_nd(0) .AND. nd < mipas_nd(1) ) then ! winter-summer
930 
931  indexd1 = i_polarwin
932  indexd2 = i_polarsum
933  factd = ( nd-mipas_nd(0) ) / ( mipas_nd(1)-mipas_nd(0) )
934 
935  elseif( nd >= mipas_nd(1) .AND. nd < mipas_nd(2) ) then ! summer-winter
936 
937  indexd1 = i_polarsum
938  indexd2 = i_polarwin
939  factd = ( nd-mipas_nd(1) ) / ( mipas_nd(2)-mipas_nd(1) )
940 
941  elseif( nd >= mipas_nd(2) .AND. nd < mipas_nd(3) ) then ! winter-summer
942 
943  indexd1 = i_polarwin
944  indexd2 = i_polarsum
945  factd = ( nd-mipas_nd(2) ) / ( mipas_nd(3)-mipas_nd(2) )
946 
947  endif
948 
949  ! latitude interpolation
950  if ( lat < mipas_lat(1) ) then ! south pole
951 
952  interp_gas(:,:) = mipas_gas(:,:,indexd1 ) * ( 1.0_rp-factd ) &
953  + mipas_gas(:,:,indexd2 ) * ( factd )
954 
955  interp_z(:) = mipas_z(:,indexd1 ) * ( 1.0_rp-factd ) &
956  + mipas_z(:,indexd2 ) * ( factd )
957 
958  elseif( lat >= mipas_lat(1) .AND. lat < mipas_lat(2) ) then ! south pole-SH mid
959 
960  factlat = ( lat-mipas_lat(1) ) / ( mipas_lat(2)-mipas_lat(1) )
961 
962  interp_gas(:,:) = mipas_gas(:,:,indexd1) * ( 1.0_rp-factd ) * ( 1.0_rp-factlat ) &
963  + mipas_gas(:,:,indexd2) * ( factd ) * ( 1.0_rp-factlat ) &
964  + mipas_gas(:,:,i_midlat) * ( factlat )
965 
966  interp_z(:) = mipas_z(:,indexd1) * ( 1.0_rp-factd ) * ( 1.0_rp-factlat ) &
967  + mipas_z(:,indexd2) * ( factd ) * ( 1.0_rp-factlat ) &
968  + mipas_z(:,i_midlat) * ( factlat )
969 
970  elseif( lat >= mipas_lat(2) .AND. lat < mipas_lat(3) ) then ! SH mid-EQ
971 
972  factlat = ( lat-mipas_lat(2) ) / ( mipas_lat(3)-mipas_lat(2) )
973 
974  interp_gas(:,:) = mipas_gas(:,:,i_midlat) * ( 1.0_rp-factlat ) &
975  + mipas_gas(:,:,i_tropic) * ( factlat )
976 
977  interp_z(:) = mipas_z(:,i_midlat) * ( 1.0_rp-factlat ) &
978  + mipas_z(:,i_tropic) * ( factlat )
979 
980  elseif( lat >= mipas_lat(3) .AND. lat < mipas_lat(4) ) then ! EQ-NH mid
981 
982  factlat = ( lat-mipas_lat(3) ) / ( mipas_lat(4)-mipas_lat(3) )
983 
984  interp_gas(:,:) = mipas_gas(:,:,i_tropic) * ( 1.0_rp-factlat ) &
985  + mipas_gas(:,:,i_midlat) * ( factlat )
986 
987  interp_z(:) = mipas_z(:,i_tropic) * ( 1.0_rp-factlat ) &
988  + mipas_z(:,i_midlat) * ( factlat )
989 
990  elseif( lat >= mipas_lat(4) .AND. lat < mipas_lat(5) ) then ! NH mid-north pole
991 
992  factlat = ( lat-mipas_lat(4) ) / ( mipas_lat(5)-mipas_lat(4) )
993 
994  interp_gas(:,:) = mipas_gas(:,:,i_midlat) * ( 1.0_rp-factlat ) &
995  + mipas_gas(:,:,indexd2) * ( 1.0_rp-factd ) * ( factlat ) &
996  + mipas_gas(:,:,indexd1) * ( factd ) * ( factlat )
997 
998  interp_z(:) = mipas_z(:,i_midlat) * ( 1.0_rp-factlat ) &
999  + mipas_z(:,indexd2) * ( 1.0_rp-factd ) * ( factlat ) &
1000  + mipas_z(:,indexd1) * ( factd ) * ( factlat )
1001 
1002  elseif( lat >= mipas_lat(5) ) then ! north pole
1003 
1004  interp_gas(:,:) = mipas_gas(:,:,indexd2) * ( 1.0_rp-factd ) &
1005  + mipas_gas(:,:,indexd1) * ( factd )
1006 
1007  interp_z(:) = mipas_z(:,indexd2) * ( 1.0_rp-factd ) &
1008  + mipas_z(:,indexd1) * ( factd )
1009 
1010  endif
1011 
1012  gas(:,:) = 0.0_rp
1013  cfc(:,:) = 0.0_rp
1014 
1015  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_h2o ), kmax, z(:), gas(:,1) )
1016  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_co2 ), kmax, z(:), gas(:,2) )
1017  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_o3 ), kmax, z(:), gas(:,3) )
1018  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_n2o ), kmax, z(:), gas(:,4) )
1019  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_co ), kmax, z(:), gas(:,5) )
1020  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_ch4 ), kmax, z(:), gas(:,6) )
1021  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_o2 ), kmax, z(:), gas(:,7) )
1022 
1023  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f11 ), kmax, z(:), cfc(:, 1) )
1024  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f12 ), kmax, z(:), cfc(:, 2) )
1025  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f14 ), kmax, z(:), cfc(:, 4) )
1026  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f22 ), kmax, z(:), cfc(:, 9) )
1027  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_sf6 ), kmax, z(:), cfc(:,22) )
1028  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_clono2), kmax, z(:), cfc(:,23) )
1029  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_ccl4 ), kmax, z(:), cfc(:,24) )
1030  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_n2o5 ), kmax, z(:), cfc(:,25) )
1031  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_hno4 ), kmax, z(:), cfc(:,27) )
1032 
1033  return
1034  end subroutine profile_read_mipas2001
1035 
1036  !-----------------------------------------------------------------------------
1038  subroutine profile_interp( imax1, x1, y1, imax2, x2, y2 )
1039  implicit none
1040 
1041  integer, intent(in) :: imax1
1042  real(rp), intent(in) :: x1(imax1)
1043  real(rp), intent(in) :: y1(imax1)
1044  integer, intent(in) :: imax2
1045  real(rp), intent(in) :: x2(imax2)
1046  real(rp), intent(out) :: y2(imax2)
1047 
1048  real(rp) :: fact
1049 
1050  integer :: i1, i2
1051  !---------------------------------------------------------------------------
1052 
1053  !$omp parallel do default(none) &
1054  !$omp shared(imax1,x2,x1,y2,y1,imax2) &
1055  !$omp private(i1,fact,i2) OMP_SCHEDULE_
1056  do i2 = 1, imax2
1057 
1058  if ( x2(i2) > x1(1) ) then ! extrapolation
1059 
1060  fact = ( x1(1) - x2(i2) ) / ( x1(2) - x1(1) )
1061 
1062  y2(i2) = y1(1) * ( 1.0_rp-fact ) &
1063  + y1(2) * ( fact )
1064 
1065  elseif( x2(i2) <= x1(imax1) ) then ! extrapolation
1066 
1067  fact = ( x1(imax1) - x2(i2) ) / ( x1(imax1) - x1(imax1-1) )
1068 
1069  y2(i2) = y1(imax1-1) * ( fact ) &
1070  + y1(imax1 ) * ( 1.0_rp-fact )
1071 
1072  else
1073  do i1 = 1, imax1-1
1074  if ( x2(i2) <= x1(i1 ) &
1075  .AND. x2(i2) > x1(i1+1) ) then ! interpolation
1076 
1077  fact = ( x2(i2) - x1(i1) ) / ( x1(i1+1) - x1(i1) )
1078 
1079  y2(i2) = y1(i1 ) * ( 1.0_rp-fact ) &
1080  + y1(i1+1) * ( fact )
1081 
1082  exit
1083  endif
1084  enddo
1085  endif
1086 
1087  enddo
1088 
1089  return
1090  end subroutine profile_interp
1091 
1092  !-----------------------------------------------------------------------------
1094  subroutine atmos_phy_rd_profile_setup_zgrid( &
1095  KA, KS, KE, &
1096  KMAX, KADD, &
1097  toa, CZ, FZ, &
1098  zh, z )
1099  implicit none
1100  integer, intent(in) :: ka, ks, ke
1101  integer, intent(in) :: kmax, kadd
1102  real(rp), intent(in) :: toa
1103  real(rp), intent(in) :: cz( ka)
1104  real(rp), intent(in) :: fz(0:ka)
1105  real(rp), intent(out) :: zh(kmax+1)
1106  real(rp), intent(out) :: z (kmax)
1107 
1108  real(rp) :: dz
1109  integer :: k, rd_k
1110  !---------------------------------------------------------------------------
1111 
1112  if ( kadd > 0 ) then
1113  !--- additional layer over the computational domain
1114  dz = ( toa - fz(ke)*1.e-3_rp ) / real( kadd, kind=rp )
1115 
1116  zh(1) = toa
1117  do k = 2, kadd
1118  zh(k) = zh(k-1) - dz
1119  enddo
1120  zh(kadd+1) = fz(ke)*1.e-3_rp
1121 
1122  ! linear interpolation for center of layers
1123  do k = 1, kadd
1124  z(k) = 0.5_rp * ( zh(k+1) + zh(k) )
1125  enddo
1126  endif
1127 
1128  !--- in computational domain [NOTE] layer index is reversed
1129  do k = ks-1, ke
1130  rd_k = kmax - ( k - ks )
1131  zh(rd_k) = fz(k)*1.e-3_rp
1132  enddo
1133  do k = ks, ke
1134  rd_k = kmax - ( k - ks )
1135  z(rd_k) = cz(k)*1.e-3_rp
1136  enddo
1137 
1138  !----- report data -----
1139  if ( debug ) then
1140 
1141  log_newline
1142  log_info("ATMOS_PHY_RD_PROFILE_setup_zgrid",'(1x,A)') 'Vertical Coordinate'
1143  log_info_cont('(1x,A)') '| -GRID CENTER- -GRID INTERFACE- |'
1144  log_info_cont('(1x,A)') '| RD_k z k CZ FZ k zh RD_k |'
1145  if ( kadd > 0 ) then
1146  rd_k = 1
1147  log_info_cont('(1x,A,F8.3,I5,A)') '| ',zh(rd_k),rd_k, ' | TOA'
1148  log_info_cont('(1x,A,I5,F8.3,A)') '|',rd_k,z(rd_k), ' | '
1149  do rd_k = 2, kadd-1
1150  log_info_cont('(1x,A,F8.3,I5,A)') '| ',zh(rd_k),rd_k, ' | '
1151  log_info_cont('(1x,A,I5,F8.3,A)') '|',rd_k,z(rd_k), ' | '
1152  enddo
1153  rd_k = kadd
1154  log_info_cont('(1x,A,F8.3,I5,A)') '| ',zh(rd_k),rd_k, ' | '
1155  log_info_cont('(1x,A,I5,F8.3,A)') '|',rd_k,z(rd_k), ' | KADD'
1156  rd_k = kadd+1
1157  k = kmax - rd_k + ks
1158  log_info_cont('(1x,A,F8.3,I5,F8.3,I5,A)') '| ',fz(k)*1.e-3_rp,k,zh(rd_k),rd_k,' | '
1159  log_info_cont('(1x,A,I5,F8.3,I5,F8.3,A)') '|',rd_k,z(rd_k),k,cz(k)*1.e-3_rp, ' | KADD+1=KE'
1160  else
1161  rd_k = 1
1162  k = kmax - rd_k + ks
1163  log_info_cont('(1x,A,F8.3,I5,F8.3,I5,A)') '| ',fz(k)*1.e-3_rp,k,zh(rd_k),rd_k,' | TOA=KE'
1164  log_info_cont('(1x,A,I5,F8.3,I5,F8.3,A)') '|',rd_k,z(rd_k),k,cz(k)*1.e-3_rp, ' | '
1165  endif
1166  do rd_k = kadd+2, kmax-1
1167  k = kmax - rd_k + ks
1168  log_info_cont('(1x,A,F8.3,I5,F8.3,I5,A)') '| ',fz(k)*1.e-3_rp,k,zh(rd_k),rd_k,' | '
1169  log_info_cont('(1x,A,I5,F8.3,I5,F8.3,A)') '|',rd_k,z(rd_k),k,cz(k)*1.e-3_rp, ' | '
1170  enddo
1171  rd_k = kmax
1172  k = kmax - rd_k + ks
1173  log_info_cont('(1x,A,F8.3,I5,F8.3,I5,A)') '| ',fz(k),k,zh(rd_k),rd_k, ' | '
1174  log_info_cont('(1x,A,I5,F8.3,I5,F8.3,A)') '|',rd_k,z(rd_k),k,cz(k)*1.e-3_rp, ' | RD_KMAX=KS'
1175  rd_k = kmax+1
1176  k = kmax - rd_k + ks
1177  log_info_cont('(1x,A,F8.3,I5,F8.3,I5,A)') '| ',fz(k)*1.e-3_rp,k,zh(rd_k),rd_k,' | Ground'
1178  log_info_cont('(1x,A)') '|=====================================================|'
1179 
1180  endif
1181 
1182  return
1183  end subroutine atmos_phy_rd_profile_setup_zgrid
1184 
1185  !-----------------------------------------------------------------------------
1187  subroutine profile_read_user( &
1188  kmax, &
1189  ngas, &
1190  ncfc, &
1191  naero, &
1192  zh, &
1193  z, &
1194  pres, &
1195  presh, &
1196  temp, &
1197  temph, &
1198  gas, &
1199  cfc )
1200  use scale_prc, only: &
1201  prc_abort
1202  use scale_const, only: &
1203  mdry => const_mdry, &
1204  mvap => const_mvap, &
1205  ppm => const_ppm
1206  implicit none
1207 
1208  integer, intent(in) :: kmax
1209  integer, intent(in) :: ngas
1210  integer, intent(in) :: ncfc
1211  integer, intent(in) :: naero
1212  real(rp), intent(in) :: zh (kmax+1)
1213  real(rp), intent(in) :: z (kmax)
1214  real(rp), intent(out) :: pres (kmax)
1215  real(rp), intent(out) :: presh(kmax+1)
1216  real(rp), intent(out) :: temp (kmax)
1217  real(rp), intent(out) :: temph(kmax+1)
1218  real(rp), intent(out) :: gas (kmax,ngas)
1219  real(rp), intent(out) :: cfc (kmax,ncfc)
1220 
1221  integer, parameter :: user_klim = 500
1222  integer :: user_kmax
1223  real(rp) :: user_z (user_klim)
1224  real(rp) :: user_pres(user_klim)
1225  real(rp) :: user_temp(user_klim)
1226  real(rp) :: user_qv (user_klim)
1227  real(rp) :: user_o3 (user_klim)
1228 
1229  real(rp), allocatable :: work_z(:)
1230  real(rp), allocatable :: work (:)
1231 
1232  real(rp) :: plog (kmax)
1233  real(rp) :: plogh(kmax+1)
1234 
1235  character(len=H_LONG) :: dummy, fname
1236 
1237  integer :: fid, ierr
1238  integer :: k
1239  !---------------------------------------------------------------------------
1240 
1241  log_info("PROFILE_read_user",*) 'user-defined profile'
1242 
1243  gas(:,:) = 0.0_rp
1244  cfc(:,:) = 0.0_rp
1245 
1246  call io_get_fname(fname, profile_user_fname)
1247  log_info("PROFILE_read_user",*) 'FILENAME:', trim(fname)
1248 
1249  fid = io_get_available_fid()
1250  open( unit = fid, &
1251  file = fname, &
1252  form = 'formatted', &
1253  status = 'old', &
1254  iostat = ierr )
1255 
1256  if ( ierr /= 0 ) then !--- missing
1257  log_error("PROFILE_read_user",*) 'File not found. check!'
1258  call prc_abort
1259  endif
1260 
1261  read(fid,*) dummy
1262 
1263  do k = 1, user_klim
1264  read(fid,*,iostat=ierr) user_z(k), user_pres(k), user_temp(k), user_qv(k), user_o3(k)
1265  if ( ierr /= 0 ) exit
1266  enddo
1267  user_kmax = k - 1
1268  close(fid)
1269 
1270  allocate( work_z(user_kmax) )
1271  allocate( work(user_kmax) )
1272 
1273  do k = 1, user_kmax
1274  work_z(k) = user_z(k) / 1000.0_rp ! [m->km]
1275  work(k) = log( user_pres(k)/100.0_rp ) ! log[Pa->hPA]
1276  enddo
1277 
1278  ! pressure interpolation
1279  call profile_interp( user_kmax, & ! [IN]
1280  work_z(:), & ! [IN]
1281  work(:), & ! [IN]
1282  kmax+1, & ! [IN]
1283  zh(:), & ! [IN]
1284  plogh(:) ) ! [OUT]
1285 
1286  presh(:) = exp( plogh(:) )
1287 
1288  call profile_interp( kmax+1, zh(:), plogh(:), kmax, z(:), plog(:) )
1289  pres(:) = exp( plog(:) )
1290 
1291  do k = 1, user_kmax
1292  work(k) = user_temp(k)
1293  enddo
1294 
1295  call profile_interp( user_kmax, & ! [IN]
1296  work_z(:), & ! [IN]
1297  work(:), & ! [IN]
1298  kmax, & ! [IN]
1299  z(:), & ! [IN]
1300  temp(:) ) ! [OUT]
1301 
1302  call profile_interp( user_kmax, & ! [IN]
1303  work_z(:), & ! [IN]
1304  work(:), & ! [IN]
1305  kmax+1, & ! [IN]
1306  zh(:), & ! [IN]
1307  temph(:) ) ! [OUT]
1308 
1309  do k = 1, user_kmax
1310  work(k) = user_qv(k) / mvap * mdry / ppm ! [kg/kg->PPM]
1311  enddo
1312 
1313  call profile_interp( user_kmax, & ! [IN]
1314  work_z(:), & ! [IN]
1315  work(:), & ! [IN]
1316  kmax, & ! [IN]
1317  z(:), & ! [IN]
1318  gas(:,1) ) ! [OUT]
1319 
1320  do k = 1, user_kmax
1321  work(k) = user_o3(k) / 48.0_rp * mdry / ppm ! [kg/kg->PPM]
1322  enddo
1323 
1324  call profile_interp( user_kmax, & ! [IN]
1325  work_z(:), & ! [IN]
1326  work(:), & ! [IN]
1327  kmax, & ! [IN]
1328  z(:), & ! [IN]
1329  gas(:,3) ) ! [OUT]
1330 
1331  return
1332  end subroutine profile_read_user
1333 
1334 end module scale_atmos_phy_rd_profile
scale_const::const_grav
real(rp), public const_grav
standard acceleration of gravity [m/s2]
Definition: scale_const.F90:49
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:350
scale_file::file_open
subroutine, public file_open(basename, fid, mode, single, allnodes, aggregate, rankid, postfix)
Definition: scale_file.F90:536
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:39
scale_precision
module PRECISION
Definition: scale_precision.F90:14
scale_atmos_solarins::atmos_solarins_fixedlatlon
logical, public atmos_solarins_fixedlatlon
Definition: scale_atmos_solarins.F90:56
scale_io::io_get_available_fid
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_io.F90:373
scale_calendar
module CALENDAR
Definition: scale_calendar.F90:13
scale_atmos_solarins
module atmosphere / SOLARINS
Definition: scale_atmos_solarins.F90:14
scale_file
module file
Definition: scale_file.F90:15
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:41
scale_io::io_get_fname
subroutine, public io_get_fname(outstr, instr, rank, ext, len)
generate process specific filename
Definition: scale_io.F90:421
scale_io
module STDIO
Definition: scale_io.F90:10
scale_file::file_close
subroutine, public file_close(fid, abort)
Definition: scale_file.F90:6234
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_const::const_mdry
real(rp), public const_mdry
mass weight (dry air) [g/mol]
Definition: scale_const.F90:58
scale_const::const_ppm
real(rp), parameter, public const_ppm
parts par million
Definition: scale_const.F90:100
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_setup_zgrid
subroutine, public atmos_phy_rd_profile_setup_zgrid(KA, KS, KE, KMAX, KADD, toa, CZ, FZ, zh, z)
Setup vertical grid for radiation.
Definition: scale_atmos_phy_rd_profile.F90:1099
scale_calendar::calendar_date2daysec
subroutine, public calendar_date2daysec(absday, abssec, ymdhms, subsec, offset_year)
Convert from gregorian date to absolute day/second.
Definition: scale_calendar.F90:192
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_solarins::atmos_solarins_lat
real(rp), public atmos_solarins_lat
Definition: scale_atmos_solarins.F90:58
scale_precision::dp
integer, parameter, public dp
Definition: scale_precision.F90:32
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_finalize
subroutine, public atmos_phy_rd_profile_finalize
finalize
Definition: scale_atmos_phy_rd_profile.F90:201
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:546
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:33
scale_const::const_mvap
real(rp), public const_mvap
mass weight (water vapor) [g/mol]
Definition: scale_const.F90:67
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:135
scale_atmos_solarins::atmos_solarins_date
integer, dimension(6), public atmos_solarins_date
Definition: scale_atmos_solarins.F90:61
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
scale_atmos_phy_rd_profile
module atmosphere / physics/ radiation / profile
Definition: scale_atmos_phy_rd_profile.F90:15
scale_atmos_solarins::atmos_solarins_fixeddate
logical, public atmos_solarins_fixeddate
Definition: scale_atmos_solarins.F90:60
Generated by   doxygen 1.8.17