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