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