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