SCALE-RM
mod_cnvtopo.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
10 module mod_cnvtopo
11  !-----------------------------------------------------------------------------
12  !
13  !++ used modules
14  !
15  use scale_precision
16  use scale_stdio
17  use scale_prof
18  use scale_grid_index
19  use scale_tracer
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: cnvtopo_setup
28  public :: cnvtopo
29 
30  !-----------------------------------------------------------------------------
31  !
32  !++ Public parameters & variables
33  !
34  logical, public :: cnvtopo_donothing
35  logical, public :: cnvtopo_usegtopo30 = .false.
36  logical, public :: cnvtopo_usegmted2010 = .false.
37  logical, public :: cnvtopo_usedem50m = .false.
38 
39  !-----------------------------------------------------------------------------
40  !
41  !++ Private procedure
42  !
43  private :: cnvtopo_gtopo30
44  private :: cnvtopo_gmted2010
45  private :: cnvtopo_dem50m
46  private :: cnvtopo_smooth
47 
48  !-----------------------------------------------------------------------------
49  !
50  !++ Private parameters & variables
51  !
52  character(len=H_SHORT), private :: cnvtopo_smooth_type = 'LAPLACIAN'
53  integer, private :: cnvtopo_smooth_hypdiff_niter = 20
54  logical, private :: cnvtopo_smooth_local = .true.
55  integer, private :: cnvtopo_smooth_itelim = 10000
56 
57  logical, private :: cnvtopo_copy_parent = .false.
58 
59  real(RP), private :: cnvtopo_unittile_ddeg = 0.0_rp ! dx for unit tile [deg]
60  real(RP), private :: cnvtopo_oversampling_factor = 2.0_rp ! factor of min. dx against the unit tile
61  real(RP), private :: cnvtopo_smooth_maxslope_ratio = 1.0_rp ! ratio of DZDX, DZDY
62  real(RP), private :: cnvtopo_smooth_maxslope = -1.0_rp ! [deg]
63 
64  real(RP), private :: cnvtopo_smooth_maxslope_limit
65 
66  !-----------------------------------------------------------------------------
67 contains
68  !-----------------------------------------------------------------------------
70  subroutine cnvtopo_setup
71  use scale_process, only: &
72  prc_mpistop
73  use scale_const, only: &
74  d2r => const_d2r, &
75  huge => const_huge
76  use scale_comm, only: &
77  comm_horizontal_min
78  use scale_grid, only: &
79  dx, &
80  dy, &
81  grid_cdz, &
82  grid_fdx, &
83  grid_fdy
84  use scale_grid_real, only: &
85  real_dlat, &
86  real_dlon
87  implicit none
88 
89  character(len=H_SHORT) :: cnvtopo_name = 'NONE' ! keep backward compatibility
90 
91  namelist / param_cnvtopo / &
92  cnvtopo_name, &
93  cnvtopo_usegtopo30, &
94  cnvtopo_usegmted2010, &
95  cnvtopo_usedem50m, &
96  cnvtopo_unittile_ddeg, &
97  cnvtopo_oversampling_factor, &
98  cnvtopo_smooth_hypdiff_niter, &
99  cnvtopo_smooth_maxslope_ratio, &
100  cnvtopo_smooth_maxslope, &
101  cnvtopo_smooth_local, &
102  cnvtopo_smooth_itelim, &
103  cnvtopo_smooth_type, &
104  cnvtopo_copy_parent
105 
106  real(RP) :: minslope(ia,ja)
107  real(RP) :: dxl(ia-1)
108  real(RP) :: dyl(ja-1)
109  real(RP) :: dzdx, dzdy
110 
111  real(RP) :: drad(ia,ja)
112  real(RP) :: drad_min
113 
114  integer :: ierr
115  integer :: k, i, j
116  !---------------------------------------------------------------------------
117 
118  if( io_l ) write(io_fid_log,*)
119  if( io_l ) write(io_fid_log,*) '++++++ Module[convert topo] / Categ[preprocess] / Origin[SCALE-RM]'
120 
121  if ( cnvtopo_smooth_local ) then
122  dxl(:) = dx
123  dyl(:) = dy
124  else
125  dxl(:) = grid_fdx(:)
126  dyl(:) = grid_fdy(:)
127  endif
128 
129  !--- read namelist
130  rewind(io_fid_conf)
131  read(io_fid_conf,nml=param_cnvtopo,iostat=ierr)
132  if( ierr < 0 ) then !--- missing
133  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
134  elseif( ierr > 0 ) then !--- fatal error
135  write(*,*) 'xxx Not appropriate names in namelist PARAM_CNVTOPO. Check!'
136  call prc_mpistop
137  endif
138  if( io_nml ) write(io_fid_nml,nml=param_cnvtopo)
139 
140  select case(cnvtopo_name)
141  case('NONE')
142  ! do nothing
143  case('GTOPO30')
144  cnvtopo_usegtopo30 = .true.
145  cnvtopo_usegmted2010 = .false.
146  cnvtopo_usedem50m = .false.
147  case('GMTED2010')
148  cnvtopo_usegtopo30 = .false.
149  cnvtopo_usegmted2010 = .true.
150  cnvtopo_usedem50m = .false.
151  case('DEM50M')
152  cnvtopo_usegtopo30 = .false.
153  cnvtopo_usegmted2010 = .false.
154  cnvtopo_usedem50m = .true.
155  case('COMBINE')
156  cnvtopo_usegtopo30 = .true.
157  cnvtopo_usegmted2010 = .true.
158  cnvtopo_usedem50m = .true.
159  case default
160  write(*,*) 'xxx Unsupported TYPE:', trim(cnvtopo_name)
161  call prc_mpistop
162  endselect
163 
164  cnvtopo_donothing = .true.
165 
166  if ( cnvtopo_usegtopo30 ) then
167  cnvtopo_donothing = .false.
168  if( io_l ) write(io_fid_log,*) '*** Use GTOPO, global 30 arcsec. data'
169  if ( cnvtopo_usegmted2010 ) then
170  if( io_l ) write(io_fid_log,*) '*** Use GMTED2010, new global 5 arcsec. data'
171  if( io_l ) write(io_fid_log,*) '*** Overwrite Existing region'
172  endif
173  if ( cnvtopo_usedem50m ) then
174  if( io_l ) write(io_fid_log,*) '*** Use DEM 50m data for Japan region'
175  if( io_l ) write(io_fid_log,*) '*** Overwrite Japan region'
176  endif
177  elseif ( cnvtopo_usegmted2010 ) then
178  cnvtopo_donothing = .false.
179  if( io_l ) write(io_fid_log,*) '*** Use GMTED2010, new global 5 arcsec. data'
180  if ( cnvtopo_usedem50m ) then
181  if( io_l ) write(io_fid_log,*) '*** Use DEM 50m data for Japan region'
182  if( io_l ) write(io_fid_log,*) '*** Overwrite Japan region'
183  endif
184  elseif ( cnvtopo_usedem50m ) then
185  cnvtopo_donothing = .false.
186  if( io_l ) write(io_fid_log,*) '*** Use DEM 50m data, Japan region only'
187  endif
188 
189  if ( cnvtopo_donothing ) then
190  if( io_l ) write(io_fid_log,*) '*** Do nothing for topography data'
191  else
192  drad(:,:) = min( real_dlat(:,:), real_dlon(:,:) )
193  call comm_horizontal_min( drad_min, drad(:,:) )
194 
195  if ( cnvtopo_unittile_ddeg > 0.0_rp ) then
196  cnvtopo_oversampling_factor = ( drad_min / d2r ) / cnvtopo_unittile_ddeg
197  endif
198  cnvtopo_oversampling_factor = max( 1.0_rp, cnvtopo_oversampling_factor )
199  cnvtopo_unittile_ddeg = ( drad_min / d2r ) / cnvtopo_oversampling_factor
200 
201  if( io_l ) write(io_fid_log,*) '*** The size of tile [deg] = ', cnvtopo_unittile_ddeg
202  if( io_l ) write(io_fid_log,*) '*** oversampling factor = ', cnvtopo_oversampling_factor
203  endif
204 
205  if( cnvtopo_smooth_maxslope > 0.0_rp ) then
206 
207  cnvtopo_smooth_maxslope_limit = cnvtopo_smooth_maxslope
208 
209  else
210  minslope(:,:) = huge
211 
212  j = js-1
213  i = is-1
214  do k = ks, ke
215  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dxl(i) ) / d2r
216  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dyl(j) ) / d2r
217  minslope(is,js) = min( minslope(is,js), dzdx, dzdy )
218  enddo
219 
220  j = js-1
221  do i = is, ie
222  do k = ks, ke
223  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dxl(i) ) / d2r
224  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dyl(j) ) / d2r
225  minslope(i,js) = min( minslope(i,js), dzdx, dzdy )
226  enddo
227  enddo
228 
229  i = is-1
230  do j = js, je
231  do k = ks, ke
232  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dxl(i) ) / d2r
233  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dyl(j) ) / d2r
234  minslope(is,j) = min( minslope(is,j), dzdx, dzdy )
235  enddo
236  enddo
237 
238  do j = js, je
239  do i = is, ie
240  do k = ks, ke
241  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dxl(i) ) / d2r
242  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * grid_cdz(k), dyl(j) ) / d2r
243  minslope(i,j) = min( minslope(i,j), dzdx, dzdy )
244  enddo
245  enddo
246  enddo
247 
248  call comm_horizontal_min( cnvtopo_smooth_maxslope_limit, minslope(:,:) )
249  end if
250 
251  return
252  end subroutine cnvtopo_setup
253 
254  !-----------------------------------------------------------------------------
256  subroutine cnvtopo
257  use scale_process, only: &
258  prc_mpistop
259  use scale_topography, only: &
260  topo_fillhalo, &
261  topo_zsfc, &
262  topo_write
263  use mod_copytopo, only: &
264  copytopo
265  implicit none
266  !---------------------------------------------------------------------------
267 
268  if ( cnvtopo_donothing ) then
269  if( io_l ) write(io_fid_log,*)
270  if( io_l ) write(io_fid_log,*) '++++++ SKIP CONVERT TOPOGRAPHY DATA ++++++'
271  else
272  if( io_l ) write(io_fid_log,*)
273  if( io_l ) write(io_fid_log,*) '++++++ START CONVERT TOPOGRAPHY DATA ++++++'
274 
275  if ( cnvtopo_usegtopo30 ) then
276  call cnvtopo_gtopo30
277  endif
278 
279  if ( cnvtopo_usegmted2010 ) then
280  call cnvtopo_gmted2010
281  endif
282 
283  if ( cnvtopo_usedem50m ) then
284  call cnvtopo_dem50m
285  endif
286 
287  call cnvtopo_smooth( topo_zsfc(:,:) ) ! (inout)
288  call topo_fillhalo
289 
290  if( cnvtopo_copy_parent ) call copytopo( topo_zsfc )
291 
292  if( io_l ) write(io_fid_log,*) '++++++ END CONVERT TOPOGRAPHY DATA ++++++'
293 
294  ! output topography file
295  call topo_write
296  endif
297 
298  return
299  end subroutine cnvtopo
300 
301  !-----------------------------------------------------------------------------
303  subroutine cnvtopo_gtopo30
304  use scale_process, only: &
305  prc_mpistop
306  use scale_const, only: &
307  radius => const_radius, &
308  pi => const_pi, &
309  eps => const_eps, &
310  d2r => const_d2r
311  use scale_topography, only: &
312  topo_zsfc
313  use scale_grid_real, only: &
314  real_laty, &
315  real_lonx
316  implicit none
317 
318  character(len=H_LONG) :: gtopo30_in_catalogue = ''
319  character(len=H_LONG) :: gtopo30_in_dir = ''
320 
321  namelist / param_cnvtopo_gtopo30 / &
322  gtopo30_in_catalogue, &
323  gtopo30_in_dir
324 
325  ! data catalogue list
326  integer, parameter :: tile_nlim = 100
327  integer :: tile_nmax
328  real(RP) :: tile_lats (tile_nlim)
329  real(RP) :: tile_late (tile_nlim)
330  real(RP) :: tile_lons (tile_nlim)
331  real(RP) :: tile_lone (tile_nlim)
332  character(len=H_LONG) :: tile_fname(tile_nlim)
333 
334  ! GTOPO30 data
335  integer, parameter :: isize_orig = 4800
336  integer, parameter :: jsize_orig = 6000
337  integer(2) :: tile_height_orig(isize_orig,jsize_orig)
338  real(RP) :: tile_dlat_orig, tile_dlon_orig
339 
340  ! GTOPO30 data (oversampling)
341  integer :: ios
342  integer :: jos
343  integer :: isize
344  integer :: jsize
345  integer(2), allocatable :: tile_height(:,:)
346  real(RP), allocatable :: tile_lath (:)
347  real(RP), allocatable :: tile_lonh (:)
348  real(RP) :: tile_dlat, tile_dlon
349  real(RP) :: area, area_fraction
350 
351  integer :: iloc
352  integer :: jloc
353  real(RP) :: ifrac_l ! fraction for iloc
354  real(RP) :: jfrac_b ! fraction for jloc
355 
356  real(RP) :: real_lonx_mod(0:ia,ja)
357  real(RP) :: domain_lats, domain_late
358  real(RP) :: domain_lons, domain_lone
359  integer :: domain_lonsloc(2), domain_loneloc(2)
360  logical :: check_idl
361 
362  real(RP) :: topo_sum(ia,ja)
363  real(RP) :: area_sum(ia,ja)
364  real(RP) :: topo, mask
365 
366  character(len=H_LONG) :: fname
367 
368  real(RP) :: zerosw
369  logical :: hit_lat, hit_lon
370  integer :: index
371  integer :: fid, ierr
372  integer :: i, j, ii, jj, iii, jjj, t
373  !---------------------------------------------------------------------------
374 
375  if( io_l ) write(io_fid_log,*)
376  if( io_l ) write(io_fid_log,*) '++++++ Module[convert GTOPO30] / Categ[preprocess] / Origin[SCALE-RM]'
377 
378  !--- read namelist
379  rewind(io_fid_conf)
380  read(io_fid_conf,nml=param_cnvtopo_gtopo30,iostat=ierr)
381  if( ierr < 0 ) then !--- missing
382  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
383  elseif( ierr > 0 ) then !--- fatal error
384  write(*,*) 'xxx Not appropriate names in namelist PARAM_CNVTOPO_GTOPO30. Check!'
385  call prc_mpistop
386  endif
387  if( io_nml ) write(io_fid_nml,nml=param_cnvtopo_gtopo30)
388 
389  do j = 1, ja
390  do i = 1, ia
391  area_sum(i,j) = 0.0_rp
392  topo_sum(i,j) = 0.0_rp
393  enddo
394  enddo
395 
396  real_lonx_mod(:,:) = mod( real_lonx(:,:)+3.0_dp*pi, 2.0_dp*pi ) - pi
397 
398  domain_lats = minval(real_laty(:,:))
399  domain_late = maxval(real_laty(:,:))
400  domain_lons = minval(real_lonx_mod(:,:))
401  domain_lone = maxval(real_lonx_mod(:,:))
402 
403  domain_lonsloc = minloc(real_lonx_mod(:,:))
404  domain_loneloc = maxloc(real_lonx_mod(:,:))
405 
406  check_idl = .false.
407  if ( domain_lons < real_lonx_mod(0 ,domain_lonsloc(2)) &
408  .OR. domain_lone > real_lonx_mod(ia,domain_loneloc(2)) ) then
409  check_idl = .true.
410  domain_lons = minval(real_lonx_mod(:,:),mask=(real_lonx_mod>0.0_rp))
411  domain_lone = maxval(real_lonx_mod(:,:),mask=(real_lonx_mod<0.0_rp))
412  endif
413 
414  jos = nint( 30.0_rp / 60.0_rp / 60.0_rp / cnvtopo_unittile_ddeg - 0.5_rp ) + 1
415  ios = nint( 30.0_rp / 60.0_rp / 60.0_rp / cnvtopo_unittile_ddeg - 0.5_rp ) + 1
416  jsize = jsize_orig * jos
417  isize = isize_orig * ios
418 
419  allocate( tile_height(isize,jsize) )
420  allocate( tile_lath(0:jsize) )
421  allocate( tile_lonh(0:isize) )
422 
423  if( io_l ) write(io_fid_log,*) '*** Oversampling (j) orig = ', jsize_orig, ', use = ', jsize
424  if( io_l ) write(io_fid_log,*) '*** Oversampling (i) orig = ', isize_orig, ', use = ', isize
425 
426  tile_dlat_orig = 30.0_rp / 60.0_rp / 60.0_rp * d2r
427  tile_dlon_orig = 30.0_rp / 60.0_rp / 60.0_rp * d2r
428  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT :', tile_dlat_orig/d2r
429  if( io_l ) write(io_fid_log,*) '*** TILE_DLON :', tile_dlon_orig/d2r
430 
431  tile_dlat = tile_dlat_orig / jos
432  tile_dlon = tile_dlon_orig / ios
433  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT (OS) :', tile_dlat/d2r
434  if( io_l ) write(io_fid_log,*) '*** TILE_DLON (OS) :', tile_dlon/d2r
435 
436  !---< READ from external files >---
437 
438  ! catalogue file
439  fname = trim(gtopo30_in_dir)//'/'//trim(gtopo30_in_catalogue)
440 
441  if( io_l ) write(io_fid_log,*)
442  if( io_l ) write(io_fid_log,*) '+++ Input catalogue file:', trim(fname)
443 
444  fid = io_get_available_fid()
445  open( fid, &
446  file = trim(fname), &
447  form = 'formatted', &
448  status = 'old', &
449  iostat = ierr )
450 
451  if ( ierr /= 0 ) then
452  write(*,*) 'xxx catalogue file not found!', trim(fname)
453  call prc_mpistop
454  endif
455 
456  do t = 1, tile_nlim
457  read(fid,*,iostat=ierr) index, tile_lats(t), tile_late(t), & ! South->North
458  tile_lons(t), tile_lone(t), & ! WEST->EAST
459  tile_fname(t)
460  if ( ierr /= 0 ) exit
461 
462  if ( tile_lons(t) >= 180.0_rp ) then
463  tile_lons(t) = tile_lons(t) - 360.0_rp
464  tile_lone(t) = tile_lone(t) - 360.0_rp
465  endif
466  if ( tile_lons(t) < -180.0_rp ) tile_lons(t) = tile_lons(t) + 360.0_rp
467  if ( tile_lone(t) < -180.0_rp ) tile_lone(t) = tile_lone(t) + 360.0_rp
468 
469  enddo
470 
471  tile_nmax = t - 1
472  close(fid)
473 
474  ! data file
475  do t = 1, tile_nmax
476  hit_lat = .false.
477  hit_lon = .false.
478 
479  if ( ( tile_lats(t)*d2r >= domain_lats .AND. tile_lats(t)*d2r < domain_late ) &
480  .OR. ( tile_late(t)*d2r >= domain_lats .AND. tile_late(t)*d2r < domain_late ) ) then
481  hit_lat = .true.
482  endif
483 
484  if ( ( domain_lats >= tile_lats(t)*d2r .AND. domain_lats < tile_late(t)*d2r ) &
485  .OR. ( domain_late >= tile_lats(t)*d2r .AND. domain_late < tile_late(t)*d2r ) ) then
486  hit_lat = .true.
487  endif
488 
489  if ( check_idl ) then
490  if ( ( tile_lons(t)*d2r >= domain_lons .AND. tile_lons(t)*d2r < pi ) &
491  .OR. ( tile_lons(t)*d2r >= -pi .AND. tile_lons(t)*d2r < domain_lone ) &
492  .OR. ( tile_lone(t)*d2r >= domain_lons .AND. tile_lone(t)*d2r < pi ) &
493  .OR. ( tile_lone(t)*d2r >= -pi .AND. tile_lone(t)*d2r < domain_lone ) ) then
494  hit_lon = .true.
495  endif
496  else
497  if ( ( tile_lons(t)*d2r >= domain_lons .AND. tile_lons(t)*d2r < domain_lone ) &
498  .OR. ( tile_lone(t)*d2r >= domain_lons .AND. tile_lone(t)*d2r < domain_lone ) ) then
499  hit_lon = .true.
500  endif
501  endif
502 
503  if ( ( domain_lons >= tile_lons(t)*d2r .AND. domain_lons < tile_lone(t)*d2r ) &
504  .OR. ( domain_lone >= tile_lons(t)*d2r .AND. domain_lone < tile_lone(t)*d2r ) ) then
505  hit_lon = .true.
506  endif
507 
508  if ( hit_lat .AND. hit_lon ) then
509  fname = trim(gtopo30_in_dir)//'/'//trim(tile_fname(t))
510 
511  if( io_l ) write(io_fid_log,*)
512  if( io_l ) write(io_fid_log,*) '+++ Input data file :', trim(fname)
513  if( io_l ) write(io_fid_log,*) '*** Domain (LAT) :', domain_lats/d2r, domain_late/d2r
514  if( io_l ) write(io_fid_log,*) '*** (LON) :', domain_lons/d2r, domain_lone/d2r
515  if ( check_idl ) then
516  if( io_l ) write(io_fid_log,*) '*** (Date line exists within the domain)'
517  endif
518  if( io_l ) write(io_fid_log,*) '*** Tile (LAT) :', tile_lats(t), tile_late(t)
519  if( io_l ) write(io_fid_log,*) '*** (LON) :', tile_lons(t), tile_lone(t)
520 
521  fid = io_get_available_fid()
522  open( fid, &
523  file = trim(fname), &
524  form = 'unformatted', &
525  access = 'direct', &
526  status = 'old', &
527  recl = isize_orig*jsize_orig*2, &
528  iostat = ierr )
529 
530  if ( ierr /= 0 ) then
531  write(*,*) 'xxx data file not found!'
532  call prc_mpistop
533  endif
534 
535  read(fid,rec=1) tile_height_orig(:,:)
536  close(fid)
537 
538  ! oversampling
539  do jj = 1, jsize_orig
540  do ii = 1, isize_orig
541  do j = 1, jos
542  do i = 1, ios
543  jjj = (jj-1) * jos + j
544  iii = (ii-1) * ios + i
545 
546  tile_height(iii,jjj) = tile_height_orig(ii,jsize_orig-jj+1) ! reverse y-axis
547  enddo
548  enddo
549  enddo
550  enddo
551 
552  tile_lath(0) = tile_lats(t) * d2r
553  do jj = 1, jsize
554  tile_lath(jj) = tile_lath(jj-1) + tile_dlat
555 ! if( IO_L ) write(IO_FID_LOG,*) jj, TILE_LATH(jj)
556  enddo
557 
558  tile_lonh(0) = tile_lons(t) * d2r
559  do ii = 1, isize
560  tile_lonh(ii) = tile_lonh(ii-1) + tile_dlon
561 ! if( IO_L ) write(IO_FID_LOG,*) ii, TILE_LONH(ii)
562  enddo
563 
564  ! match and calc fraction
565  do jj = 1, jsize
566  do ii = 1, isize
567 
568  iloc = 1 ! Z_sfc(1,1) is used for dummy grid
569  ifrac_l = 1.0_rp
570 
571  jloc = 1 ! Z_sfc(1,1) is used for dummy grid
572  jfrac_b = 1.0_rp
573 
574  if ( tile_lath(jj ) < domain_lats &
575  .OR. tile_lath(jj-1) > domain_late ) then
576  cycle
577  endif
578 
579  if ( check_idl ) then
580  if ( tile_lonh(ii ) < domain_lons &
581  .AND. tile_lonh(ii-1) > domain_lone ) then
582  cycle
583  endif
584  else
585  if ( tile_lonh(ii ) < domain_lons &
586  .OR. tile_lonh(ii-1) > domain_lone ) then
587  cycle
588  endif
589  endif
590 
591  jloop: do j = js-1, je+1
592  iloop: do i = is-1, ie+1
593  if ( tile_lonh(ii-1) >= real_lonx_mod(i-1,j ) &
594  .AND. tile_lonh(ii-1) < real_lonx_mod(i ,j ) &
595  .AND. tile_lath(jj-1) >= real_laty(i ,j-1) &
596  .AND. tile_lath(jj-1) < real_laty(i ,j ) ) then
597 
598  iloc = i
599  ifrac_l = min( real_lonx_mod(i,j)-tile_lonh(ii-1), tile_dlon ) / tile_dlon
600 
601  jloc = j
602  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
603  exit jloop
604 
605  endif
606 
607  if ( real_lonx_mod(i-1,j) >= real_lonx_mod(i ,j ) &
608  .AND. tile_lath(jj-1) >= real_laty(i ,j-1) &
609  .AND. tile_lath(jj-1) < real_laty(i ,j ) ) then ! across the IDL
610 
611  if ( tile_lonh(ii-1) >= real_lonx_mod(i-1,j) &
612  .AND. tile_lonh(ii-1) < pi ) then
613 
614  iloc = i
615  ifrac_l = min( real_lonx_mod(i,j)-tile_lonh(ii-1)+2.0_rp*pi, tile_dlon ) / tile_dlon
616 
617  jloc = j
618  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
619  exit jloop
620 
621  elseif( tile_lonh(ii-1) >= -pi &
622  .AND. tile_lonh(ii-1) < real_lonx_mod(i ,j) ) then
623 
624  iloc = i
625  ifrac_l = min( real_lonx_mod(i,j)-tile_lonh(ii-1), tile_dlon ) / tile_dlon
626 
627  jloc = j
628  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
629  exit jloop
630 
631  endif
632 
633  endif
634  enddo iloop
635  enddo jloop
636 
637  if( iloc == 1 .AND. jloc == 1 ) cycle
638 
639  topo = real( TILE_HEIGHT(ii,jj), kind=rp )
640  mask = 0.5_rp - sign( 0.5_rp, topo ) ! if Height is negative, mask = 1
641 
642  area = radius * radius * tile_dlon * ( sin(tile_lath(jj))-sin(tile_lath(jj-1)) ) * ( 1.0_rp - mask )
643 
644 ! if( IO_L ) write(IO_FID_LOG,*) ii, jj, area, iloc, jloc, ifrac_l, jfrac_b, TILE_HEIGHT(ii,jj)
645 
646  area_fraction = ( ifrac_l) * ( jfrac_b) * area
647  area_sum(iloc ,jloc ) = area_sum(iloc ,jloc ) + area_fraction
648  topo_sum(iloc ,jloc ) = topo_sum(iloc ,jloc ) + area_fraction * topo
649 
650  area_fraction = (1.0_rp-ifrac_l) * ( jfrac_b) * area
651  area_sum(iloc+1,jloc ) = area_sum(iloc+1,jloc ) + area_fraction
652  topo_sum(iloc+1,jloc ) = topo_sum(iloc+1,jloc ) + area_fraction * topo
653 
654  area_fraction = ( ifrac_l) * (1.0_rp-jfrac_b) * area
655  area_sum(iloc ,jloc+1) = area_sum(iloc ,jloc+1) + area_fraction
656  topo_sum(iloc ,jloc+1) = topo_sum(iloc ,jloc+1) + area_fraction * topo
657 
658  area_fraction = (1.0_rp-ifrac_l) * (1.0_rp-jfrac_b) * area
659  area_sum(iloc+1,jloc+1) = area_sum(iloc+1,jloc+1) + area_fraction
660  topo_sum(iloc+1,jloc+1) = topo_sum(iloc+1,jloc+1) + area_fraction * topo
661  enddo
662  enddo
663 
664  endif
665  enddo ! tile loop
666 
667  do j = js, je
668  do i = is, ie
669  zerosw = 0.5_rp - sign( 0.5_rp, area_sum(i,j)-eps )
670  topo_zsfc(i,j) = topo_sum(i,j) * ( 1.0_rp-zerosw ) / ( area_sum(i,j)-zerosw )
671  enddo
672  enddo
673 
674  return
675  end subroutine cnvtopo_gtopo30
676 
677  !-----------------------------------------------------------------------------
679  subroutine cnvtopo_gmted2010
680  implicit none
681  !---------------------------------------------------------------------------
682 
683  return
684  end subroutine cnvtopo_gmted2010
685 
686  !-----------------------------------------------------------------------------
688  subroutine cnvtopo_dem50m
689  use scale_process, only: &
690  prc_mpistop
691  use scale_const, only: &
692  radius => const_radius, &
693  pi => const_pi, &
694  eps => const_eps, &
695  d2r => const_d2r
696  use scale_topography, only: &
697  topo_zsfc
698  use scale_grid_real, only: &
699  real_laty, &
700  real_lonx
701  implicit none
702 
703  character(len=H_LONG) :: dem50m_in_catalogue = ''
704  character(len=H_LONG) :: dem50m_in_dir = ''
705 
706  namelist / param_cnvtopo_dem50m / &
707  dem50m_in_catalogue, &
708  dem50m_in_dir
709 
710  ! data catalogue list
711  integer, parameter :: tile_nlim = 1000
712  integer :: tile_nmax
713  real(RP) :: tile_lats (tile_nlim)
714  real(RP) :: tile_late (tile_nlim)
715  real(RP) :: tile_lons (tile_nlim)
716  real(RP) :: tile_lone (tile_nlim)
717  character(len=H_LONG) :: tile_fname(tile_nlim)
718 
719  ! DEM50M data
720  integer, parameter :: isize_orig = 1600
721  integer, parameter :: jsize_orig = 1600
722  real(SP) :: tile_height_orig(isize_orig,jsize_orig)
723  real(RP) :: tile_dlat_orig, tile_dlon_orig
724 
725  ! DEM50M data (oversampling)
726  integer :: ios
727  integer :: jos
728  integer :: isize
729  integer :: jsize
730  real(SP), allocatable :: tile_height(:,:)
731  real(RP), allocatable :: tile_lath (:)
732  real(RP), allocatable :: tile_lonh (:)
733  real(RP) :: tile_dlat, tile_dlon
734  real(RP) :: area, area_fraction
735 
736  integer :: iloc
737  integer :: jloc
738  real(RP) :: ifrac_l ! fraction for iloc
739  real(RP) :: jfrac_b ! fraction for jloc
740 
741  real(RP) :: real_lonx_mod(0:ia,ja)
742  real(RP) :: domain_lats, domain_late
743  real(RP) :: domain_lons, domain_lone
744  integer :: domain_lonsloc(2), domain_loneloc(2)
745  logical :: check_idl
746 
747  real(RP) :: topo_sum(ia,ja)
748  real(RP) :: area_sum(ia,ja)
749  real(RP) :: topo, mask
750 
751  character(len=H_LONG) :: fname
752 
753  real(RP) :: zerosw
754  logical :: hit_lat, hit_lon
755  integer :: index
756  integer :: fid, ierr
757  integer :: i, j, ii, jj, iii, jjj, t
758  !---------------------------------------------------------------------------
759 
760  if( io_l ) write(io_fid_log,*)
761  if( io_l ) write(io_fid_log,*) '++++++ Module[convert DEM50M] / Categ[preprocess] / Origin[SCALE-RM]'
762 
763  !--- read namelist
764  rewind(io_fid_conf)
765  read(io_fid_conf,nml=param_cnvtopo_dem50m,iostat=ierr)
766  if( ierr < 0 ) then !--- missing
767  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
768  elseif( ierr > 0 ) then !--- fatal error
769  write(*,*) 'xxx Not appropriate names in namelist PARAM_CNVTOPO_DEM50M. Check!'
770  call prc_mpistop
771  endif
772  if( io_nml ) write(io_fid_nml,nml=param_cnvtopo_dem50m)
773 
774  do j = 1, ja
775  do i = 1, ia
776  area_sum(i,j) = 0.0_rp
777  topo_sum(i,j) = 0.0_rp
778  enddo
779  enddo
780 
781  real_lonx_mod(:,:) = mod( real_lonx(:,:)+3.0_dp*pi, 2.0_dp*pi ) - pi
782 
783  domain_lats = minval(real_laty(:,:))
784  domain_late = maxval(real_laty(:,:))
785  domain_lons = minval(real_lonx_mod(:,:))
786  domain_lone = maxval(real_lonx_mod(:,:))
787 
788  domain_lonsloc = minloc(real_lonx_mod(:,:))
789  domain_loneloc = maxloc(real_lonx_mod(:,:))
790 
791  check_idl = .false.
792  if ( domain_lons < real_lonx_mod(0 ,domain_lonsloc(2)) &
793  .OR. domain_lone > real_lonx_mod(ia,domain_loneloc(2)) ) then
794  check_idl = .true.
795  domain_lons = minval(real_lonx_mod(:,:),mask=(real_lonx_mod>0.0_rp))
796  domain_lone = maxval(real_lonx_mod(:,:),mask=(real_lonx_mod<0.0_rp))
797  endif
798 
799  jos = nint( 5.0_rp / 60.0_rp / 200.0_rp / cnvtopo_unittile_ddeg - 0.5_rp ) + 1
800  ios = nint( 7.5_rp / 60.0_rp / 200.0_rp / cnvtopo_unittile_ddeg - 0.5_rp ) + 1
801  jsize = jsize_orig * jos
802  isize = isize_orig * ios
803 
804  allocate( tile_height(isize,jsize) )
805  allocate( tile_lath(0:jsize) )
806  allocate( tile_lonh(0:isize) )
807 
808  if( io_l ) write(io_fid_log,*) '*** Oversampling (j) orig = ', jsize_orig, ', use = ', jsize
809  if( io_l ) write(io_fid_log,*) '*** Oversampling (i) orig = ', isize_orig, ', use = ', isize
810 
811  tile_dlat_orig = 5.0_rp / 60.0_rp / 200.0_rp * d2r
812  tile_dlon_orig = 7.5_rp / 60.0_rp / 200.0_rp * d2r
813  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT :', tile_dlat_orig/d2r
814  if( io_l ) write(io_fid_log,*) '*** TILE_DLON :', tile_dlon_orig/d2r
815 
816  tile_dlat = tile_dlat_orig / jos
817  tile_dlon = tile_dlon_orig / ios
818  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT (OS) :', tile_dlat/d2r
819  if( io_l ) write(io_fid_log,*) '*** TILE_DLON (OS) :', tile_dlon/d2r
820 
821  !---< READ from external files >---
822 
823  ! catalogue file
824  fname = trim(dem50m_in_dir)//'/'//trim(dem50m_in_catalogue)
825 
826  if( io_l ) write(io_fid_log,*)
827  if( io_l ) write(io_fid_log,*) '+++ Input catalogue file:', trim(fname)
828 
829  fid = io_get_available_fid()
830  open( fid, &
831  file = trim(fname), &
832  form = 'formatted', &
833  status = 'old', &
834  iostat = ierr )
835 
836  if ( ierr /= 0 ) then
837  write(*,*) 'xxx catalogue file not found!', trim(fname)
838  call prc_mpistop
839  endif
840 
841  do t = 1, tile_nlim
842  read(fid,*,iostat=ierr) index, tile_lats(t), tile_late(t), & ! South->North
843  tile_lons(t), tile_lone(t), & ! WEST->EAST
844  tile_fname(t)
845  if ( ierr /= 0 ) exit
846 
847  if ( tile_lons(t) >= 180.0_rp ) then
848  tile_lons(t) = tile_lons(t) - 360.0_rp
849  tile_lone(t) = tile_lone(t) - 360.0_rp
850  endif
851  if ( tile_lons(t) < -180.0_rp ) tile_lons(t) = tile_lons(t) + 360.0_rp
852  if ( tile_lone(t) < -180.0_rp ) tile_lone(t) = tile_lone(t) + 360.0_rp
853 
854  enddo
855 
856  tile_nmax = t - 1
857  close(fid)
858 
859  ! data file
860  do t = 1, tile_nmax
861  hit_lat = .false.
862  hit_lon = .false.
863 
864  if ( ( tile_lats(t)*d2r >= domain_lats .AND. tile_lats(t)*d2r < domain_late ) &
865  .OR. ( tile_late(t)*d2r >= domain_lats .AND. tile_late(t)*d2r < domain_late ) ) then
866  hit_lat = .true.
867  endif
868 
869  if ( ( domain_lats >= tile_lats(t)*d2r .AND. domain_lats < tile_late(t)*d2r ) &
870  .OR. ( domain_late >= tile_lats(t)*d2r .AND. domain_late < tile_late(t)*d2r ) ) then
871  hit_lat = .true.
872  endif
873 
874  if ( check_idl ) then
875  if ( ( tile_lons(t)*d2r >= domain_lons .AND. tile_lons(t)*d2r < pi ) &
876  .OR. ( tile_lons(t)*d2r >= -pi .AND. tile_lons(t)*d2r < domain_lone ) &
877  .OR. ( tile_lone(t)*d2r >= domain_lons .AND. tile_lone(t)*d2r < pi ) &
878  .OR. ( tile_lone(t)*d2r >= -pi .AND. tile_lone(t)*d2r < domain_lone ) ) then
879  hit_lon = .true.
880  endif
881  else
882  if ( ( tile_lons(t)*d2r >= domain_lons .AND. tile_lons(t)*d2r < domain_lone ) &
883  .OR. ( tile_lone(t)*d2r >= domain_lons .AND. tile_lone(t)*d2r < domain_lone ) ) then
884  hit_lon = .true.
885  endif
886  endif
887 
888  if ( ( domain_lons >= tile_lons(t)*d2r .AND. domain_lons < tile_lone(t)*d2r ) &
889  .OR. ( domain_lone >= tile_lons(t)*d2r .AND. domain_lone < tile_lone(t)*d2r ) ) then
890  hit_lon = .true.
891  endif
892 
893  if ( hit_lat .AND. hit_lon ) then
894  fname = trim(dem50m_in_dir)//'/'//trim(tile_fname(t))
895 
896  if( io_l ) write(io_fid_log,*)
897  if( io_l ) write(io_fid_log,*) '+++ Input data file :', trim(fname)
898  if( io_l ) write(io_fid_log,*) '*** Domain (LAT) :', domain_lats/d2r, domain_late/d2r
899  if( io_l ) write(io_fid_log,*) '*** (LON) :', domain_lons/d2r, domain_lone/d2r
900  if ( check_idl ) then
901  if( io_l ) write(io_fid_log,*) '*** (Date line exists within the domain)'
902  endif
903  if( io_l ) write(io_fid_log,*) '*** Tile (LAT) :', tile_lats(t), tile_late(t)
904  if( io_l ) write(io_fid_log,*) '*** (LON) :', tile_lons(t), tile_lone(t)
905 
906  fid = io_get_available_fid()
907  open( fid, &
908  file = trim(fname), &
909  form = 'unformatted', &
910  access = 'direct', &
911  status = 'old', &
912  recl = isize_orig*jsize_orig*4, &
913  iostat = ierr )
914 
915  if ( ierr /= 0 ) then
916  write(*,*) 'xxx data file not found!'
917  call prc_mpistop
918  endif
919 
920  read(fid,rec=1) tile_height_orig(:,:)
921  close(fid)
922 
923  ! oversampling
924  do jj = 1, jsize_orig
925  do ii = 1, isize_orig
926  do j = 1, jos
927  do i = 1, ios
928  jjj = (jj-1) * jos + j
929  iii = (ii-1) * ios + i
930 
931  tile_height(iii,jjj) = tile_height_orig(ii,jj)
932  enddo
933  enddo
934  enddo
935  enddo
936 
937  tile_lath(0) = tile_lats(t) * d2r
938  do jj = 1, jsize
939  tile_lath(jj) = tile_lath(jj-1) + tile_dlat
940 ! if( IO_L ) write(IO_FID_LOG,*) jj, TILE_LATH(jj)
941  enddo
942 
943  tile_lonh(0) = tile_lons(t) * d2r
944  do ii = 1, isize
945  tile_lonh(ii) = tile_lonh(ii-1) + tile_dlon
946 ! if( IO_L ) write(IO_FID_LOG,*) ii, TILE_LONH(ii)
947  enddo
948 
949  ! match and calc fraction
950  do jj = 1, jsize
951  do ii = 1, isize
952 
953  iloc = 1 ! Z_sfc(1,1) is used for dummy grid
954  ifrac_l = 1.0_rp
955 
956  jloc = 1 ! Z_sfc(1,1) is used for dummy grid
957  jfrac_b = 1.0_rp
958 
959  if ( tile_lath(jj ) < domain_lats &
960  .OR. tile_lath(jj-1) > domain_late ) then
961  cycle
962  endif
963 
964  if ( check_idl ) then
965  if ( tile_lonh(ii ) < domain_lons &
966  .AND. tile_lonh(ii-1) > domain_lone ) then
967  cycle
968  endif
969  else
970  if ( tile_lonh(ii ) < domain_lons &
971  .OR. tile_lonh(ii-1) > domain_lone ) then
972  cycle
973  endif
974  endif
975 
976  jloop: do j = js-1, je+1
977  iloop: do i = is-1, ie+1
978  if ( tile_lonh(ii-1) >= real_lonx_mod(i-1,j ) &
979  .AND. tile_lonh(ii-1) < real_lonx_mod(i ,j ) &
980  .AND. tile_lath(jj-1) >= real_laty(i ,j-1) &
981  .AND. tile_lath(jj-1) < real_laty(i ,j ) ) then
982 
983  iloc = i
984  ifrac_l = min( real_lonx_mod(i,j)-tile_lonh(ii-1), tile_dlon ) / tile_dlon
985 
986  jloc = j
987  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
988  exit jloop
989 
990  endif
991 
992  if ( real_lonx_mod(i-1,j) >= real_lonx_mod(i ,j ) &
993  .AND. tile_lath(jj-1) >= real_laty(i ,j-1) &
994  .AND. tile_lath(jj-1) < real_laty(i ,j ) ) then ! across the IDL
995 
996  if ( tile_lonh(ii-1) >= real_lonx_mod(i-1,j) &
997  .AND. tile_lonh(ii-1) < pi ) then
998 
999  iloc = i
1000  ifrac_l = min( real_lonx_mod(i,j)-tile_lonh(ii-1)+2.0_rp*pi, tile_dlon ) / tile_dlon
1001 
1002  jloc = j
1003  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
1004  exit jloop
1005 
1006  elseif( tile_lonh(ii-1) >= -pi &
1007  .AND. tile_lonh(ii-1) < real_lonx_mod(i ,j) ) then
1008 
1009  iloc = i
1010  ifrac_l = min( real_lonx_mod(i,j)-tile_lonh(ii-1), tile_dlon ) / tile_dlon
1011 
1012  jloc = j
1013  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
1014  exit jloop
1015 
1016  endif
1017 
1018  endif
1019  enddo iloop
1020  enddo jloop
1021 
1022  if( iloc == 1 .AND. jloc == 1 ) cycle
1023 
1024  topo = real( TILE_HEIGHT(ii,jj), kind=rp )
1025  mask = 0.5_rp - sign( 0.5_rp,topo ) ! if Height is negative, mask = 1
1026 
1027  area = radius * radius * tile_dlon * ( sin(tile_lath(jj))-sin(tile_lath(jj-1)) ) * ( 1.0_rp - mask )
1028 
1029 ! if( IO_L ) write(IO_FID_LOG,*) ii, jj, area, iloc, jloc, ifrac_l, jfrac_b, TILE_HEIGHT(ii,jj)
1030 
1031  area_fraction = ( ifrac_l) * ( jfrac_b) * area
1032  area_sum(iloc ,jloc ) = area_sum(iloc ,jloc ) + area_fraction
1033  topo_sum(iloc ,jloc ) = topo_sum(iloc ,jloc ) + area_fraction * topo
1034 
1035  area_fraction = (1.0_rp-ifrac_l) * ( jfrac_b) * area
1036  area_sum(iloc+1,jloc ) = area_sum(iloc+1,jloc ) + area_fraction
1037  topo_sum(iloc+1,jloc ) = topo_sum(iloc+1,jloc ) + area_fraction * topo
1038 
1039  area_fraction = ( ifrac_l) * (1.0_rp-jfrac_b) * area
1040  area_sum(iloc ,jloc+1) = area_sum(iloc ,jloc+1) + area_fraction
1041  topo_sum(iloc ,jloc+1) = topo_sum(iloc ,jloc+1) + area_fraction * topo
1042 
1043  area_fraction = (1.0_rp-ifrac_l) * (1.0_rp-jfrac_b) * area
1044  area_sum(iloc+1,jloc+1) = area_sum(iloc+1,jloc+1) + area_fraction
1045  topo_sum(iloc+1,jloc+1) = topo_sum(iloc+1,jloc+1) + area_fraction * topo
1046  enddo
1047  enddo
1048 
1049  endif
1050  enddo ! tile loop
1051 
1052  do j = js, je
1053  do i = is, ie
1054  mask = 0.5_rp + sign( 0.5_rp, area_sum(i,j)-eps ) ! if any data is found, overwrite
1055  zerosw = 0.5_rp - sign( 0.5_rp, area_sum(i,j)-eps )
1056  topo = topo_sum(i,j) * ( 1.0_rp-zerosw ) / ( area_sum(i,j)-zerosw )
1057  topo_zsfc(i,j) = ( mask ) * topo & ! overwrite
1058  + ( 1.0_rp-mask ) * topo_zsfc(i,j) ! keep existing value
1059  enddo
1060  enddo
1061 
1062  return
1063  end subroutine cnvtopo_dem50m
1064 
1065  !-----------------------------------------------------------------------------
1067  subroutine cnvtopo_smooth( &
1068  Zsfc )
1069  use scale_const, only: &
1070  d2r => const_d2r
1071  use scale_process, only: &
1072  prc_mpistop
1073  use scale_grid, only: &
1074  dx, &
1075  dy, &
1076  grid_fdx, &
1077  grid_fdy
1078  use scale_comm, only: &
1079  comm_horizontal_max
1080  use scale_rm_statistics, only: &
1081  stat_detail
1082  use scale_topography, only: &
1083  topo_fillhalo
1084  implicit none
1085 
1086  real(RP), intent(inout) :: zsfc(ia,ja)
1087 
1088  real(RP) :: dzsfc_dx(1,ia,ja,1) ! d(Zsfc)/dx at u-position
1089  real(RP) :: dzsfc_dy(1,ia,ja,1) ! d(Zsfc)/dy at v-position
1090 
1091  real(RP) :: dxl(ia-1)
1092  real(RP) :: dyl(ja-1)
1093 
1094  real(RP) :: flx_x(ia,ja)
1095  real(RP) :: flx_y(ia,ja)
1096 
1097  real(RP) :: slope(ia,ja)
1098  real(RP) :: maxslope
1099  real(RP) :: flag
1100 
1101  character(len=H_SHORT) :: varname(1)
1102 
1103  integer :: ite
1104  integer :: i, j
1105  !---------------------------------------------------------------------------
1106 
1107  if( io_l ) write(io_fid_log,*)
1108  if( io_l ) write(io_fid_log,*) '*** Apply smoothing. Slope limit = ', cnvtopo_smooth_maxslope_limit
1109  if( io_l ) write(io_fid_log,*) '*** Smoothing type = ', cnvtopo_smooth_type
1110  if( io_l ) write(io_fid_log,*) '*** Smoothing locally = ', cnvtopo_smooth_local
1111 
1112  if ( cnvtopo_smooth_local ) then
1113  dxl(:) = dx
1114  dyl(:) = dy
1115  else
1116  dxl(:) = grid_fdx(:)
1117  dyl(:) = grid_fdy(:)
1118  endif
1119 
1120 
1121 
1122  ! digital filter
1123  do ite = 1, cnvtopo_smooth_itelim+1
1124  if( io_l ) write(io_fid_log,*)
1125  if( io_l ) write(io_fid_log,*) '*** Smoothing itelation : ', ite
1126 
1127  call topo_fillhalo( zsfc )
1128 
1129  do j = 1, ja
1130  do i = 1, ia-1
1131  dzsfc_dx(1,i,j,1) = atan2( ( zsfc(i+1,j)-zsfc(i,j) ), dxl(i) ) / d2r
1132  enddo
1133  enddo
1134  dzsfc_dx(1,ia,:,1) = 0.0_rp
1135  do j = 1, ja-1
1136  do i = 1, ia
1137  dzsfc_dy(1,i,j,1) = atan2( ( zsfc(i,j+1)-zsfc(i,j) ), dyl(j) ) / d2r
1138  enddo
1139  enddo
1140  dzsfc_dy(1,:,ja,1) = 0.0_rp
1141 
1142  slope(:,:) = max( abs(dzsfc_dx(1,:,:,1)), abs(dzsfc_dy(1,:,:,1)) )
1143  call comm_horizontal_max( maxslope, slope(:,:) )
1144 
1145  if( io_l ) write(io_fid_log,*) '*** maximum slope [deg] : ', maxslope
1146 
1147  if( maxslope < cnvtopo_smooth_maxslope_limit ) exit
1148 
1149  varname(1) = "DZsfc_DX"
1150  call stat_detail( dzsfc_dx(:,:,:,:), varname(:) )
1151  varname(1) = "DZsfc_DY"
1152  call stat_detail( dzsfc_dy(:,:,:,:), varname(:) )
1153 
1154  select case( cnvtopo_smooth_type )
1155  case( 'GAUSSIAN' )
1156 
1157  ! 3 by 3 gaussian filter
1158  do j = js, je
1159  do i = is, ie
1160  zsfc(i,j) = ( 0.2500_rp * zsfc(i ,j ) &
1161  + 0.1250_rp * zsfc(i-1,j ) &
1162  + 0.1250_rp * zsfc(i+1,j ) &
1163  + 0.1250_rp * zsfc(i ,j-1) &
1164  + 0.1250_rp * zsfc(i ,j+1) &
1165  + 0.0625_rp * zsfc(i-1,j-1) &
1166  + 0.0625_rp * zsfc(i+1,j-1) &
1167  + 0.0625_rp * zsfc(i-1,j+1) &
1168  + 0.0625_rp * zsfc(i+1,j+1) )
1169  enddo
1170  enddo
1171 
1172  case( 'LAPLACIAN' )
1173 
1174  do j = js , je
1175  do i = is-1, ie
1176  flx_x(i,j) = zsfc(i+1,j) - zsfc(i,j)
1177 ! FLX_TMP(i,j) = Zsfc(i+1,j) - Zsfc(i,j)
1178  enddo
1179  enddo
1180 !!$ call TOPO_fillhalo( FLX_TMP )
1181 !!$ do j = JS , JE
1182 !!$ do i = IS-1, IE
1183 !!$ FLX_X(i,j) = - ( FLX_TMP(i+1,j) - FLX_TMP(i,j) )
1184 !!$ enddo
1185 !!$ enddo
1186 
1187  do j = js-1, je
1188  do i = is , ie
1189  flx_y(i,j) = zsfc(i,j+1) - zsfc(i,j)
1190 ! FLX_TMP(i,j) = Zsfc(i,j+1) - Zsfc(i,j)
1191  enddo
1192  enddo
1193 !!$ call TOPO_fillhalo( FLX_TMP )
1194 !!$ do j = JS-1, JE
1195 !!$ do i = IS , IE
1196 !!$ FLX_Y(i,j) = - ( FLX_TMP(i,j+1) - FLX_TMP(i,j) )
1197 !!$ enddo
1198 !!$ enddo
1199 
1200 
1201  if ( cnvtopo_smooth_local ) then
1202  do j = js , je
1203  do i = is-1, ie
1204  flag = 0.5_rp &
1205  + sign(0.5_rp, max( abs(dzsfc_dx(1,i+1,j ,1)), &
1206  abs(dzsfc_dx(1,i ,j ,1)), &
1207  abs(dzsfc_dx(1,i-1,j ,1)), &
1208  abs(dzsfc_dy(1,i+1,j ,1)), &
1209  abs(dzsfc_dy(1,i+1,j-1,1)), &
1210  abs(dzsfc_dy(1,i ,j ,1)), &
1211  abs(dzsfc_dy(1,i ,j-1,1)) &
1212  ) - cnvtopo_smooth_maxslope_limit )
1213  flx_x(i,j) = flx_x(i,j) &
1214  * dxl(i) / grid_fdx(i) &
1215  * flag
1216  enddo
1217  enddo
1218  do j = js-1, je
1219  do i = is , ie
1220  flag = 0.5_rp &
1221  + sign(0.5_rp, max( abs(dzsfc_dy(1,i ,j+1,1)), &
1222  abs(dzsfc_dy(1,i ,j ,1)), &
1223  abs(dzsfc_dy(1,i ,j-1,1)), &
1224  abs(dzsfc_dx(1,i ,j+1,1)), &
1225  abs(dzsfc_dx(1,i-1,j+1,1)), &
1226  abs(dzsfc_dx(1,i ,j ,1)), &
1227  abs(dzsfc_dx(1,i-1,j ,1)) &
1228  ) - cnvtopo_smooth_maxslope_limit )
1229  flx_y(i,j) = flx_y(i,j) &
1230  * dyl(j) / grid_fdy(j) &
1231  * flag
1232  enddo
1233  enddo
1234  endif
1235 
1236  do j = js, je
1237  do i = is, ie
1238  zsfc(i,j) = max( 0.0_rp, &
1239  zsfc(i,j) &
1240  + 0.1_rp * ( ( flx_x(i,j) - flx_x(i-1,j) ) &
1241  + ( flx_y(i,j) - flx_y(i,j-1) ) ) )
1242  enddo
1243  enddo
1244 
1245  case default
1246  write(*,*) 'xxx Invalid smoothing type'
1247  call prc_mpistop
1248  end select
1249 
1250  enddo
1251 
1252  if ( ite > cnvtopo_smooth_itelim ) then
1253  write(*,*) 'xxx not converged'
1254  call prc_mpistop
1255  else
1256  if( io_l ) write(io_fid_log,*) '*** smoothing complete.'
1257  endif
1258 
1259 
1260 
1261  if ( cnvtopo_smooth_hypdiff_niter > 0 ) then
1262 
1263  if( io_l ) write(io_fid_log,*)
1264  if( io_l ) write(io_fid_log,*) '*** Apply hyperdiffusion.'
1265 
1266  call cnvtopo_hypdiff( zsfc(:,:), cnvtopo_smooth_hypdiff_niter )
1267 
1268  do j = 1, ja
1269  do i = 1, ia-1
1270  dzsfc_dx(1,i,j,1) = atan2( ( zsfc(i+1,j)-zsfc(i,j) ), dxl(i) ) / d2r
1271  enddo
1272  enddo
1273  dzsfc_dx(1,ia,:,1) = 0.0_rp
1274  do j = 1, ja-1
1275  do i = 1, ia
1276  dzsfc_dy(1,i,j,1) = atan2( ( zsfc(i,j+1)-zsfc(i,j) ), dyl(j) ) / d2r
1277  enddo
1278  enddo
1279  dzsfc_dy(1,:,ja,1) = 0.0_rp
1280 
1281  slope(:,:) = max( abs(dzsfc_dx(1,:,:,1)), abs(dzsfc_dy(1,:,:,1)) )
1282  call comm_horizontal_max( maxslope, slope(:,:) )
1283 
1284  if( io_l ) write(io_fid_log,*) '*** maximum slope [deg] : ', maxslope
1285 
1286  end if
1287 
1288  call topo_fillhalo( zsfc )
1289 
1290  varname(1) = "DZsfc_DX"
1291  call stat_detail( dzsfc_dx(:,:,:,:), varname(:) )
1292  varname(1) = "DZsfc_DY"
1293  call stat_detail( dzsfc_dy(:,:,:,:), varname(:) )
1294 
1295  if( io_l ) write(io_fid_log,*)
1296 
1297  return
1298  end subroutine cnvtopo_smooth
1299 
1300  subroutine cnvtopo_hypdiff( Zsfc, nite )
1301  use scale_topography, only: &
1302  topo_fillhalo
1303  real(RP), intent(inout) :: Zsfc(IA,JA)
1304  integer, intent(in) :: nite
1305 
1306  real(RP), pointer :: p1(:,:)
1307  real(RP), pointer :: p2(:,:)
1308  real(RP), target :: work1(IA,JA)
1309  real(RP), target :: work2(IA,JA)
1310 
1311  integer :: i, j
1312  integer :: ite, n
1313 
1314  ! reduce grid-scale variation
1315  do ite = 1, nite
1316  call topo_fillhalo( zsfc )
1317  work2(:,:) = zsfc(:,:)
1318  p1 => work2
1319  p2 => work1
1320  do n = 1, 4 ! 8th derivative
1321 ! do n = 1, 2 ! 4th derivative
1322  do j = js, je
1323  do i = is, ie
1324  p2(i,j) = ( - p1(i+1,j) + p1(i,j)*2.0_rp - p1(i-1,j) &
1325  - p1(i,j+1) + p1(i,j)*2.0_rp - p1(i,j-1) ) / 8.0_rp
1326  end do
1327  end do
1328  call topo_fillhalo( p2 )
1329  if ( mod(n,2) == 0 ) then
1330  p1 => work2
1331  p2 => work1
1332  else
1333  p1 => work1
1334  p2 => work2
1335  end if
1336  end do
1337  do j = js, je
1338  do i = is, ie
1339  zsfc(i,j) = max( zsfc(i,j) - p1(i,j), 0.0_rp )
1340  end do
1341  end do
1342  end do
1343 
1344  return
1345  end subroutine cnvtopo_hypdiff
1346 
1347 end module mod_cnvtopo
scale_topography::topo_write
subroutine, public topo_write
Write topography.
Definition: scale_topography.F90:164
scale_grid_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_grid_index.F90:60
scale_grid_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_grid_index.F90:63
scale_topography::topo_fillhalo
subroutine, public topo_fillhalo(Zsfc)
HALO Communication.
Definition: scale_topography.F90:98
scale_grid::dy
real(rp), public dy
length in the main region [m]: y
Definition: scale_grid_cartesian.F90:43
scale_const::const_huge
real(rp), public const_huge
huge number
Definition: scale_const.F90:38
scale_process::prc_mpistop
subroutine, public prc_mpistop
Abort MPI.
Definition: scale_process.F90:267
scale_grid::grid_fdy
real(rp), dimension(:), allocatable, public grid_fdy
y-length of grid(j+1) to grid(j) [m]
Definition: scale_grid_cartesian.F90:71
scale_grid::dx
real(rp), public dx
length in the main region [m]: x
Definition: scale_grid_cartesian.F90:42
scale_stdio::io_l
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:61
mod_cnvtopo
module Convert topography
Definition: mod_cnvtopo.f90:10
mod_cnvtopo::cnvtopo_usedem50m
logical, public cnvtopo_usedem50m
Definition: mod_cnvtopo.f90:37
scale_const::const_radius
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:46
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_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:35
mod_copytopo
module Copy topography
Definition: mod_copytopo.f90:10
scale_rm_statistics
module Statistics
Definition: scale_rm_statistics.F90:14
scale_grid_index
module grid index
Definition: scale_grid_index.F90:14
mod_cnvtopo::cnvtopo_donothing
logical, public cnvtopo_donothing
Definition: mod_cnvtopo.f90:34
scale_stdio::io_nml
logical, public io_nml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:62
scale_tracer
module TRACER
Definition: scale_tracer.F90:15
scale_grid_index::ia
integer, public ia
of whole cells: x, local, with HALO
Definition: scale_grid_index.F90:55
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_real
module GRID (real space)
Definition: scale_grid_real.F90:14
scale_grid_real::real_dlon
real(rp), dimension(:,:), allocatable, public real_dlon
delta longitude
Definition: scale_grid_real.F90:52
scale_comm
module COMMUNICATION
Definition: scale_comm.F90:23
scale_grid_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_grid_index.F90:62
mod_cnvtopo::cnvtopo_hypdiff
subroutine cnvtopo_hypdiff(Zsfc, nite)
Definition: mod_cnvtopo.f90:1301
scale_process
module PROCESS
Definition: scale_process.F90:14
scale_grid_real::real_dlat
real(rp), dimension(:,:), allocatable, public real_dlat
delta latitude
Definition: scale_grid_real.F90:53
mod_cnvtopo::cnvtopo_setup
subroutine, public cnvtopo_setup
Setup.
Definition: mod_cnvtopo.f90:71
scale_rm_statistics::stat_detail
subroutine, public stat_detail(var, varname, supress_globalcomm)
Search global maximum & minimum value.
Definition: scale_rm_statistics.F90:243
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
mod_cnvtopo::cnvtopo
subroutine, public cnvtopo
Driver.
Definition: mod_cnvtopo.f90:257
mod_cnvtopo::cnvtopo_usegtopo30
logical, public cnvtopo_usegtopo30
Definition: mod_cnvtopo.f90:35
mod_copytopo::copytopo
subroutine, public copytopo(topo_cd)
Setup and Main.
Definition: mod_copytopo.f90:84
scale_prof
module profiler
Definition: scale_prof.F90:10
scale_grid_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_grid_index.F90:61
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:36
scale_precision
module PRECISION
Definition: scale_precision.F90:13
scale_topography::topo_zsfc
real(rp), dimension(:,:), allocatable, public topo_zsfc
absolute ground height [m]
Definition: scale_topography.F90:36
scale_grid::grid_cdz
real(rp), dimension(:), allocatable, public grid_cdz
z-length of control volume [m]
Definition: scale_grid_cartesian.F90:59
scale_grid::grid_fdx
real(rp), dimension(:), allocatable, public grid_fdx
x-length of grid(i+1) to grid(i) [m]
Definition: scale_grid_cartesian.F90:70
scale_topography
module TOPOGRAPHY
Definition: scale_topography.F90:10
scale_const::const_pi
real(rp), public const_pi
pi
Definition: scale_const.F90:34
scale_stdio::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
scale_stdio::io_fid_log
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
scale_precision::rp
integer, parameter, public rp
Definition: scale_precision.F90:43
scale_grid_real::real_lonx
real(rp), dimension(:,:), allocatable, public real_lonx
longitude at staggered point (uy) [rad,0-2pi]
Definition: scale_grid_real.F90:46
scale_stdio::io_fid_nml
integer, public io_fid_nml
Log file ID (only for output namelist)
Definition: scale_stdio.F90:57
mod_cnvtopo::cnvtopo_usegmted2010
logical, public cnvtopo_usegmted2010
Definition: mod_cnvtopo.f90:36
scale_grid_real::real_laty
real(rp), dimension(:,:), allocatable, public real_laty
latitude at staggered point (xv) [rad,-pi,pi]
Definition: scale_grid_real.F90:50
scale_grid_index::ja
integer, public ja
of whole cells: y, local, with HALO
Definition: scale_grid_index.F90:56
Generated by   doxygen 1.8.14