SCALE-RM
mod_cnvtopo.F90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
10 #include "scalelib.h"
11 module mod_cnvtopo
12  !-----------------------------------------------------------------------------
13  !
14  !++ used modules
15  !
16  use scale_precision
17  use scale_io
18  use scale_prof
19  use scale_atmos_grid_cartesc_index
20  use scale_tracer
21  !-----------------------------------------------------------------------------
22  implicit none
23  private
24  !-----------------------------------------------------------------------------
25  !
26  !++ Public procedure
27  !
28  public :: cnvtopo_setup
29  public :: cnvtopo
30 
31  !-----------------------------------------------------------------------------
32  !
33  !++ Public parameters & variables
34  !
35  !-----------------------------------------------------------------------------
36  !
37  !++ Private procedure
38  !
39  private :: cnvtopo_gtopo30
40  private :: cnvtopo_gmted2010
41  private :: cnvtopo_dem50m
42  private :: cnvtopo_userfile
43  private :: cnvtopo_smooth
44 
45  !-----------------------------------------------------------------------------
46  !
47  !++ Private parameters & variables
48  !
49  character(len=H_SHORT), private :: CNVTOPO_smooth_type = 'LAPLACIAN' ! smoothing type
50  ! 'OFF' Do not apply smoothing
51  ! 'LAPLACIAN' Laplacian filter
52  ! 'GAUSSIAN' Gaussian filter
53  logical, private :: CNVTOPO_DoNothing
54  logical, private :: CNVTOPO_UseGTOPO30 = .false.
55  logical, private :: CNVTOPO_UseGMTED2010 = .false.
56  logical, private :: CNVTOPO_UseDEM50M = .false.
57  logical, private :: CNVTOPO_UseUSERFILE = .false.
58 
59  integer, private :: CNVTOPO_smooth_hypdiff_order = 4
60  integer, private :: CNVTOPO_smooth_hypdiff_niter = 20
61  logical, private :: CNVTOPO_smooth_local = .true.
62  integer, private :: CNVTOPO_smooth_itelim = 10000
63  logical, private :: CNVTOPO_smooth_trim_ocean = .true.
64  real(RP), private :: CNVTOPO_smooth_maxslope_ratio = 5.0_rp ! ratio of DZDX, DZDY
65  real(RP), private :: CNVTOPO_smooth_maxslope = -1.0_rp ! [deg]
66  real(RP), private :: CNVTOPO_smooth_maxslope_limit
67 
68  logical, private :: CNVTOPO_copy_parent = .false.
69 
70 
71  !-----------------------------------------------------------------------------
72 contains
73  !-----------------------------------------------------------------------------
75  subroutine cnvtopo_setup
76  use scale_prc, only: &
77  prc_abort
78  use scale_const, only: &
79  d2r => const_d2r, &
80  huge => const_huge
81  use scale_statistics, only: &
82  statistics_horizontal_min
83  use scale_atmos_grid_cartesc, only: &
84  cdz => atmos_grid_cartesc_cdz, &
85  fdx => atmos_grid_cartesc_fdx, &
86  fdy => atmos_grid_cartesc_fdy
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_useuserfile, &
97  cnvtopo_smooth_trim_ocean, &
98  cnvtopo_smooth_hypdiff_order, &
99  cnvtopo_smooth_hypdiff_niter, &
100  cnvtopo_smooth_maxslope_ratio, &
101  cnvtopo_smooth_maxslope, &
102  cnvtopo_smooth_local, &
103  cnvtopo_smooth_itelim, &
104  cnvtopo_smooth_type, &
105  cnvtopo_copy_parent
106 
107  real(rp) :: minslope(ia,ja)
108  real(rp) :: dxl(ia-1)
109  real(rp) :: dyl(ja-1)
110  real(rp) :: dzdx, dzdy
111 
112  integer :: ierr
113  integer :: k, i, j
114  !---------------------------------------------------------------------------
115 
116  log_newline
117  log_info("CNVTOPO_setup",*) 'Setup'
118 
119  dxl(:) = fdx(:)
120  dyl(:) = fdy(:)
121 
122  !--- read namelist
123  rewind(io_fid_conf)
124  read(io_fid_conf,nml=param_cnvtopo,iostat=ierr)
125  if( ierr < 0 ) then !--- missing
126  log_info("CNVTOPO_setup",*) 'Not found namelist. Default used.'
127  elseif( ierr > 0 ) then !--- fatal error
128  log_error("CNVTOPO_setup",*) 'Not appropriate names in namelist PARAM_CNVTOPO. Check!'
129  call prc_abort
130  endif
131  log_nml(param_cnvtopo)
132 
133  select case(cnvtopo_name)
134  case('NONE')
135  ! do nothing
136  case('GTOPO30')
137  cnvtopo_usegtopo30 = .true.
138  cnvtopo_usegmted2010 = .false.
139  cnvtopo_usedem50m = .false.
140  cnvtopo_useuserfile = .false.
141 !!$ case('GMTED2010')
142 !!$ CNVTOPO_UseGTOPO30 = .false.
143 !!$ CNVTOPO_UseGMTED2010 = .true.
144 !!$ CNVTOPO_UseDEM50M = .false.
145 !!$ CNVTOPO_UseUSERFILE = .false.
146  case('DEM50M')
147  cnvtopo_usegtopo30 = .false.
148  cnvtopo_usegmted2010 = .false.
149  cnvtopo_usedem50m = .true.
150  cnvtopo_useuserfile = .false.
151  case('COMBINE')
152  cnvtopo_usegtopo30 = .true.
153  cnvtopo_usegmted2010 = .true.
154  cnvtopo_usedem50m = .true.
155  cnvtopo_useuserfile = .false.
156  case('USERFILE')
157  ! You can use GTOPO30, GMTED2010, DEM50M and combine User-defined file as you like
158  cnvtopo_useuserfile = .true.
159  case default
160  log_error("CNVTOPO_setup",*) 'Unsupported TYPE: ', trim(cnvtopo_name)
161  call prc_abort
162  endselect
163 
164  cnvtopo_donothing = .true.
165 
166  if ( cnvtopo_usegtopo30 ) then
167  cnvtopo_donothing = .false.
168  log_info("CNVTOPO_setup",*) 'Use GTOPO, global 30 arcsec. data'
169  if ( cnvtopo_usegmted2010 ) then
170  log_info("CNVTOPO_setup",*) 'Use GMTED2010, new global 5 arcsec. data'
171  log_info("CNVTOPO_setup",*) 'Overwrite Existing region'
172  endif
173  if ( cnvtopo_usedem50m ) then
174  log_info("CNVTOPO_setup",*) 'Use DEM 50m data for Japan region'
175  log_info("CNVTOPO_setup",*) 'Overwrite Japan region'
176  endif
177  elseif ( cnvtopo_usegmted2010 ) then
178  cnvtopo_donothing = .false.
179  log_info("CNVTOPO_setup",*) 'Use GMTED2010, new global 5 arcsec. data'
180  if ( cnvtopo_usedem50m ) then
181  log_info("CNVTOPO_setup",*) 'Use DEM 50m data for Japan region'
182  log_info("CNVTOPO_setup",*) 'Overwrite Japan region'
183  endif
184  elseif ( cnvtopo_usedem50m ) then
185  cnvtopo_donothing = .false.
186  log_info("CNVTOPO_setup",*) 'Use DEM 50m data, Japan region only'
187  elseif ( cnvtopo_useuserfile ) then
188  cnvtopo_donothing = .false.
189  log_info("CNVTOPO_setup",*) 'Use user-defined file'
190  endif
191 
192  if ( cnvtopo_donothing ) then
193  log_info("CNVTOPO_setup",*) 'Do nothing for topography data'
194  endif
195 
196  if( cnvtopo_smooth_maxslope > 0.0_rp ) then
197 
198  cnvtopo_smooth_maxslope_limit = cnvtopo_smooth_maxslope
199 
200  else
201  minslope(:,:) = huge
202 
203  j = js-1
204  i = is-1
205  do k = ks, ke
206  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
207  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
208  minslope(is,js) = min( minslope(is,js), dzdx, dzdy )
209  enddo
210 
211  j = js-1
212  do i = is, ie
213  do k = ks, ke
214  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
215  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
216  minslope(i,js) = min( minslope(i,js), dzdx, dzdy )
217  enddo
218  enddo
219 
220  i = is-1
221  !$omp parallel do &
222  !$omp private(DZDX,DZDY)
223  do j = js, je
224  do k = ks, ke
225  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
226  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
227  minslope(is,j) = min( minslope(is,j), dzdx, dzdy )
228  enddo
229  enddo
230 
231  !$omp parallel do &
232  !$omp private(DZDX,DZDY)
233  do j = js, je
234  do i = is, ie
235  do k = ks, ke
236  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
237  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
238  minslope(i,j) = min( minslope(i,j), dzdx, dzdy )
239  enddo
240  enddo
241  enddo
242 
243  call statistics_horizontal_min( ia, is, ie, ja, js, je, &
244  minslope(:,:), cnvtopo_smooth_maxslope_limit )
245  end if
246 
247  return
248  end subroutine cnvtopo_setup
249 
250  !-----------------------------------------------------------------------------
252  subroutine cnvtopo
253  use scale_const, only: &
254  undef => const_undef, &
255  d2r => const_d2r
256  use scale_prc, only: &
257  prc_abort
258  use scale_topography, only: &
259  topography_fillhalo, &
260  topography_zsfc
261  use mod_copytopo, only: &
262  copytopo
263  use scale_atmos_grid_cartesc_real, only: &
264  latxv => atmos_grid_cartesc_real_latxv, &
265  lonuy => atmos_grid_cartesc_real_lonuy
266  implicit none
267 
268  integer :: i, j
269  !---------------------------------------------------------------------------
270 
271  if ( cnvtopo_donothing ) then
272  log_newline
273  log_progress(*) 'skip convert topography data'
274  else
275  log_newline
276  log_progress(*) 'start convert topography data'
277 
278  !$omp parallel do
279 !OCL XFILL
280  do j = 1, ja
281  do i = 1, ia
282  topography_zsfc(i,j) = 0.0_rp
283  end do
284  end do
285 
286  if ( cnvtopo_usegtopo30 ) then
287  call cnvtopo_gtopo30( topography_zsfc(:,:) ) ! [INOUT]
288  endif
289 
290  if ( cnvtopo_usegmted2010 ) then
291  call cnvtopo_gmted2010( topography_zsfc(:,:) ) ! [INOUT]
292  endif
293 
294  if ( cnvtopo_usedem50m ) then
295  call cnvtopo_dem50m( topography_zsfc(:,:) ) ! [INOUT]
296  endif
297 
298  if ( cnvtopo_useuserfile ) then
299  call cnvtopo_userfile( topography_zsfc(:,:) ) ! [INOUT]
300  endif
301 
302  call cnvtopo_smooth( topography_zsfc(:,:) ) ! (inout)
303  call topography_fillhalo( fill_bnd=.true. )
304 
305  if( cnvtopo_copy_parent ) call copytopo( topography_zsfc )
306 
307  log_progress(*) 'end convert topography data'
308 
309  endif
310 
311  return
312  end subroutine cnvtopo
313 
314  !-----------------------------------------------------------------------------
316  subroutine cnvtopo_gtopo30( TOPO_Zsfc )
317  use scale_prc, only: &
318  prc_abort
319  use scale_const, only: &
320  undef => const_undef
321  use mod_cnv2d, only: &
322  cnv2d_tile_init, &
323  cnv2d_exec
324  implicit none
325  real(rp), intent(inout) :: topo_zsfc(ia,ja)
326 
327  character(len=H_LONG) :: gtopo30_in_dir = '.'
328  character(len=H_LONG) :: gtopo30_in_catalogue = ''
329 
330  namelist / param_cnvtopo_gtopo30 / &
331  gtopo30_in_dir, &
332  gtopo30_in_catalogue
333 
334  ! GTOPO30 data
335  real(rp), parameter :: gtopo30_dlat = 30.0_rp / 60.0_rp / 60.0_rp ! 30 arc sec.
336  real(rp), parameter :: gtopo30_dlon = 30.0_rp / 60.0_rp / 60.0_rp ! 30 arc sec.
337 
338  ! topo
339  real(rp) :: zsfc(ia,ja)
340 
341  integer :: ierr
342  integer :: i, j
343  !---------------------------------------------------------------------------
344 
345  !--- read namelist
346  rewind(io_fid_conf)
347  read(io_fid_conf,nml=param_cnvtopo_gtopo30,iostat=ierr)
348  if( ierr < 0 ) then !--- missing
349  log_info("CNVTOPO_GTOPO30",*) 'Not found namelist. Default used.'
350  elseif( ierr > 0 ) then !--- fatal error
351  log_error("CNVTOPO_GTOPO30",*) 'Not appropriate names in namelist PARAM_CNVTOPO_GTOPO30. Check!'
352  call prc_abort
353  endif
354  log_nml(param_cnvtopo_gtopo30)
355 
356 
357  call cnv2d_tile_init( 'INT2', & ! [IN]
358  gtopo30_dlat, gtopo30_dlon, & ! [IN]
359  gtopo30_in_dir, gtopo30_in_catalogue, & ! [IN]
360  'LINEAR' ) ! [IN]
361 
362  call cnv2d_exec( zsfc(:,:), min_value = -9000.0_rp, yrevers = .true. ) ! [OUT]
363 
364  !$omp parallel do collapse(2)
365  do j = 1, ja
366  do i = 1, ia
367  if ( zsfc(i,j) /= undef ) topo_zsfc(i,j) = zsfc(i,j) ! replace data
368  end do
369  end do
370 
371  return
372  end subroutine cnvtopo_gtopo30
373 
374  !-----------------------------------------------------------------------------
376  subroutine cnvtopo_gmted2010( TOPO_Zsfc )
377  implicit none
378  real(rp), intent(inout) :: topo_zsfc(ia,ja)
379  !---------------------------------------------------------------------------
380 
381  return
382  end subroutine cnvtopo_gmted2010
383 
384  !-----------------------------------------------------------------------------
386  subroutine cnvtopo_dem50m( TOPO_Zsfc )
387  use scale_prc, only: &
388  prc_abort
389  use scale_const, only: &
390  undef => const_undef
391  use mod_cnv2d, only: &
392  cnv2d_tile_init, &
393  cnv2d_exec
394  implicit none
395  real(rp), intent(inout) :: topo_zsfc(ia,ja)
396 
397  character(len=H_LONG) :: dem50m_in_dir = '.'
398  character(len=H_LONG) :: dem50m_in_catalogue = ''
399 
400  namelist / param_cnvtopo_dem50m / &
401  dem50m_in_dir, &
402  dem50m_in_catalogue
403 
404  real(rp), parameter :: dem50m_dlat = 5.0_rp / 60.0_rp / 200.0_rp ! 30 arc sec.
405  real(rp), parameter :: dem50m_dlon = 7.5_rp / 60.0_rp / 200.0_rp ! 30 arc sec.
406 
407  ! topo
408  real(rp) :: zsfc(ia,ja)
409 
410  integer :: ierr
411  integer :: i, j
412  !---------------------------------------------------------------------------
413 
414  !--- read namelist
415  rewind(io_fid_conf)
416  read(io_fid_conf,nml=param_cnvtopo_dem50m,iostat=ierr)
417  if( ierr < 0 ) then !--- missing
418  log_info("CNVTOPO_DEM50M",*) 'Not found namelist. Default used.'
419  elseif( ierr > 0 ) then !--- fatal error
420  log_error("CNVTOPO_DEM50M",*) 'Not appropriate names in namelist PARAM_CNVTOPO_DEM50M. Check!'
421  call prc_abort
422  endif
423  log_nml(param_cnvtopo_dem50m)
424 
425  call cnv2d_tile_init( 'REAL4', & ! [IN]
426  dem50m_dlat, dem50m_dlon, & ! [IN]
427  dem50m_in_dir, dem50m_in_catalogue, & ! [IN]
428  'LINEAR' ) ! [IN]
429 
430  call cnv2d_exec( zsfc(:,:), min_value = -900.0_rp ) ! [OUT]
431 
432  !$omp parallel do collapse(2)
433  do j = 1, ja
434  do i = 1, ia
435  if ( zsfc(i,j) /= undef ) topo_zsfc(i,j) = zsfc(i,j) ! replace data
436  end do
437  end do
438 
439  return
440  end subroutine cnvtopo_dem50m
441 
442  !-----------------------------------------------------------------------------
444  subroutine cnvtopo_userfile( TOPO_Zsfc )
445  use scale_prc, only: &
446  prc_abort
447  use scale_const, only: &
448  undef => const_undef
449  use mod_cnv2d, only: &
450  cnv2d_tile_init, &
451  cnv2d_grads_init, &
452  cnv2d_exec
453  implicit none
454  real(rp), intent(inout) :: topo_zsfc(ia,ja)
455 
456  character(len=H_SHORT) :: userfile_type = '' ! "TILE" or "GrADS"
457 
458  ! TILE data
459  character(len=H_SHORT) :: userfile_dtype = 'REAL4' ! datatype (REAL4,REAL8,INT2,INT4)
460  real(rp) :: userfile_dlat = -1.0_rp ! width of latitude tile [deg.]
461  real(rp) :: userfile_dlon = -1.0_rp ! width of longitude tile [deg.]
462  character(len=H_LONG) :: userfile_dir = '.' ! directory contains data files (GrADS format)
463  character(len=H_LONG) :: userfile_catalogue = '' ! catalogue file
464  logical :: userfile_yrevers = .false. ! data of the latitude direction is stored in ordar of North->South?
465  real(rp) :: userfile_minval = 0.0_rp
466 
467  ! GrADS data
468  character(len=H_LONG) :: userfile_grads_filename = '' ! single data file (GrADS format)
469  character(len=H_SHORT) :: userfile_grads_varname = 'topo'
470  character(len=H_SHORT) :: userfile_grads_latname = 'lat'
471  character(len=H_SHORT) :: userfile_grads_lonname = 'lon'
472  character(len=H_SHORT) :: userfile_interp_type = 'LINEAR'
473  integer :: userfile_interp_level = 5
474 
475 
476  namelist / param_cnvtopo_userfile / &
477  userfile_type, &
478  userfile_dtype, &
479  userfile_dlat, &
480  userfile_dlon, &
481  userfile_catalogue, &
482  userfile_dir, &
483  userfile_yrevers, &
484  userfile_minval, &
485  userfile_grads_filename, &
486  userfile_grads_varname, &
487  userfile_grads_latname, &
488  userfile_grads_lonname, &
489  userfile_interp_type, &
490  userfile_interp_level
491 
492  ! topo
493  real(rp) :: zsfc(ia,ja)
494 
495  character(len=H_LONG) :: fname
496 
497  integer :: ierr
498  integer :: i, j
499  !---------------------------------------------------------------------------
500 
501  !--- read namelist
502 
503  !--- read namelist
504  rewind(io_fid_conf)
505  read(io_fid_conf,nml=param_cnvtopo_userfile,iostat=ierr)
506  if( ierr < 0 ) then !--- missing
507  log_info("CNVTOPO_USERFILE",*) 'Not found namelist. Default used.'
508  elseif( ierr > 0 ) then !--- fatal error
509  log_error("CNVTOPO_USERFILE",*) 'Not appropriate names in namelist PARAM_CNVTOPO_USERFILE. Check!'
510  call prc_abort
511  endif
512  log_nml(param_cnvtopo_userfile)
513 
514 
515  select case ( userfile_type )
516  case ( 'TILE' )
517 
518  if ( userfile_dlat <= 0.0_rp ) then
519  log_error("CNVTOPO_USERFILE",*) 'USERFILE_DLAT (width (deg.) of latitude tile) should be positive. Check! ', userfile_dlat
520  call prc_abort
521  endif
522  if ( userfile_dlon <= 0.0_rp ) then
523  log_error("CNVTOPO_USERFILE",*) 'USERFILE_DLON (width (deg.) of longitude tile) should be positive. Check! ', userfile_dlon
524  call prc_abort
525  endif
526  if ( userfile_catalogue == '' ) then
527  log_error("CNVTOPO_USERFILE",*) 'Catalogue file does not specified. Check!'
528  call prc_abort
529  endif
530 
531  call cnv2d_tile_init( userfile_dtype, & ! [IN]
532  userfile_dlat, userfile_dlon, & ! [IN]
533  userfile_dir, userfile_catalogue, & ! [IN]
534  'LINEAR' ) ! [IN]
535 
536  case ( 'GrADS' )
537 
538  if ( userfile_grads_filename == '' ) then
539  log_error("CNVTOPO_USERFILE",*) 'GrADS file name does not specified. Check!'
540  call prc_abort
541  endif
542 
543  call cnv2d_grads_init( userfile_grads_filename, &
544  userfile_grads_varname, &
545  userfile_grads_latname, &
546  userfile_grads_lonname, &
547  userfile_interp_type, &
548  interp_level = userfile_interp_level )
549 
550  case default
551  log_error("CNVTOPO_USERFILE",*) 'USERFILE_TYPE is invalid: ',trim(userfile_type)
552  log_error_cont(*) 'It must be "TILE" or "GrADS"'
553  call prc_abort
554  end select
555 
556  call cnv2d_exec( zsfc(:,:), & ! [OUT]
557  min_value = userfile_minval, yrevers = userfile_yrevers ) ! [IN]
558 
559  !$omp parallel do collapse(2)
560  do j = 1, ja
561  do i = 1, ia
562  if ( zsfc(i,j) /= undef ) topo_zsfc(i,j) = zsfc(i,j) ! replace data
563  end do
564  end do
565 
566  return
567  end subroutine cnvtopo_userfile
568 
569  !-----------------------------------------------------------------------------
571  subroutine cnvtopo_smooth( &
572  Zsfc )
573  use scale_const, only: &
574  eps => const_eps, &
575  d2r => const_d2r
576  use scale_prc, only: &
577  prc_abort
578  use scale_atmos_grid_cartesc, only: &
579  fdx => atmos_grid_cartesc_fdx, &
580  fdy => atmos_grid_cartesc_fdy
581  use scale_statistics, only: &
582  statistics_detail, &
583  statistics_horizontal_max
584  use scale_topography, only: &
585  topography_fillhalo
586  use scale_filter, only: &
587  filter_hyperdiff
588  use scale_landuse, only: &
589  landuse_fact_ocean
590  implicit none
591 
592  real(rp), intent(inout) :: zsfc(ia,ja)
593 
594  real(rp) :: dzsfc_dxy(ia,ja,2) ! d(Zsfc)/dx at u-position and d(Zsfc)/dy at v-position
595 
596  real(rp) :: dxl(ia-1)
597  real(rp) :: dyl(ja-1)
598 
599  real(rp) :: flx_x(ia,ja)
600  real(rp) :: flx_y(ia,ja)
601 
602  real(rp) :: slope(ia,ja)
603  real(rp) :: maxslope
604  real(rp), pointer :: topo_sign(:,:)
605  real(rp) :: flag, ocean_flag
606 
607  character(len=8), parameter :: varname(2) = (/ "DZsfc_DX", "DZsfc_DY" /)
608 
609 
610  integer :: ite
611  integer :: i, j
612  !---------------------------------------------------------------------------
613 
614  if ( cnvtopo_smooth_type == 'OFF' ) then
615  log_newline
616  log_info("CNVTOPO_smooth",*) 'Do not apply smoothing.'
617 
618  return
619  else
620  log_newline
621  log_info("CNVTOPO_smooth",*) 'Apply smoothing.'
622  log_info_cont(*) 'Slope limit = ', cnvtopo_smooth_maxslope_limit
623  log_info_cont(*) 'Smoothing type = ', cnvtopo_smooth_type
624  log_info_cont(*) 'Smoothing locally = ', cnvtopo_smooth_local
625  log_newline
626  endif
627 
628  dxl(:) = fdx(:)
629  dyl(:) = fdy(:)
630 
631  if ( cnvtopo_smooth_trim_ocean ) then
632  allocate( topo_sign(ia,ja) )
633  !$omp parallel do &
634  !$omp private(ocean_flag)
635  do j = 1, ja
636  do i = 1, ia
637  ocean_flag = ( 0.5_rp + sign( 0.5_rp, landuse_fact_ocean(i,j) - 1.0_rp + eps ) ) & ! fact_ocean==1
638  * ( 0.5_rp + sign( 0.5_rp, eps - abs(zsfc(i,j)) ) ) ! |Zsfc| < EPS
639  topo_sign(i,j) = sign( 1.0_rp, zsfc(i,j) ) * ( 1.0_rp - ocean_flag )
640  end do
641  end do
642  else
643  topo_sign => null()
644  end if
645 
646  ! digital filter
647  do ite = 1, cnvtopo_smooth_itelim+1
648  log_progress(*) 'smoothing itelation : ', ite
649 
650  call topography_fillhalo( zsfc=zsfc(:,:), fill_bnd=.true. )
651 
652  !$omp parallel do
653  do j = 1, ja
654  do i = 1, ia-1
655  dzsfc_dxy(i,j,1) = atan2( ( zsfc(i+1,j)-zsfc(i,j) ), dxl(i) ) / d2r
656  enddo
657  enddo
658  dzsfc_dxy(ia,:,1) = 0.0_rp
659  !$omp parallel do
660  do j = 1, ja-1
661  do i = 1, ia
662  dzsfc_dxy(i,j,2) = atan2( ( zsfc(i,j+1)-zsfc(i,j) ), dyl(j) ) / d2r
663  enddo
664  enddo
665  dzsfc_dxy(:,ja,2) = 0.0_rp
666 
667  slope(:,:) = max( abs(dzsfc_dxy(:,:,1)), abs(dzsfc_dxy(:,:,2)) )
668  call statistics_horizontal_max( ia, is, ie, ja, js, je, &
669  slope(:,:), maxslope )
670 
671  log_progress(*) 'maximum slope [deg] : ', maxslope
672 
673  if( maxslope < cnvtopo_smooth_maxslope_limit ) exit
674 
675  call statistics_detail( ia, is, ie, ja, js, je, 2, &
676  varname(:), dzsfc_dxy(:,:,:) )
677 
678  select case( cnvtopo_smooth_type )
679  case( 'GAUSSIAN' )
680 
681  ! 3 by 3 gaussian filter
682  !$omp parallel do
683  do j = js, je
684  do i = is, ie
685  zsfc(i,j) = ( 0.2500_rp * zsfc(i ,j ) &
686  + 0.1250_rp * zsfc(i-1,j ) &
687  + 0.1250_rp * zsfc(i+1,j ) &
688  + 0.1250_rp * zsfc(i ,j-1) &
689  + 0.1250_rp * zsfc(i ,j+1) &
690  + 0.0625_rp * zsfc(i-1,j-1) &
691  + 0.0625_rp * zsfc(i+1,j-1) &
692  + 0.0625_rp * zsfc(i-1,j+1) &
693  + 0.0625_rp * zsfc(i+1,j+1) )
694  enddo
695  enddo
696 
697  case( 'LAPLACIAN' )
698 
699  !$omp parallel do
700  do j = js , je
701  do i = is-1, ie
702  flx_x(i,j) = zsfc(i+1,j) - zsfc(i,j)
703 ! FLX_TMP(i,j) = Zsfc(i+1,j) - Zsfc(i,j)
704  enddo
705  enddo
706 !!$ call TOPOGRAPHY_fillhalo( FLX_TMP )
707 !!$ do j = JS , JE
708 !!$ do i = IS-1, IE
709 !!$ FLX_X(i,j) = - ( FLX_TMP(i+1,j) - FLX_TMP(i,j) )
710 !!$ enddo
711 !!$ enddo
712 
713  !$omp parallel do
714  do j = js-1, je
715  do i = is , ie
716  flx_y(i,j) = zsfc(i,j+1) - zsfc(i,j)
717 ! FLX_TMP(i,j) = Zsfc(i,j+1) - Zsfc(i,j)
718  enddo
719  enddo
720 !!$ call TOPOGRAPHY_fillhalo( FLX_TMP )
721 !!$ do j = JS-1, JE
722 !!$ do i = IS , IE
723 !!$ FLX_Y(i,j) = - ( FLX_TMP(i,j+1) - FLX_TMP(i,j) )
724 !!$ enddo
725 !!$ enddo
726 
727 
728  if ( cnvtopo_smooth_local ) then
729  !$omp parallel do &
730  !$omp private(flag)
731  do j = js , je
732  do i = is-1, ie
733  flag = 0.5_rp &
734  + sign(0.5_rp, max( abs(dzsfc_dxy(i+1,j ,1)), &
735  abs(dzsfc_dxy(i ,j ,1)), &
736  abs(dzsfc_dxy(i-1,j ,1)), &
737  abs(dzsfc_dxy(i+1,j ,2)), &
738  abs(dzsfc_dxy(i+1,j-1,2)), &
739  abs(dzsfc_dxy(i ,j ,2)), &
740  abs(dzsfc_dxy(i ,j-1,2)) &
741  ) - cnvtopo_smooth_maxslope_limit )
742  flx_x(i,j) = flx_x(i,j) * flag
743  enddo
744  enddo
745  !$omp parallel do &
746  !$omp private(flag)
747  do j = js-1, je
748  do i = is , ie
749  flag = 0.5_rp &
750  + sign(0.5_rp, max( abs(dzsfc_dxy(i ,j+1,2)), &
751  abs(dzsfc_dxy(i ,j ,2)), &
752  abs(dzsfc_dxy(i ,j-1,2)), &
753  abs(dzsfc_dxy(i ,j+1,1)), &
754  abs(dzsfc_dxy(i-1,j+1,1)), &
755  abs(dzsfc_dxy(i ,j ,1)), &
756  abs(dzsfc_dxy(i-1,j ,1)) &
757  ) - cnvtopo_smooth_maxslope_limit )
758  flx_y(i,j) = flx_y(i,j) * flag
759  enddo
760  enddo
761  endif
762 
763  !$omp parallel do
764  do j = js, je
765  do i = is, ie
766  zsfc(i,j) = zsfc(i,j) &
767  + 0.1_rp * ( ( flx_x(i,j) - flx_x(i-1,j) ) &
768  + ( flx_y(i,j) - flx_y(i,j-1) ) )
769  enddo
770  enddo
771 
772  case default
773  log_error("CNVTOPO_smooth",*) 'Invalid smoothing type'
774  call prc_abort
775  end select
776 
777  if ( cnvtopo_smooth_trim_ocean ) then
778  !$omp parallel do
779  do j = js, je
780  do i = is, ie
781  zsfc(i,j) = sign( max( zsfc(i,j) * topo_sign(i,j), 0.0_rp ), topo_sign(i,j) )
782  end do
783  end do
784  end if
785 
786  enddo
787 
788  if ( ite > cnvtopo_smooth_itelim ) then
789  log_error("CNVTOPO_smooth",*) 'Smoothing did not converge until ', cnvtopo_smooth_itelim,' times of iteration.'
790 
791  log_error_cont(*) 'Please try different parameters of PARAM_CNVTOPO.'
792  log_error_cont(*) '- Number limit of iteration (CNVTOPO_smooth_itelim) = ', cnvtopo_smooth_itelim
793  log_error_cont(*) '- Maximum ratio of slope dZ/dX, dZ/dY (CNVTOPO_smooth_maxslope_ratio) = ', cnvtopo_smooth_maxslope_ratio
794  log_error_cont(*) ' Or, Maximum of slope with degree (CNVTOPO_smooth_maxslope) = ', cnvtopo_smooth_maxslope
795  log_error_cont(*) '- Smoothing type LAPLACIAN/GAUSSIAN/OFF (CNVTOPO_smooth_type) = ', trim(cnvtopo_smooth_type)
796  call prc_abort
797  else
798  log_newline
799  log_info("CNVTOPO_smooth",*) 'smoothing complete.'
800  endif
801 
802 
803 
804  if ( cnvtopo_smooth_hypdiff_niter > 0 ) then
805 
806  log_newline
807  log_info("CNVTOPO_smooth",*) 'Apply hyperdiffusion.'
808 
809  call filter_hyperdiff( ia, is, ie, ja, js, je, &
810  zsfc(:,:), & ! (inout)
811  cnvtopo_smooth_hypdiff_order, cnvtopo_smooth_hypdiff_niter, & ! (in)
812  limiter_sign = topo_sign(:,:) ) ! (in)
813 
814  call topography_fillhalo( zsfc=zsfc(:,:), fill_bnd=.true. )
815 
816  !$omp parallel do
817  do j = 1, ja
818  do i = 1, ia-1
819  dzsfc_dxy(i,j,1) = atan2( ( zsfc(i+1,j)-zsfc(i,j) ), dxl(i) ) / d2r
820  enddo
821  enddo
822  dzsfc_dxy(ia,:,1) = 0.0_rp
823  !$omp parallel do
824  do j = 1, ja-1
825  do i = 1, ia
826  dzsfc_dxy(i,j,2) = atan2( ( zsfc(i,j+1)-zsfc(i,j) ), dyl(j) ) / d2r
827  enddo
828  enddo
829  dzsfc_dxy(:,ja,2) = 0.0_rp
830 
831  slope(:,:) = max( abs(dzsfc_dxy(:,:,1)), abs(dzsfc_dxy(:,:,2)) )
832  call statistics_horizontal_max( ia, is, ie, ja, js, je, &
833  slope(:,:), maxslope )
834 
835  log_info("CNVTOPO_smooth",*) 'maximum slope [deg] : ', maxslope
836 
837  end if
838 
839  call topography_fillhalo( zsfc=zsfc(:,:), fill_bnd=.true. )
840 
841  call statistics_detail( ia, is, ie, ja, js, je, 2, &
842  varname(:), dzsfc_dxy(:,:,:) )
843 
844  log_newline
845 
846  return
847  end subroutine cnvtopo_smooth
848 
849 end module mod_cnvtopo
mod_cnv2d::cnv2d_grads_init
subroutine, public cnv2d_grads_init(FILE_NAME, VAR_NAME, LAT_NAME, LON_NAME, interp_type, interp_level, search_limit, POSTFIX)
Definition: mod_cnv2d.F90:227
scale_statistics
module Statistics
Definition: scale_statistics.F90:11
scale_landuse::landuse_fact_ocean
real(rp), dimension(:,:), allocatable, public landuse_fact_ocean
ocean factor
Definition: scale_landuse.F90:44
scale_atmos_grid_cartesc_index::ke
integer, public ke
end point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:52
mod_cnv2d
module Convert 2D data
Definition: mod_cnv2d.F90:11
scale_prc::prc_abort
subroutine, public prc_abort
Abort Process.
Definition: scale_prc.F90:342
mod_cnvtopo::cnvtopo
subroutine, public cnvtopo
Driver.
Definition: mod_cnvtopo.F90:253
scale_atmos_grid_cartesc::atmos_grid_cartesc_cdz
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_cdz
z-length of control volume [m]
Definition: scale_atmos_grid_cartesC.F90:42
scale_precision
module PRECISION
Definition: scale_precision.F90:14
mod_cnvtopo
module Convert topography
Definition: mod_cnvtopo.F90:11
scale_const::const_eps
real(rp), public const_eps
small number
Definition: scale_const.F90:33
scale_atmos_grid_cartesc::atmos_grid_cartesc_fdy
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fdy
y-length of grid(j+1) to grid(j) [m]
Definition: scale_atmos_grid_cartesC.F90:63
mod_copytopo::copytopo
subroutine, public copytopo(TOPO_child)
Setup and Main.
Definition: mod_copytopo.F90:74
scale_topography
module TOPOGRAPHY
Definition: scale_topography.F90:11
scale_atmos_grid_cartesc_real
module Atmosphere GRID CartesC Real(real space)
Definition: scale_atmos_grid_cartesC_real.F90:11
scale_prc
module PROCESS
Definition: scale_prc.F90:11
scale_atmos_grid_cartesc_index::ie
integer, public ie
end point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:54
scale_io
module STDIO
Definition: scale_io.F90:10
scale_tracer::k
real(rp), public k
Definition: scale_tracer.F90:44
scale_atmos_grid_cartesc_index
module atmosphere / grid / cartesC index
Definition: scale_atmos_grid_cartesC_index.F90:12
scale_const
module CONSTANT
Definition: scale_const.F90:11
scale_atmos_grid_cartesc_index::ia
integer, public ia
Definition: scale_atmos_grid_cartesC_index.F90:48
scale_const::const_huge
real(rp), public const_huge
huge number
Definition: scale_const.F90:35
scale_prof
module profiler
Definition: scale_prof.F90:11
scale_atmos_grid_cartesc_index::is
integer, public is
start point of inner domain: x, local
Definition: scale_atmos_grid_cartesC_index.F90:53
scale_atmos_grid_cartesc_index::ja
integer, public ja
Definition: scale_atmos_grid_cartesC_index.F90:49
mod_cnv2d::cnv2d_exec
subroutine, public cnv2d_exec(var, step, min_value, yrevers)
Definition: mod_cnv2d.F90:330
scale_atmos_grid_cartesc::atmos_grid_cartesc_fdx
real(rp), dimension(:), allocatable, public atmos_grid_cartesc_fdx
x-length of grid(i+1) to grid(i) [m]
Definition: scale_atmos_grid_cartesC.F90:62
mod_copytopo
module Copy topography
Definition: mod_copytopo.F90:11
scale_tracer
module TRACER
Definition: scale_tracer.F90:12
scale_atmos_grid_cartesc_real::atmos_grid_cartesc_real_latxv
real(rp), dimension(:,:), allocatable, public atmos_grid_cartesc_real_latxv
latitude at staggered point (xv) [rad,-pi,pi]
Definition: scale_atmos_grid_cartesC_real.F90:54
scale_atmos_grid_cartesc_index::ks
integer, public ks
start point of inner domain: z, local
Definition: scale_atmos_grid_cartesC_index.F90:51
scale_topography::topography_zsfc
real(rp), dimension(:,:), allocatable, public topography_zsfc
absolute ground height [m]
Definition: scale_topography.F90:38
mod_cnv2d::cnv2d_tile_init
subroutine, public cnv2d_tile_init(dtype, dlat, dlon, dir, catalogue, interp_type, interp_level, nmax)
Definition: mod_cnv2d.F90:126
scale_atmos_grid_cartesc_index::js
integer, public js
start point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:55
scale_filter
module FILTER
Definition: scale_filter.F90:11
scale_const::const_d2r
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:32
scale_atmos_grid_cartesc
module atmosphere / grid / cartesC
Definition: scale_atmos_grid_cartesC.F90:12
scale_landuse
module LANDUSE
Definition: scale_landuse.F90:19
mod_cnvtopo::cnvtopo_setup
subroutine, public cnvtopo_setup
Setup.
Definition: mod_cnvtopo.F90:76
scale_const::const_undef
real(rp), public const_undef
Definition: scale_const.F90:41
scale_topography::topography_fillhalo
subroutine, public topography_fillhalo(Zsfc, FILL_BND)
HALO Communication.
Definition: scale_topography.F90:117
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:56
scale_atmos_grid_cartesc_index::je
integer, public je
end point of inner domain: y, local
Definition: scale_atmos_grid_cartesC_index.F90:56
scale_atmos_grid_cartesc_real::atmos_grid_cartesc_real_lonuy
real(rp), dimension(:,:), allocatable, public atmos_grid_cartesc_real_lonuy
longitude at staggered point (uy) [rad,0-2pi]
Definition: scale_atmos_grid_cartesC_real.F90:49
Generated by   doxygen 1.8.17