SCALE-RM
scale_atmos_phy_rd_profile.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
17 !-------------------------------------------------------------------------------
18 #include "inc_openmp.h"
19 module scale_atmos_phy_rd_profile
20  !-----------------------------------------------------------------------------
21  !
22  !++ used modules
23  !
24  use scale_precision
25  use scale_stdio
26  use scale_prof
27  use scale_grid_index
28  !-----------------------------------------------------------------------------
29  implicit none
30  private
31  !-----------------------------------------------------------------------------
32  !
33  !++ Public procedure
34  !
35  public :: atmos_phy_rd_profile_setup
36  public :: atmos_phy_rd_profile_setup_zgrid
37  public :: atmos_phy_rd_profile_read
38 
39  !-----------------------------------------------------------------------------
40  !
41  !++ Public parameters & variables
42  !
43  logical, public :: atmos_phy_rd_profile_use_climatology = .true.
44 
45  !-----------------------------------------------------------------------------
46  !
47  !++ Private procedure
48  !
49  private :: profile_setup_cira86
50  private :: profile_setup_mipas2001
51  private :: readfile_mipas2001
52  private :: profile_read_climatology
53  private :: profile_read_cira86
54  private :: profile_read_mipas2001
55  private :: profile_read_user
56  private :: profile_interp
57 
58  !-----------------------------------------------------------------------------
59  !
60  !++ Private parameters & variables
61  !
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_h2o = .true.
66  logical, private :: atmos_phy_rd_profile_use_co2 = .true.
67  logical, private :: atmos_phy_rd_profile_use_o3 = .true.
68  logical, private :: atmos_phy_rd_profile_use_n2o = .true.
69  logical, private :: atmos_phy_rd_profile_use_co = .true.
70  logical, private :: atmos_phy_rd_profile_use_ch4 = .true.
71  logical, private :: atmos_phy_rd_profile_use_o2 = .true.
72  logical, private :: atmos_phy_rd_profile_use_cfc = .true.
73  logical, private :: debug = .false.
74 
75  integer, private :: cira_ntime
76  integer, private :: cira_nplev
77  integer, private :: cira_nlat
78  real(RP), private, allocatable :: cira_nd (:) ! [day]
79  real(RP), private, allocatable :: cira_plog(:) ! log([hPa])
80  real(RP), private, allocatable :: cira_lat (:) ! [rad]
81  real(RP), private, allocatable :: cira_temp(:,:,:) ! [K]
82  real(RP), private, allocatable :: cira_z (:,:,:) ! [km]
83 
84  real(RP), private, allocatable :: interp_temp(:) ! [K]
85  real(RP), private, allocatable :: interp_z (:) ! [km]
86 
87  integer, private, parameter :: mipas_kmax = 121
88  integer, private, parameter :: mipas_ntime = 2
89  real(RP), private :: mipas_nd (0:mipas_ntime+1) ! [day]
90  real(RP), private :: mipas_lat (5) ! [rad]
91  real(RP), private :: mipas_z (mipas_kmax,4) ! [km]
92  real(RP), private :: mipas_pres(mipas_kmax,4) ! (not used) [hPa]
93  real(RP), private :: mipas_temp(mipas_kmax,4) ! (not used) [K]
94  real(RP), private :: mipas_gas (mipas_kmax,30,4) ! [ppmv]
95 
96  integer, private, parameter :: i_tropic = 1
97  integer, private, parameter :: i_midlat = 2
98  integer, private, parameter :: i_polarsum = 3
99  integer, private, parameter :: i_polarwin = 4
100 
101  integer, private, parameter :: i_n2 = 1
102  integer, private, parameter :: i_o2 = 2
103  integer, private, parameter :: i_co2 = 3
104  integer, private, parameter :: i_o3 = 4
105  integer, private, parameter :: i_h2o = 5
106  integer, private, parameter :: i_ch4 = 6
107  integer, private, parameter :: i_n2o = 7
108  integer, private, parameter :: i_hno3 = 8
109  integer, private, parameter :: i_co = 9
110  integer, private, parameter :: i_no2 = 10
111  integer, private, parameter :: i_n2o5 = 11
112  integer, private, parameter :: i_clo = 12
113  integer, private, parameter :: i_hocl = 13
114  integer, private, parameter :: i_clono2 = 14
115  integer, private, parameter :: i_no = 15
116  integer, private, parameter :: i_hno4 = 16
117  integer, private, parameter :: i_hcn = 17
118  integer, private, parameter :: i_nh3 = 18
119  integer, private, parameter :: i_f11 = 19
120  integer, private, parameter :: i_f12 = 20
121  integer, private, parameter :: i_f14 = 21
122  integer, private, parameter :: i_f22 = 22
123  integer, private, parameter :: i_ccl4 = 23
124  integer, private, parameter :: i_cof2 = 24
125  integer, private, parameter :: i_h2o2 = 25
126  integer, private, parameter :: i_c2h2 = 26
127  integer, private, parameter :: i_c2h6 = 27
128  integer, private, parameter :: i_ocs = 28
129  integer, private, parameter :: i_so2 = 29
130  integer, private, parameter :: i_sf6 = 30
131 
132  logical, private :: report_firsttime = .true.
133 
134  !-----------------------------------------------------------------------------
135 contains
136  !-----------------------------------------------------------------------------
138  subroutine atmos_phy_rd_profile_setup
139  use scale_process, only: &
140  prc_mpistop
141  implicit none
142 
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_use_climatology, &
149  atmos_phy_rd_profile_cira86_in_filename, &
150  atmos_phy_rd_profile_mipas2001_in_basename, &
151  atmos_phy_rd_profile_user_in_filename, &
152  atmos_phy_rd_profile_use_h2o, &
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[RADIATION PROFILE] / Categ[ATMOS PHYSICS] / Origin[SCALElib]'
167 
168  atmos_phy_rd_profile_cira86_in_filename = profile_cira86_fname
169  atmos_phy_rd_profile_mipas2001_in_basename = profile_mipas2001_dir
170  atmos_phy_rd_profile_user_in_filename = profile_user_fname
171 
172  !--- read namelist
173  rewind(io_fid_conf)
174  read(io_fid_conf,nml=param_atmos_phy_rd_profile,iostat=ierr)
175 
176  if( ierr < 0 ) then !--- missing
177  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
178  elseif( ierr > 0 ) then !--- fatal error
179  write(*,*) 'xxx Not appropriate names in namelist PARAM_ATMOS_PHY_RD_PROFILE. Check!'
180  call prc_mpistop
181  endif
182  if( io_nml ) write(io_fid_nml,nml=param_atmos_phy_rd_profile)
183 
184  profile_cira86_fname = atmos_phy_rd_profile_cira86_in_filename
185  profile_mipas2001_dir = atmos_phy_rd_profile_mipas2001_in_basename
186  profile_user_fname = atmos_phy_rd_profile_user_in_filename
187 
188  if( io_l ) write(io_fid_log,*)
189  if( io_l ) write(io_fid_log,*) '*** Climatological profile for radiation'
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,*)
230  if( io_l ) write(io_fid_log,*) '*** Read CIRA86 climatology, filename : ', trim(profile_cira86_fname)
231 
232  inquire( file=trim(profile_cira86_fname), exist=exist )
233  if ( .NOT. exist ) then !--- missing
234  write(*,*) '*** [PROFILE_setup_CIRA86] 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,*)
421  if( io_l ) write(io_fid_log,*) '*** Read MIPAS2001 climatology ***'
422 
423  mipas_date(:, 0) = (/ 2000, 12, 22, 12, 0, 0 /)
424  mipas_date(:, 1) = (/ 2001, 6, 21, 12, 0, 0 /)
425  mipas_date(:, 2) = (/ 2001, 12, 22, 12, 0, 0 /)
426  mipas_date(:, 3) = (/ 2002, 6, 21, 12, 0, 0 /)
427 
428  do t = 0, mipas_ntime+1
429  call calendar_date2daysec( nday, nsec, & ! [OUT]
430  mipas_date(:,t), subsec, & ! [IN]
431  offset_year ) ! [IN]
432 
433  mipas_nd(t) = real(nday,kind=RP) + nsec / 86400.0_rp
434  enddo
435 
436  mipas_lat(1) = -90.0_rp * const_d2r
437  mipas_lat(2) = -45.0_rp * const_d2r
438  mipas_lat(3) = 0.0_rp * const_d2r
439  mipas_lat(4) = 45.0_rp * const_d2r
440  mipas_lat(5) = 90.0_rp * const_d2r
441 
442  do rgn = i_tropic, i_polarwin
443  fname = trim(profile_mipas2001_dir)//'/'//mipas_fname(rgn)
444  if( io_l ) write(io_fid_log,*) '*** filename : ', trim(fname)
445 
446  fid = io_get_available_fid()
447  open( unit = fid, &
448  file = trim(fname), &
449  form = 'formatted', &
450  status = 'old', &
451  iostat = ierr )
452 
453  if ( ierr /= 0 ) then !--- missing
454  write(*,*) '*** [PROFILE_setup_MIPAS2001] File not found. check!'
455  call prc_mpistop
456  endif
457 
458  do l = 1, 24
459  read(fid,*) dummy
460  enddo
461 
462  call readfile_mipas2001( fid, mipas_z(:,rgn) )
463  call readfile_mipas2001( fid, mipas_pres(:,rgn) )
464  call readfile_mipas2001( fid, mipas_temp(:,rgn) )
465 
466  call readfile_mipas2001( fid, mipas_gas(:,i_n2 ,rgn) )
467  call readfile_mipas2001( fid, mipas_gas(:,i_o2 ,rgn) )
468  call readfile_mipas2001( fid, mipas_gas(:,i_co2 ,rgn) )
469  call readfile_mipas2001( fid, mipas_gas(:,i_o3 ,rgn) )
470  call readfile_mipas2001( fid, mipas_gas(:,i_h2o ,rgn) )
471  call readfile_mipas2001( fid, mipas_gas(:,i_ch4 ,rgn) )
472  call readfile_mipas2001( fid, mipas_gas(:,i_n2o ,rgn) )
473  call readfile_mipas2001( fid, mipas_gas(:,i_hno3 ,rgn) )
474  call readfile_mipas2001( fid, mipas_gas(:,i_co ,rgn) )
475  call readfile_mipas2001( fid, mipas_gas(:,i_no2 ,rgn) )
476  call readfile_mipas2001( fid, mipas_gas(:,i_n2o5 ,rgn) )
477  call readfile_mipas2001( fid, mipas_gas(:,i_clo ,rgn) )
478  call readfile_mipas2001( fid, mipas_gas(:,i_hocl ,rgn) )
479  call readfile_mipas2001( fid, mipas_gas(:,i_clono2,rgn) )
480  call readfile_mipas2001( fid, mipas_gas(:,i_no ,rgn) )
481  call readfile_mipas2001( fid, mipas_gas(:,i_hno4 ,rgn) )
482  call readfile_mipas2001( fid, mipas_gas(:,i_hcn ,rgn) )
483  call readfile_mipas2001( fid, mipas_gas(:,i_nh3 ,rgn) )
484  call readfile_mipas2001( fid, mipas_gas(:,i_f11 ,rgn) )
485  call readfile_mipas2001( fid, mipas_gas(:,i_f12 ,rgn) )
486  call readfile_mipas2001( fid, mipas_gas(:,i_f14 ,rgn) )
487  call readfile_mipas2001( fid, mipas_gas(:,i_f22 ,rgn) )
488  call readfile_mipas2001( fid, mipas_gas(:,i_ccl4 ,rgn) )
489  call readfile_mipas2001( fid, mipas_gas(:,i_cof2 ,rgn) )
490  call readfile_mipas2001( fid, mipas_gas(:,i_h2o2 ,rgn) )
491  call readfile_mipas2001( fid, mipas_gas(:,i_c2h2 ,rgn) )
492  call readfile_mipas2001( fid, mipas_gas(:,i_c2h6 ,rgn) )
493  call readfile_mipas2001( fid, mipas_gas(:,i_ocs ,rgn) )
494  call readfile_mipas2001( fid, mipas_gas(:,i_so2 ,rgn) )
495  call readfile_mipas2001( fid, mipas_gas(:,i_sf6 ,rgn) )
496  close(fid)
497  enddo
498 
499  return
500  end subroutine profile_setup_mipas2001
501 
502  !-----------------------------------------------------------------------------
504  subroutine readfile_mipas2001( &
505  fid, &
506  var )
507  implicit none
508 
509  integer, intent(in) :: fid
510  real(RP), intent(out) :: var(121)
511 
512  character(len=H_LONG) :: dummy
513  real(RP) :: tmp5(5), tmp1
514 
515  integer :: nstr, l
516  !---------------------------------------------------------------------------
517 
518  read(fid,*) dummy
519  !if( IO_L ) write(IO_FID_LOG,*) dummy
520 
521  nstr = mipas_kmax
522  do l = 1, 24
523  read(fid,*) tmp5(:)
524  !if( IO_L ) write(IO_FID_LOG,*) l, tmp5(:)
525  var(nstr ) = tmp5(1)
526  var(nstr-1) = tmp5(2)
527  var(nstr-2) = tmp5(3)
528  var(nstr-3) = tmp5(4)
529  var(nstr-4) = tmp5(5)
530 
531  nstr = nstr - 5
532  enddo
533 
534  read(fid,*) tmp1
535  var(nstr) = tmp1
536 
537  return
538  end subroutine readfile_mipas2001
539 
540  !-----------------------------------------------------------------------------
542  subroutine atmos_phy_rd_profile_read( &
543  kmax, &
544  ngas, &
545  ncfc, &
546  naero, &
547  real_lat, &
548  now_date, &
549  zh, &
550  z, &
551  rhodz, &
552  pres, &
553  presh, &
554  temp, &
555  temph, &
556  gas, &
557  cfc, &
558  aerosol_conc, &
559  aerosol_radi, &
560  cldfrac )
561  use scale_const, only: &
562  grav => const_grav
563  use scale_atmos_solarins, only: &
564  solarins_fixedlatlon => atmos_solarins_fixedlatlon, &
565  solarins_fixeddate => atmos_solarins_fixeddate, &
566  solarins_lat => atmos_solarins_lat, &
567  solarins_date => atmos_solarins_date
568  implicit none
569 
570  integer, intent(in) :: kmax
571  integer, intent(in) :: ngas
572  integer, intent(in) :: ncfc
573  integer, intent(in) :: naero
574  real(RP), intent(in) :: real_lat
575  integer, intent(in) :: now_date(6)
576  real(RP), intent(in) :: zh(kmax+1)
577  real(RP), intent(in) :: z (kmax)
578  real(RP), intent(out) :: rhodz (kmax)
579  real(RP), intent(out) :: pres (kmax)
580  real(RP), intent(out) :: presh (kmax+1)
581  real(RP), intent(out) :: temp (kmax)
582  real(RP), intent(out) :: temph (kmax+1)
583  real(RP), intent(out) :: gas (kmax,ngas)
584  real(RP), intent(out) :: cfc (kmax,ncfc)
585  real(RP), intent(out) :: aerosol_conc(kmax,naero)
586  real(RP), intent(out) :: aerosol_radi(kmax,naero)
587  real(RP), intent(out) :: cldfrac (kmax)
588 
589  real(RP) :: lat
590  integer :: date(6)
591 
592  integer :: k
593  !---------------------------------------------------------------------------
594 
595  if ( atmos_phy_rd_profile_use_climatology ) then
596 
597  if ( solarins_fixedlatlon ) then
598  lat = solarins_lat
599  else
600  lat = real_lat
601  endif
602 
603  date = now_date
604  if ( solarins_fixeddate ) then
605  if( solarins_date(1) >= 0 ) date(1) = solarins_date(1)
606  if( solarins_date(2) >= 1 ) date(2) = solarins_date(2)
607  if( solarins_date(3) >= 1 ) date(3) = solarins_date(3)
608  if( solarins_date(4) >= 0 ) date(4) = solarins_date(4)
609  if( solarins_date(5) >= 0 ) date(5) = solarins_date(5)
610  if( solarins_date(6) >= 0 ) date(6) = solarins_date(6)
611  endif
612 
613  call profile_read_climatology( kmax, & ! [IN]
614  ngas, & ! [IN]
615  ncfc, & ! [IN]
616  naero, & ! [IN]
617  lat, & ! [IN], tentative treatment
618  date, & ! [IN]
619  zh, & ! [IN]
620  z, & ! [IN]
621  pres, & ! [OUT]
622  presh, & ! [OUT]
623  temp, & ! [OUT]
624  temph, & ! [OUT]
625  gas, & ! [OUT]
626  cfc ) ! [OUT]
627 
628  else
629 
630  call profile_read_user( kmax, & ! [IN]
631  ngas, & ! [IN]
632  ncfc, & ! [IN]
633  naero, & ! [IN]
634  zh, & ! [IN]
635  z, & ! [IN]
636  pres, & ! [OUT]
637  presh, & ! [OUT]
638  temp, & ! [OUT]
639  temph, & ! [OUT]
640  gas, & ! [OUT]
641  cfc ) ! [OUT]
642 
643  endif
644 
645 
646  do k = 1, kmax
647  rhodz(k) = ( presh(k+1) - presh(k) ) * 100.0_rp / grav
648  enddo
649 
650  ! no cloud/aerosol
651  aerosol_conc(:,:) = 0.0_rp
652  aerosol_radi(:,:) = 0.0_rp
653  cldfrac(:) = 0.0_rp
654 
655  if ( .NOT. atmos_phy_rd_profile_use_h2o ) gas(:,1) = 0.0_rp
656  if ( .NOT. atmos_phy_rd_profile_use_co2 ) gas(:,2) = 0.0_rp
657  if ( .NOT. atmos_phy_rd_profile_use_o3 ) gas(:,3) = 0.0_rp
658  if ( .NOT. atmos_phy_rd_profile_use_n2o ) gas(:,4) = 0.0_rp
659  if ( .NOT. atmos_phy_rd_profile_use_co ) gas(:,5) = 0.0_rp
660  if ( .NOT. atmos_phy_rd_profile_use_ch4 ) gas(:,6) = 0.0_rp
661  if ( .NOT. atmos_phy_rd_profile_use_o2 ) gas(:,7) = 0.0_rp
662  if ( .NOT. atmos_phy_rd_profile_use_cfc ) cfc(:,:) = 0.0_rp
663 
664  !----- report data -----
665  if ( debug .AND. report_firsttime ) then
666  report_firsttime = .false.
667 
668  if( io_l ) write(io_fid_log,*)
669  if( io_l ) write(io_fid_log,'(1x,A)') &
670  '|=============== Vertical Coordinate =================|'
671  if( io_l ) write(io_fid_log,'(1x,A)') &
672  '| -GRID CENTER- -GRID INTERFACE- |'
673  if( io_l ) write(io_fid_log,'(1x,A)') &
674  '| k z pres temp zh pres temp k |'
675  if( io_l ) write(io_fid_log,'(1x,A)') &
676  '| [km] [hPa] [K] [km] [hPa] [K] |'
677  k = 1
678  if( io_l ) write(io_fid_log,'(1x,A,F8.3,F10.4,F8.2,I5,A)') &
679  '| ',zh(k),presh(k),temph(k),k, ' | TOA'
680  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,F10.4,F8.2,A)') &
681  '|',k,z(k),pres(k),temp(k), ' | '
682  do k = 2, kmax-1
683  if( io_l ) write(io_fid_log,'(1x,A,F8.3,F10.4,F8.2,I5,A)') &
684  '| ',zh(k),presh(k),temph(k),k, ' | '
685  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,F10.4,F8.2,A)') &
686  '|',k,z(k),pres(k),temp(k), ' | '
687  enddo
688  k = kmax
689  if( io_l ) write(io_fid_log,'(1x,A,F8.3,F10.4,F8.2,I5,A)') &
690  '| ',zh(k),presh(k),temph(k),k, ' | '
691  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,F10.4,F8.2,A)') &
692  '|',k,z(k),pres(k),temp(k), ' | '
693  k = kmax+1
694  if( io_l ) write(io_fid_log,'(1x,A,F8.3,F10.4,F8.2,I5,A)') &
695  '| ',zh(k),presh(k),temph(k),k, ' | Ground'
696  if( io_l ) write(io_fid_log,'(1x,A)') &
697  '|================================================================|'
698 
699  if( io_l ) write(io_fid_log,*)
700  if( io_l ) write(io_fid_log,'(1x,A)') &
701  '|=====================================================================================|'
702  if( io_l ) write(io_fid_log,'(1x,A)') &
703  '| -Gas concetrations [ppmv]- |'
704  if( io_l ) write(io_fid_log,'(1x,A)') &
705  '| k z H2O CO2 O3 N2O CO CH4 O2 |'
706  do k = 1, kmax
707  if( io_l ) write(io_fid_log,'(1x,A,I5,1F9.3,7ES10.3,A)') '|',k,z(k),gas(k,:),' | '
708  enddo
709  if( io_l ) write(io_fid_log,'(1x,A)') &
710  '|=====================================================================================|'
711 
712  endif
713 
714  return
715  end subroutine atmos_phy_rd_profile_read
716 
717  !-----------------------------------------------------------------------------
719  subroutine profile_read_climatology( &
720  kmax, &
721  ngas, &
722  ncfc, &
723  naero, &
724  lat, &
725  now_date, &
726  zh, &
727  z, &
728  pres, &
729  presh, &
730  temp, &
731  temph, &
732  gas, &
733  cfc )
734  implicit none
735 
736  integer, intent(in) :: kmax
737  integer, intent(in) :: ngas
738  integer, intent(in) :: ncfc
739  integer, intent(in) :: naero
740  real(RP), intent(in) :: lat
741  integer, intent(in) :: now_date(6)
742  real(RP), intent(in) :: zh (kmax+1)
743  real(RP), intent(in) :: z (kmax)
744  real(RP), intent(out) :: pres (kmax)
745  real(RP), intent(out) :: presh(kmax+1)
746  real(RP), intent(out) :: temp (kmax)
747  real(RP), intent(out) :: temph(kmax+1)
748  real(RP), intent(out) :: gas (kmax,ngas)
749  real(RP), intent(out) :: cfc (kmax,ncfc)
750  !---------------------------------------------------------------------------
751 
752  if( io_l ) write(io_fid_log,*) '*** Update climatological profile for radiation'
753 
754  call profile_read_cira86( kmax, & ! [IN]
755  lat, & ! [IN]
756  now_date(:), & ! [IN]
757  zh(:), & ! [IN]
758  z(:), & ! [IN]
759  presh(:), & ! [OUT]
760  temph(:), & ! [OUT]
761  pres(:), & ! [OUT]
762  temp(:) ) ! [OUT]
763 
764  call profile_read_mipas2001( kmax, & ! [IN]
765  ngas, & ! [IN]
766  ncfc, & ! [IN]
767  lat, & ! [IN]
768  now_date(:), & ! [IN]
769  z(:), & ! [IN]
770  gas(:,:), & ! [OUT]
771  cfc(:,:) ) ! [OUT]
772 
773  ! no vapor
774  gas(:,1) = 0.0_rp
775 
776  return
777  end subroutine profile_read_climatology
778 
779  !-----------------------------------------------------------------------------
781  subroutine profile_read_cira86( &
782  kmax, &
783  lat, &
784  now_date, &
785  zh, &
786  z, &
787  presh, &
788  temph, &
789  pres, &
790  temp )
791  use scale_calendar, only: &
792  calendar_date2daysec
793  implicit none
794 
795  integer, intent(in) :: kmax
796  real(RP), intent(in) :: lat
797  integer, intent(in) :: now_date(6)
798  real(RP), intent(in) :: zh (kmax+1)
799  real(RP), intent(in) :: z (kmax)
800  real(RP), intent(out) :: presh(kmax+1)
801  real(RP), intent(out) :: temph(kmax+1)
802  real(RP), intent(out) :: pres (kmax)
803  real(RP), intent(out) :: temp (kmax)
804 
805  real(RP) :: plogh(kmax+1)
806  real(RP) :: plog (kmax)
807 
808  integer :: now_date_mod(6), nday
809  real(RP) :: nd
810  real(DP) :: nsec
811  real(DP) :: subsec = 0.0_dp
812  integer :: offset_year = 0
813 
814  integer :: nplev_mod
815  integer :: indexlat, indexd
816  real(RP) :: factlat, factd
817 
818  integer :: n, t
819  !---------------------------------------------------------------------------
820 
821  ! latitude interpolation
822  if ( lat < cira_lat(1) ) then ! extrapolation
823  indexlat = 1
824  factlat = 0.0_rp
825  elseif( lat >= cira_lat(cira_nlat) ) then ! extrapolation
826  indexlat = cira_nlat - 1
827  factlat = 1.0_rp
828  else
829  do n = 1, cira_nlat-1
830  if ( lat >= cira_lat(n ) &
831  .AND. lat < cira_lat(n+1) ) then ! interpolation
832  indexlat = n
833  factlat = ( lat-cira_lat(n) ) / ( cira_lat(n+1)-cira_lat(n) )
834  endif
835  enddo
836  endif
837 
838  ! time interpolation
839  now_date_mod(2:6) = now_date(2:6)
840  now_date_mod(1) = 1986
841 
842  call calendar_date2daysec( nday, nsec, & ! [OUT]
843  now_date_mod(:), subsec, & ! [IN]
844  offset_year ) ! [IN]
845 
846  nd = real(nday,kind=RP) + nsec / 86400.0_rp
847 
848  do t = 0, cira_ntime
849  if ( nd >= cira_nd(t ) &
850  .AND. nd < cira_nd(t+1) ) then ! interpolation
851  indexd = t
852  factd = ( nd-cira_nd(t) ) / ( cira_nd(t+1)-cira_nd(t) )
853  endif
854  enddo
855 
856  interp_z(:) = cira_z(:,indexlat ,indexd ) * ( 1.0_rp-factlat ) * ( 1.0_rp-factd ) &
857  + cira_z(:,indexlat+1,indexd ) * ( factlat ) * ( 1.0_rp-factd ) &
858  + cira_z(:,indexlat ,indexd+1) * ( 1.0_rp-factlat ) * ( factd ) &
859  + cira_z(:,indexlat+1,indexd+1) * ( factlat ) * ( factd )
860 
861  interp_temp(:) = cira_temp(:,indexlat ,indexd ) * ( 1.0_rp-factlat ) * ( 1.0_rp-factd ) &
862  + cira_temp(:,indexlat+1,indexd ) * ( factlat ) * ( 1.0_rp-factd ) &
863  + cira_temp(:,indexlat ,indexd+1) * ( 1.0_rp-factlat ) * ( factd ) &
864  + cira_temp(:,indexlat+1,indexd+1) * ( factlat ) * ( factd )
865 
866  ! avoid missing value
867  nplev_mod = cira_nplev
868  do n = cira_nplev, 1, -1
869  if ( interp_temp(n) == interp_temp(n-1) ) then
870  nplev_mod = nplev_mod-1
871  else
872  exit
873  endif
874  enddo
875 
876  !do n = 1, nplev_mod
877  ! print *, n, interp_z(n), interp_temp(n), indexLAT, indexD, factLAT, factD
878  !enddo
879 
880  ! pressure interpolation
881  call profile_interp( nplev_mod, & ! [IN]
882  interp_z(1:nplev_mod), & ! [IN]
883  cira_plog(1:nplev_mod), & ! [IN]
884  kmax+1, & ! [IN]
885  zh(:), & ! [IN]
886  plogh(:) ) ! [OUT]
887 
888  presh(:) = exp( plogh(:) )
889 
890  call profile_interp( kmax+1, zh(:), plogh(:), kmax, z(:), plog(:) )
891  pres(:) = exp( plog(:) )
892 
893  call profile_interp( nplev_mod, & ! [IN]
894  interp_z(1:nplev_mod), & ! [IN]
895  interp_temp(1:nplev_mod), & ! [IN]
896  kmax, & ! [IN]
897  z(:), & ! [IN]
898  temp(:) ) ! [OUT]
899 
900  call profile_interp( nplev_mod, & ! [IN]
901  interp_z(1:nplev_mod), & ! [IN]
902  interp_temp(1:nplev_mod), & ! [IN]
903  kmax+1, & ! [IN]
904  zh(:), & ! [IN]
905  temph(:) ) ! [OUT]
906 
907  return
908  end subroutine profile_read_cira86
909 
910  !-----------------------------------------------------------------------------
912  subroutine profile_read_mipas2001( &
913  kmax, &
914  ngas, &
915  ncfc, &
916  lat, &
917  now_date, &
918  z, &
919  gas, &
920  cfc )
921  use scale_calendar, only: &
922  calendar_date2daysec
923  implicit none
924 
925  integer, intent(in) :: kmax
926  integer, intent(in) :: ngas
927  integer, intent(in) :: ncfc
928  real(RP), intent(in) :: lat
929  integer, intent(in) :: now_date(6)
930  real(RP), intent(in) :: z (kmax)
931  real(RP), intent(inout) :: gas(kmax,ngas)
932  real(RP), intent(inout) :: cfc(kmax,ncfc)
933 
934  real(RP) :: interp_gas(mipas_kmax,30)
935  real(RP) :: interp_z (mipas_kmax)
936 
937  integer :: now_date_mod(6), nday
938  real(RP) :: nd
939  real(DP) :: nsec
940  real(DP) :: subsec = 0.0_dp
941  integer :: offset_year = 0
942 
943  integer :: indexd1, indexd2
944  real(RP) :: factlat, factd
945  !---------------------------------------------------------------------------
946 
947  ! time interpolation
948  now_date_mod(2:6) = now_date(2:6)
949  now_date_mod(1) = 2001
950 
951  call calendar_date2daysec( nday, nsec, & ! [OUT]
952  now_date_mod(:), subsec, & ! [IN]
953  offset_year ) ! [IN]
954 
955  nd = real(nday,kind=RP) + nsec / 86400.0_rp
956 
957  if ( nd >= mipas_nd(0) .AND. nd < mipas_nd(1) ) then ! winter-summer
958 
959  indexd1 = i_polarwin
960  indexd2 = i_polarsum
961  factd = ( nd-mipas_nd(0) ) / ( mipas_nd(1)-mipas_nd(0) )
962 
963  elseif( nd >= mipas_nd(1) .AND. nd < mipas_nd(2) ) then ! summer-winter
964 
965  indexd1 = i_polarsum
966  indexd2 = i_polarwin
967  factd = ( nd-mipas_nd(1) ) / ( mipas_nd(2)-mipas_nd(1) )
968 
969  elseif( nd >= mipas_nd(2) .AND. nd < mipas_nd(3) ) then ! winter-summer
970 
971  indexd1 = i_polarwin
972  indexd2 = i_polarsum
973  factd = ( nd-mipas_nd(2) ) / ( mipas_nd(3)-mipas_nd(2) )
974 
975  endif
976 
977  ! latitude interpolation
978  if ( lat < mipas_lat(1) ) then ! south pole
979 
980  interp_gas(:,:) = mipas_gas(:,:,indexd1 ) * ( 1.0_rp-factd ) &
981  + mipas_gas(:,:,indexd2 ) * ( factd )
982 
983  interp_z(:) = mipas_z(:,indexd1 ) * ( 1.0_rp-factd ) &
984  + mipas_z(:,indexd2 ) * ( factd )
985 
986  elseif( lat >= mipas_lat(1) .AND. lat < mipas_lat(2) ) then ! south pole-SH mid
987 
988  factlat = ( lat-mipas_lat(1) ) / ( mipas_lat(2)-mipas_lat(1) )
989 
990  interp_gas(:,:) = mipas_gas(:,:,indexd1) * ( 1.0_rp-factd ) * ( 1.0_rp-factlat ) &
991  + mipas_gas(:,:,indexd2) * ( factd ) * ( 1.0_rp-factlat ) &
992  + mipas_gas(:,:,i_midlat) * ( factlat )
993 
994  interp_z(:) = mipas_z(:,indexd1) * ( 1.0_rp-factd ) * ( 1.0_rp-factlat ) &
995  + mipas_z(:,indexd2) * ( factd ) * ( 1.0_rp-factlat ) &
996  + mipas_z(:,i_midlat) * ( factlat )
997 
998  elseif( lat >= mipas_lat(2) .AND. lat < mipas_lat(3) ) then ! SH mid-EQ
999 
1000  factlat = ( lat-mipas_lat(2) ) / ( mipas_lat(3)-mipas_lat(2) )
1001 
1002  interp_gas(:,:) = mipas_gas(:,:,i_midlat) * ( 1.0_rp-factlat ) &
1003  + mipas_gas(:,:,i_tropic) * ( factlat )
1004 
1005  interp_z(:) = mipas_z(:,i_midlat) * ( 1.0_rp-factlat ) &
1006  + mipas_z(:,i_tropic) * ( factlat )
1007 
1008  elseif( lat >= mipas_lat(3) .AND. lat < mipas_lat(4) ) then ! EQ-NH mid
1009 
1010  factlat = ( lat-mipas_lat(3) ) / ( mipas_lat(4)-mipas_lat(3) )
1011 
1012  interp_gas(:,:) = mipas_gas(:,:,i_tropic) * ( 1.0_rp-factlat ) &
1013  + mipas_gas(:,:,i_midlat) * ( factlat )
1014 
1015  interp_z(:) = mipas_z(:,i_tropic) * ( 1.0_rp-factlat ) &
1016  + mipas_z(:,i_midlat) * ( factlat )
1017 
1018  elseif( lat >= mipas_lat(4) .AND. lat < mipas_lat(5) ) then ! NH mid-north pole
1019 
1020  factlat = ( lat-mipas_lat(4) ) / ( mipas_lat(5)-mipas_lat(4) )
1021 
1022  interp_gas(:,:) = mipas_gas(:,:,i_midlat) * ( 1.0_rp-factlat ) &
1023  + mipas_gas(:,:,indexd2) * ( 1.0_rp-factd ) * ( factlat ) &
1024  + mipas_gas(:,:,indexd1) * ( factd ) * ( factlat )
1025 
1026  interp_z(:) = mipas_z(:,i_midlat) * ( 1.0_rp-factlat ) &
1027  + mipas_z(:,indexd2) * ( 1.0_rp-factd ) * ( factlat ) &
1028  + mipas_z(:,indexd1) * ( factd ) * ( factlat )
1029 
1030  elseif( lat >= mipas_lat(5) ) then ! north pole
1031 
1032  interp_gas(:,:) = mipas_gas(:,:,indexd2) * ( 1.0_rp-factd ) &
1033  + mipas_gas(:,:,indexd1) * ( factd )
1034 
1035  interp_z(:) = mipas_z(:,indexd2) * ( 1.0_rp-factd ) &
1036  + mipas_z(:,indexd1) * ( factd )
1037 
1038  endif
1039 
1040  gas(:,:) = 0.0_rp
1041  cfc(:,:) = 0.0_rp
1042 
1043  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_h2o ), kmax, z(:), gas(:,1) )
1044  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_co2 ), kmax, z(:), gas(:,2) )
1045  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_o3 ), kmax, z(:), gas(:,3) )
1046  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_n2o ), kmax, z(:), gas(:,4) )
1047  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_co ), kmax, z(:), gas(:,5) )
1048  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_ch4 ), kmax, z(:), gas(:,6) )
1049  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_o2 ), kmax, z(:), gas(:,7) )
1050 
1051  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f11 ), kmax, z(:), cfc(:, 1) )
1052  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f12 ), kmax, z(:), cfc(:, 2) )
1053  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f14 ), kmax, z(:), cfc(:, 4) )
1054  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_f22 ), kmax, z(:), cfc(:, 9) )
1055  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_sf6 ), kmax, z(:), cfc(:,22) )
1056  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_clono2), kmax, z(:), cfc(:,23) )
1057  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_ccl4 ), kmax, z(:), cfc(:,24) )
1058  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_n2o5 ), kmax, z(:), cfc(:,25) )
1059  call profile_interp( mipas_kmax, interp_z(:), interp_gas(:,i_hno4 ), kmax, z(:), cfc(:,27) )
1060 
1061  return
1062  end subroutine profile_read_mipas2001
1063 
1064  !-----------------------------------------------------------------------------
1066  subroutine profile_interp( imax1, x1, y1, imax2, x2, y2 )
1067  implicit none
1068 
1069  integer, intent(in) :: imax1
1070  real(RP), intent(in) :: x1(imax1)
1071  real(RP), intent(in) :: y1(imax1)
1072  integer, intent(in) :: imax2
1073  real(RP), intent(in) :: x2(imax2)
1074  real(RP), intent(out) :: y2(imax2)
1075 
1076  real(RP) :: fact
1077 
1078  integer :: i1, i2
1079  !---------------------------------------------------------------------------
1080 
1081  !$omp parallel do default(none) &
1082  !$omp shared(imax1,x2,x1,y2,y1,imax2) &
1083  !$omp private(i1,fact,i2) OMP_SCHEDULE_
1084  do i2 = 1, imax2
1085 
1086  if ( x2(i2) > x1(1) ) then ! extrapolation
1087 
1088  fact = ( x1(1) - x2(i2) ) / ( x1(2) - x1(1) )
1089 
1090  y2(i2) = y1(1) * ( 1.0_rp-fact ) &
1091  + y1(2) * ( fact )
1092 
1093  elseif( x2(i2) <= x1(imax1) ) then ! extrapolation
1094 
1095  fact = ( x1(imax1) - x2(i2) ) / ( x1(imax1) - x1(imax1-1) )
1096 
1097  y2(i2) = y1(imax1-1) * ( fact ) &
1098  + y1(imax1 ) * ( 1.0_rp-fact )
1099 
1100  else
1101  do i1 = 1, imax1-1
1102  if ( x2(i2) <= x1(i1 ) &
1103  .AND. x2(i2) > x1(i1+1) ) then ! interpolation
1104 
1105  fact = ( x2(i2) - x1(i1) ) / ( x1(i1+1) - x1(i1) )
1106 
1107  y2(i2) = y1(i1 ) * ( 1.0_rp-fact ) &
1108  + y1(i1+1) * ( fact )
1109 
1110  exit
1111  endif
1112  enddo
1113  endif
1114 
1115  enddo
1116 
1117  return
1118  end subroutine profile_interp
1119 
1120  !-----------------------------------------------------------------------------
1122  subroutine atmos_phy_rd_profile_setup_zgrid( &
1123  toa, &
1124  kmax, &
1125  kadd, &
1126  zh, &
1127  z )
1128  use scale_grid, only: &
1129  cz => grid_cz, &
1130  fz => grid_fz
1131  implicit none
1132 
1133  real(RP), intent(in) :: toa
1134  integer, intent(in) :: kmax
1135  integer, intent(in) :: kadd
1136  real(RP), intent(inout) :: zh(kmax+1)
1137  real(RP), intent(inout) :: z (kmax)
1138 
1139  real(RP) :: dz
1140  integer :: k, rd_k
1141  !---------------------------------------------------------------------------
1142 
1143  if ( kadd > 0 ) then
1144  !--- additional layer over the computational domain
1145  dz = ( toa - fz(ke)*1.e-3_rp ) / real( kadd, kind=rp )
1146 
1147  zh(1) = toa
1148  do k = 2, kadd
1149  zh(k) = zh(k-1) - dz
1150  enddo
1151  zh(kadd+1) = fz(ke)*1.e-3_rp
1152 
1153  ! linear interpolation for center of layers
1154  do k = 1, kadd
1155  z(k) = 0.5_rp * ( zh(k+1) + zh(k) )
1156  enddo
1157  endif
1158 
1159  !--- in computational domain [NOTE] layer index is reversed
1160  do k = ks-1, ke
1161  rd_k = kmax - ( k - ks )
1162  zh(rd_k) = fz(k)*1.e-3_rp
1163  enddo
1164  do k = ks, ke
1165  rd_k = kmax - ( k - ks )
1166  z(rd_k) = cz(k)*1.e-3_rp
1167  enddo
1168 
1169  !----- report data -----
1170  if ( debug ) then
1171 
1172  if( io_l ) write(io_fid_log,*)
1173  if( io_l ) write(io_fid_log,'(1x,A)') &
1174  '|=============== Vertical Coordinate ===============|'
1175  if( io_l ) write(io_fid_log,'(1x,A)') &
1176  '| -GRID CENTER- -GRID INTERFACE- |'
1177  if( io_l ) write(io_fid_log,'(1x,A)') &
1178  '| RD_k z k GRID_CZ GRID_FZ k zh RD_k |'
1179  if ( kadd > 0 ) then
1180  rd_k = 1
1181  if( io_l ) write(io_fid_log,'(1x,A,F8.3,I5,A)') &
1182  '| ',zh(rd_k),rd_k, ' | TOA'
1183  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,A)') &
1184  '|',rd_k,z(rd_k), ' | '
1185  do rd_k = 2, kadd-1
1186  if( io_l ) write(io_fid_log,'(1x,A,F8.3,I5,A)') &
1187  '| ',zh(rd_k),rd_k, ' | '
1188  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,A)') &
1189  '|',rd_k,z(rd_k), ' | '
1190  enddo
1191  rd_k = kadd
1192  if( io_l ) write(io_fid_log,'(1x,A,F8.3,I5,A)') &
1193  '| ',zh(rd_k),rd_k, ' | '
1194  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,A)') &
1195  '|',rd_k,z(rd_k), ' | KADD'
1196  rd_k = kadd+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,' | '
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, ' | KADD+1=KE'
1202  else
1203  rd_k = 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,' | TOA=KE'
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  endif
1210  do rd_k = kadd+2, kmax-1
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)*1.e-3_rp,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, ' | '
1216  enddo
1217  rd_k = kmax
1218  k = kmax - rd_k + ks
1219  if( io_l ) write(io_fid_log,'(1x,A,F8.3,I5,F8.3,I5,A)') &
1220  '| ',fz(k),k,zh(rd_k),rd_k, ' | '
1221  if( io_l ) write(io_fid_log,'(1x,A,I5,F8.3,I5,F8.3,A)') &
1222  '|',rd_k,z(rd_k),k,cz(k)*1.e-3_rp, ' | RD_KMAX=KS'
1223  rd_k = kmax+1
1224  k = kmax - rd_k + ks
1225  if( io_l ) write(io_fid_log,'(1x,A,F8.3,I5,F8.3,I5,A)') &
1226  '| ',fz(k)*1.e-3_rp,k,zh(rd_k),rd_k,' | Ground'
1227  if( io_l ) write(io_fid_log,'(1x,A)') &
1228  '|=====================================================|'
1229 
1230  endif
1231 
1232  return
1233  end subroutine atmos_phy_rd_profile_setup_zgrid
1234 
1235  !-----------------------------------------------------------------------------
1237  subroutine profile_read_user( &
1238  kmax, &
1239  ngas, &
1240  ncfc, &
1241  naero, &
1242  zh, &
1243  z, &
1244  pres, &
1245  presh, &
1246  temp, &
1247  temph, &
1248  gas, &
1249  cfc )
1250  use scale_process, only: &
1251  prc_mpistop
1252  use scale_const, only: &
1253  mdry => const_mdry, &
1254  mvap => const_mvap, &
1255  ppm => const_ppm
1256  implicit none
1257 
1258  integer, intent(in) :: kmax
1259  integer, intent(in) :: ngas
1260  integer, intent(in) :: ncfc
1261  integer, intent(in) :: naero
1262  real(RP), intent(in) :: zh (kmax+1)
1263  real(RP), intent(in) :: z (kmax)
1264  real(RP), intent(out) :: pres (kmax)
1265  real(RP), intent(out) :: presh(kmax+1)
1266  real(RP), intent(out) :: temp (kmax)
1267  real(RP), intent(out) :: temph(kmax+1)
1268  real(RP), intent(out) :: gas (kmax,ngas)
1269  real(RP), intent(out) :: cfc (kmax,ncfc)
1270 
1271  integer, parameter :: user_klim = 500
1272  integer :: user_kmax
1273  real(RP) :: user_z (user_klim)
1274  real(RP) :: user_pres(user_klim)
1275  real(RP) :: user_temp(user_klim)
1276  real(RP) :: user_qv (user_klim)
1277  real(RP) :: user_o3 (user_klim)
1278 
1279  real(RP), allocatable :: work_z(:)
1280  real(RP), allocatable :: work (:)
1281 
1282  real(RP) :: plog (kmax)
1283  real(RP) :: plogh(kmax+1)
1284 
1285  character(len=H_LONG) :: dummy
1286 
1287  integer :: fid, ierr
1288  integer :: k
1289  !---------------------------------------------------------------------------
1290 
1291  if( io_l ) write(io_fid_log,*) '*** [RD_PROFILE] user-defined profile'
1292 
1293  gas(:,:) = 0.0_rp
1294  cfc(:,:) = 0.0_rp
1295 
1296  if( io_l ) write(io_fid_log,*) '*** FILENAME:', trim(profile_user_fname)
1297 
1298  fid = io_get_available_fid()
1299  open( unit = fid, &
1300  file = trim(profile_user_fname), &
1301  form = 'formatted', &
1302  status = 'old', &
1303  iostat = ierr )
1304 
1305  if ( ierr /= 0 ) then !--- missing
1306  write(*,*) '*** [PROFILE_read_user] File not found. check!'
1307  call prc_mpistop
1308  endif
1309 
1310  read(fid,*) dummy
1311 
1312  do k = 1, user_klim
1313  read(fid,*,iostat=ierr) user_z(k), user_pres(k), user_temp(k), user_qv(k), user_o3(k)
1314  if ( ierr /= 0 ) exit
1315  enddo
1316  user_kmax = k - 1
1317  close(fid)
1318 
1319  allocate( work_z(user_kmax) )
1320  allocate( work(user_kmax) )
1321 
1322  do k = 1, user_kmax
1323  work_z(k) = user_z(k) / 1000.0_rp ! [m->km]
1324  work(k) = log( user_pres(k)/100.0_rp ) ! log[Pa->hPA]
1325  enddo
1326 
1327  ! pressure interpolation
1328  call profile_interp( user_kmax, & ! [IN]
1329  work_z(:), & ! [IN]
1330  work(:), & ! [IN]
1331  kmax+1, & ! [IN]
1332  zh(:), & ! [IN]
1333  plogh(:) ) ! [OUT]
1334 
1335  presh(:) = exp( plogh(:) )
1336 
1337  call profile_interp( kmax+1, zh(:), plogh(:), kmax, z(:), plog(:) )
1338  pres(:) = exp( plog(:) )
1339 
1340  do k = 1, user_kmax
1341  work(k) = user_temp(k)
1342  enddo
1343 
1344  call profile_interp( user_kmax, & ! [IN]
1345  work_z(:), & ! [IN]
1346  work(:), & ! [IN]
1347  kmax, & ! [IN]
1348  z(:), & ! [IN]
1349  temp(:) ) ! [OUT]
1350 
1351  call profile_interp( user_kmax, & ! [IN]
1352  work_z(:), & ! [IN]
1353  work(:), & ! [IN]
1354  kmax+1, & ! [IN]
1355  zh(:), & ! [IN]
1356  temph(:) ) ! [OUT]
1357 
1358  do k = 1, user_kmax
1359  work(k) = user_qv(k) / mvap * mdry / ppm ! [kg/kg->PPM]
1360  enddo
1361 
1362  call profile_interp( user_kmax, & ! [IN]
1363  work_z(:), & ! [IN]
1364  work(:), & ! [IN]
1365  kmax, & ! [IN]
1366  z(:), & ! [IN]
1367  gas(:,1) ) ! [OUT]
1368 
1369  do k = 1, user_kmax
1370  work(k) = user_o3(k) / 48.0_rp * mdry / ppm ! [kg/kg->PPM]
1371  enddo
1372 
1373  call profile_interp( user_kmax, & ! [IN]
1374  work_z(:), & ! [IN]
1375  work(:), & ! [IN]
1376  kmax, & ! [IN]
1377  z(:), & ! [IN]
1378  gas(:,3) ) ! [OUT]
1379 
1380  return
1381  end subroutine profile_read_user
1382 
1383 end module scale_atmos_phy_rd_profile
scale_const::const_ppm
real(rp), parameter, public const_ppm
parts par million
Definition: scale_const.F90:93
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
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:19
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_stdio::io_nml
logical, public io_nml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:62
scale_atmos_phy_rd_profile::atmos_phy_rd_profile_setup
subroutine, public atmos_phy_rd_profile_setup
Setup.
Definition: scale_atmos_phy_rd_profile.F90:139
scale_stdio::io_get_available_fid
integer function, public io_get_available_fid()
search & get available file ID
Definition: scale_stdio.F90:313
scale_grid_index::kmax
integer, public kmax
of computational cells: z, local
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:1128
scale_process
module PROCESS
Definition: scale_process.F90:14
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:561
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:43
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
scale_stdio::io_fid_nml
integer, public io_fid_nml
Log file ID (only for output namelist)
Definition: scale_stdio.F90:57
Generated by   doxygen 1.8.14