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 
202  !$acc data create(minslope)
203 
204  !$acc kernels
205  minslope(:,:) = huge
206  !$acc end kernels
207 
208  j = js-1
209  i = is-1
210  !$acc kernels
211  !$acc loop independent
212  do k = ks, ke
213  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
214  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
215  minslope(is,js) = min( minslope(is,js), dzdx, dzdy )
216  enddo
217  !$accend kernels
218 
219  j = js-1
220  !$acc kernels
221  !$acc loop collapse(2) independent
222  do i = is, ie
223  do k = ks, ke
224  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
225  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
226  minslope(i,js) = min( minslope(i,js), dzdx, dzdy )
227  enddo
228  enddo
229  !$acc end kernels
230 
231  i = is-1
232  !$omp parallel do &
233  !$omp private(DZDX,DZDY)
234  !$acc kernels
235  !$acc loop collapse(2) independent
236  do j = js, je
237  do k = ks, ke
238  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
239  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
240  minslope(is,j) = min( minslope(is,j), dzdx, dzdy )
241  enddo
242  enddo
243  !$acc end kernels
244 
245  !$omp parallel do &
246  !$omp private(DZDX,DZDY)
247  !$acc kernels
248  !$acc loop collapse(3) independent
249  do j = js, je
250  do i = is, ie
251  do k = ks, ke
252  dzdx = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dxl(i) ) / d2r
253  dzdy = atan2( cnvtopo_smooth_maxslope_ratio * cdz(k), dyl(j) ) / d2r
254  minslope(i,j) = min( minslope(i,j), dzdx, dzdy )
255  enddo
256  enddo
257  enddo
258  !$acc end kernels
259 
260  call statistics_horizontal_min( ia, is, ie, ja, js, je, &
261  minslope(:,:), cnvtopo_smooth_maxslope_limit )
262 
263  !$acc end data
264 
265  end if
266 
267  return
268  end subroutine cnvtopo_setup
269 
270  !-----------------------------------------------------------------------------
272  subroutine cnvtopo
273  use scale_const, only: &
274  undef => const_undef, &
275  d2r => const_d2r
276  use scale_prc, only: &
277  prc_abort
278  use scale_topography, only: &
279  topography_fillhalo, &
280  topography_zsfc
281  use mod_copytopo, only: &
282  copytopo
283  use scale_atmos_grid_cartesc_real, only: &
284  latxv => atmos_grid_cartesc_real_latxv, &
285  lonuy => atmos_grid_cartesc_real_lonuy
286  implicit none
287 
288  integer :: i, j
289  !---------------------------------------------------------------------------
290 
291  if ( cnvtopo_donothing ) then
292  log_newline
293  log_progress(*) 'skip convert topography data'
294  else
295  log_newline
296  log_progress(*) 'start convert topography data'
297 
298  !$omp parallel do
299 !OCL XFILL
300  !$acc kernels
301  do j = 1, ja
302  do i = 1, ia
303  topography_zsfc(i,j) = 0.0_rp
304  end do
305  end do
306  !$acc end kernels
307 
308  if ( cnvtopo_usegtopo30 ) then
309  call cnvtopo_gtopo30( topography_zsfc(:,:) ) ! [INOUT]
310  endif
311 
312  if ( cnvtopo_usegmted2010 ) then
313  call cnvtopo_gmted2010( topography_zsfc(:,:) ) ! [INOUT]
314  endif
315 
316  if ( cnvtopo_usedem50m ) then
317  call cnvtopo_dem50m( topography_zsfc(:,:) ) ! [INOUT]
318  endif
319 
320  if ( cnvtopo_useuserfile ) then
321  call cnvtopo_userfile( topography_zsfc(:,:) ) ! [INOUT]
322  endif
323 
324  call cnvtopo_smooth( topography_zsfc(:,:) ) ! (inout)
325  call topography_fillhalo( fill_bnd=.true. )
326 
327  if( cnvtopo_copy_parent ) call copytopo( topography_zsfc )
328 
329  log_progress(*) 'end convert topography data'
330 
331  endif
332 
333  return
334  end subroutine cnvtopo
335 
336  !-----------------------------------------------------------------------------
338  subroutine cnvtopo_gtopo30( TOPO_Zsfc )
339  use scale_prc, only: &
340  prc_abort
341  use scale_const, only: &
342  undef => const_undef
343  use mod_cnv2d, only: &
344  cnv2d_tile_init, &
345  cnv2d_exec
346  implicit none
347  real(rp), intent(inout) :: topo_zsfc(ia,ja)
348 
349  character(len=H_LONG) :: gtopo30_in_dir = '.'
350  character(len=H_LONG) :: gtopo30_in_catalogue = ''
351 
352  namelist / param_cnvtopo_gtopo30 / &
353  gtopo30_in_dir, &
354  gtopo30_in_catalogue
355 
356  ! GTOPO30 data
357  real(rp), parameter :: gtopo30_dlat = 30.0_rp / 60.0_rp / 60.0_rp ! 30 arc sec.
358  real(rp), parameter :: gtopo30_dlon = 30.0_rp / 60.0_rp / 60.0_rp ! 30 arc sec.
359 
360  ! topo
361  real(rp) :: zsfc(ia,ja)
362 
363  integer :: ierr
364  integer :: i, j
365  !---------------------------------------------------------------------------
366 
367  !--- read namelist
368  rewind(io_fid_conf)
369  read(io_fid_conf,nml=param_cnvtopo_gtopo30,iostat=ierr)
370  if( ierr < 0 ) then !--- missing
371  log_info("CNVTOPO_GTOPO30",*) 'Not found namelist. Default used.'
372  elseif( ierr > 0 ) then !--- fatal error
373  log_error("CNVTOPO_GTOPO30",*) 'Not appropriate names in namelist PARAM_CNVTOPO_GTOPO30. Check!'
374  call prc_abort
375  endif
376  log_nml(param_cnvtopo_gtopo30)
377 
378 
379  call cnv2d_tile_init( 'INT2', & ! [IN]
380  gtopo30_dlat, gtopo30_dlon, & ! [IN]
381  gtopo30_in_dir, gtopo30_in_catalogue, & ! [IN]
382  'LINEAR' ) ! [IN]
383 
384  call cnv2d_exec( zsfc(:,:), min_value = -9000.0_rp, yrevers = .true. ) ! [OUT]
385 
386  !$omp parallel do collapse(2)
387  !$acc kernels
388  do j = 1, ja
389  do i = 1, ia
390  if ( zsfc(i,j) /= undef ) topo_zsfc(i,j) = zsfc(i,j) ! replace data
391  end do
392  end do
393  !$acc end kernels
394 
395  return
396  end subroutine cnvtopo_gtopo30
397 
398  !-----------------------------------------------------------------------------
400  subroutine cnvtopo_gmted2010( TOPO_Zsfc )
401  implicit none
402  real(rp), intent(inout) :: topo_zsfc(ia,ja)
403  !---------------------------------------------------------------------------
404 
405  return
406  end subroutine cnvtopo_gmted2010
407 
408  !-----------------------------------------------------------------------------
410  subroutine cnvtopo_dem50m( TOPO_Zsfc )
411  use scale_prc, only: &
412  prc_abort
413  use scale_const, only: &
414  undef => const_undef
415  use mod_cnv2d, only: &
416  cnv2d_tile_init, &
417  cnv2d_exec
418  implicit none
419  real(rp), intent(inout) :: topo_zsfc(ia,ja)
420 
421  character(len=H_LONG) :: dem50m_in_dir = '.'
422  character(len=H_LONG) :: dem50m_in_catalogue = ''
423 
424  namelist / param_cnvtopo_dem50m / &
425  dem50m_in_dir, &
426  dem50m_in_catalogue
427 
428  real(rp), parameter :: dem50m_dlat = 5.0_rp / 60.0_rp / 200.0_rp ! 30 arc sec.
429  real(rp), parameter :: dem50m_dlon = 7.5_rp / 60.0_rp / 200.0_rp ! 30 arc sec.
430 
431  ! topo
432  real(rp) :: zsfc(ia,ja)
433 
434  integer :: ierr
435  integer :: i, j
436  !---------------------------------------------------------------------------
437 
438  !--- read namelist
439  rewind(io_fid_conf)
440  read(io_fid_conf,nml=param_cnvtopo_dem50m,iostat=ierr)
441  if( ierr < 0 ) then !--- missing
442  log_info("CNVTOPO_DEM50M",*) 'Not found namelist. Default used.'
443  elseif( ierr > 0 ) then !--- fatal error
444  log_error("CNVTOPO_DEM50M",*) 'Not appropriate names in namelist PARAM_CNVTOPO_DEM50M. Check!'
445  call prc_abort
446  endif
447  log_nml(param_cnvtopo_dem50m)
448 
449  call cnv2d_tile_init( 'REAL4', & ! [IN]
450  dem50m_dlat, dem50m_dlon, & ! [IN]
451  dem50m_in_dir, dem50m_in_catalogue, & ! [IN]
452  'LINEAR' ) ! [IN]
453 
454  call cnv2d_exec( zsfc(:,:), min_value = -900.0_rp ) ! [OUT]
455 
456  !$omp parallel do collapse(2)
457  !$acc kernels
458  do j = 1, ja
459  do i = 1, ia
460  if ( zsfc(i,j) /= undef ) topo_zsfc(i,j) = zsfc(i,j) ! replace data
461  end do
462  end do
463  !$acc end kernels
464 
465  return
466  end subroutine cnvtopo_dem50m
467 
468  !-----------------------------------------------------------------------------
470  subroutine cnvtopo_userfile( TOPO_Zsfc )
471  use scale_prc, only: &
472  prc_abort
473  use scale_const, only: &
474  undef => const_undef
475  use mod_cnv2d, only: &
476  cnv2d_tile_init, &
477  cnv2d_grads_init, &
478  cnv2d_exec
479  implicit none
480  real(rp), intent(inout) :: topo_zsfc(ia,ja)
481 
482  character(len=H_SHORT) :: userfile_type = '' ! "TILE" or "GrADS"
483 
484  ! TILE data
485  character(len=H_SHORT) :: userfile_dtype = 'REAL4' ! datatype (REAL4,REAL8,INT2,INT4)
486  real(rp) :: userfile_dlat = -1.0_rp ! width of latitude tile [deg.]
487  real(rp) :: userfile_dlon = -1.0_rp ! width of longitude tile [deg.]
488  character(len=H_LONG) :: userfile_dir = '.' ! directory contains data files (GrADS format)
489  character(len=H_LONG) :: userfile_catalogue = '' ! catalogue file
490  logical :: userfile_yrevers = .false. ! data of the latitude direction is stored in ordar of North->South?
491  real(rp) :: userfile_minval = 0.0_rp
492 
493  ! GrADS data
494  character(len=H_LONG) :: userfile_grads_filename = '' ! single data file (GrADS format)
495  character(len=H_SHORT) :: userfile_grads_varname = 'topo'
496  character(len=H_SHORT) :: userfile_grads_latname = 'lat'
497  character(len=H_SHORT) :: userfile_grads_lonname = 'lon'
498  character(len=H_SHORT) :: userfile_interp_type = 'LINEAR'
499  integer :: userfile_interp_level = 5
500 
501 
502  namelist / param_cnvtopo_userfile / &
503  userfile_type, &
504  userfile_dtype, &
505  userfile_dlat, &
506  userfile_dlon, &
507  userfile_catalogue, &
508  userfile_dir, &
509  userfile_yrevers, &
510  userfile_minval, &
511  userfile_grads_filename, &
512  userfile_grads_varname, &
513  userfile_grads_latname, &
514  userfile_grads_lonname, &
515  userfile_interp_type, &
516  userfile_interp_level
517 
518  ! topo
519  real(rp) :: zsfc(ia,ja)
520 
521  character(len=H_LONG) :: fname
522 
523  integer :: ierr
524  integer :: i, j
525  !---------------------------------------------------------------------------
526 
527  !--- read namelist
528 
529  !--- read namelist
530  rewind(io_fid_conf)
531  read(io_fid_conf,nml=param_cnvtopo_userfile,iostat=ierr)
532  if( ierr < 0 ) then !--- missing
533  log_info("CNVTOPO_USERFILE",*) 'Not found namelist. Default used.'
534  elseif( ierr > 0 ) then !--- fatal error
535  log_error("CNVTOPO_USERFILE",*) 'Not appropriate names in namelist PARAM_CNVTOPO_USERFILE. Check!'
536  call prc_abort
537  endif
538  log_nml(param_cnvtopo_userfile)
539 
540 
541  select case ( userfile_type )
542  case ( 'TILE' )
543 
544  if ( userfile_dlat <= 0.0_rp ) then
545  log_error("CNVTOPO_USERFILE",*) 'USERFILE_DLAT (width (deg.) of latitude tile) should be positive. Check! ', userfile_dlat
546  call prc_abort
547  endif
548  if ( userfile_dlon <= 0.0_rp ) then
549  log_error("CNVTOPO_USERFILE",*) 'USERFILE_DLON (width (deg.) of longitude tile) should be positive. Check! ', userfile_dlon
550  call prc_abort
551  endif
552  if ( userfile_catalogue == '' ) then
553  log_error("CNVTOPO_USERFILE",*) 'Catalogue file does not specified. Check!'
554  call prc_abort
555  endif
556 
557  call cnv2d_tile_init( userfile_dtype, & ! [IN]
558  userfile_dlat, userfile_dlon, & ! [IN]
559  userfile_dir, userfile_catalogue, & ! [IN]
560  'LINEAR' ) ! [IN]
561 
562  case ( 'GrADS' )
563 
564  if ( userfile_grads_filename == '' ) then
565  log_error("CNVTOPO_USERFILE",*) 'GrADS file name does not specified. Check!'
566  call prc_abort
567  endif
568 
569  call cnv2d_grads_init( userfile_grads_filename, &
570  userfile_grads_varname, &
571  userfile_grads_latname, &
572  userfile_grads_lonname, &
573  userfile_interp_type, &
574  userfile_interp_level )
575 
576  case default
577  log_error("CNVTOPO_USERFILE",*) 'USERFILE_TYPE is invalid: ',trim(userfile_type)
578  log_error_cont(*) 'It must be "TILE" or "GrADS"'
579  call prc_abort
580  end select
581 
582  call cnv2d_exec( zsfc(:,:), & ! [OUT]
583  min_value = userfile_minval, yrevers = userfile_yrevers ) ! [IN]
584 
585  !$omp parallel do collapse(2)
586  !$acc kernels
587  do j = 1, ja
588  do i = 1, ia
589  if ( zsfc(i,j) /= undef ) topo_zsfc(i,j) = zsfc(i,j) ! replace data
590  end do
591  end do
592  !$acc end kernels
593 
594  return
595  end subroutine cnvtopo_userfile
596 
597  !-----------------------------------------------------------------------------
599  subroutine cnvtopo_smooth( &
600  Zsfc )
601  use scale_const, only: &
602  eps => const_eps, &
603  d2r => const_d2r
604  use scale_prc, only: &
605  prc_abort
606  use scale_atmos_grid_cartesc, only: &
607  fdx => atmos_grid_cartesc_fdx, &
608  fdy => atmos_grid_cartesc_fdy
609  use scale_statistics, only: &
610  statistics_detail, &
611  statistics_horizontal_max
612  use scale_topography, only: &
613  topography_fillhalo
614  use scale_filter, only: &
615  filter_hyperdiff
616  use scale_landuse, only: &
617  landuse_fact_ocean
618  implicit none
619 
620  real(rp), intent(inout) :: zsfc(ia,ja)
621 
622  real(rp) :: dzsfc_dxy(ia,ja,2) ! d(Zsfc)/dx at u-position and d(Zsfc)/dy at v-position
623 
624  real(rp) :: dxl(ia-1)
625  real(rp) :: dyl(ja-1)
626 
627  real(rp) :: flx_x(ia,ja)
628  real(rp) :: flx_y(ia,ja)
629 
630  real(rp) :: slope(ia,ja)
631  real(rp) :: maxslope
632  real(rp), pointer :: topo_sign(:,:)
633  real(rp), allocatable, target :: topo_sign_t(:,:)
634  real(rp) :: flag, ocean_flag
635 
636  character(len=8), parameter :: varname(2) = (/ "DZsfc_DX", "DZsfc_DY" /)
637 
638 
639  integer :: ite
640  integer :: i, j
641  !---------------------------------------------------------------------------
642 
643  if ( cnvtopo_smooth_type == 'OFF' ) then
644  log_newline
645  log_info("CNVTOPO_smooth",*) 'Do not apply smoothing.'
646 
647  return
648  else
649  log_newline
650  log_info("CNVTOPO_smooth",*) 'Apply smoothing.'
651  log_info_cont(*) 'Slope limit = ', cnvtopo_smooth_maxslope_limit
652  log_info_cont(*) 'Smoothing type = ', cnvtopo_smooth_type
653  log_info_cont(*) 'Smoothing locally = ', cnvtopo_smooth_local
654  log_newline
655  endif
656 
657  dxl(:) = fdx(:)
658  dyl(:) = fdy(:)
659 
660  if ( cnvtopo_smooth_trim_ocean ) then
661  allocate( topo_sign_t(ia,ja) )
662  topo_sign => topo_sign_t
663  !$acc enter data create(TOPO_sign)
664  !$omp parallel do &
665  !$omp private(ocean_flag)
666  !$acc kernels
667  do j = 1, ja
668  do i = 1, ia
669  ocean_flag = ( 0.5_rp + sign( 0.5_rp, landuse_fact_ocean(i,j) - 1.0_rp + eps ) ) & ! fact_ocean==1
670  * ( 0.5_rp + sign( 0.5_rp, eps - abs(zsfc(i,j)) ) ) ! |Zsfc| < EPS
671  topo_sign(i,j) = sign( 1.0_rp, zsfc(i,j) ) * ( 1.0_rp - ocean_flag )
672  end do
673  end do
674  !$acc end kernels
675  else
676  topo_sign => null()
677  end if
678 
679  !$acc data copyin(DXL, DYL) create(DZsfc_DXY, FLX_X, FLX_Y)
680 
681  ! digital filter
682  do ite = 1, cnvtopo_smooth_itelim+1
683  log_progress(*) 'smoothing itelation : ', ite
684 
685  call topography_fillhalo( zsfc=zsfc(:,:), fill_bnd=.true. )
686 
687  !$omp parallel do
688  !$acc kernels
689  do j = 1, ja
690  do i = 1, ia-1
691  dzsfc_dxy(i,j,1) = atan2( ( zsfc(i+1,j)-zsfc(i,j) ), dxl(i) ) / d2r
692  enddo
693  enddo
694  !$acc end kernels
695  !$acc kernels
696  dzsfc_dxy(ia,:,1) = 0.0_rp
697  !$acc end kernels
698  !$omp parallel do
699  !$acc kernels
700  do j = 1, ja-1
701  do i = 1, ia
702  dzsfc_dxy(i,j,2) = atan2( ( zsfc(i,j+1)-zsfc(i,j) ), dyl(j) ) / d2r
703  enddo
704  enddo
705  !$acc end kernels
706  !$acc kernels
707  dzsfc_dxy(:,ja,2) = 0.0_rp
708  !$acc end kernels
709 
710  !$acc kernels
711  slope(:,:) = max( abs(dzsfc_dxy(:,:,1)), abs(dzsfc_dxy(:,:,2)) )
712  !$acc end kernels
713  call statistics_horizontal_max( ia, is, ie, ja, js, je, &
714  slope(:,:), maxslope )
715 
716  log_progress(*) 'maximum slope [deg] : ', maxslope
717 
718  if( maxslope < cnvtopo_smooth_maxslope_limit ) exit
719 
720  call statistics_detail( ia, is, ie, ja, js, je, 2, &
721  varname(:), dzsfc_dxy(:,:,:) )
722 
723  select case( cnvtopo_smooth_type )
724  case( 'GAUSSIAN' )
725 
726  ! 3 by 3 gaussian filter
727  !$omp parallel do
728  !$acc kernels
729  !$acc loop collapse(2) independent
730  do j = js, je
731  do i = is, ie
732  zsfc(i,j) = ( 0.2500_rp * zsfc(i ,j ) &
733  + 0.1250_rp * zsfc(i-1,j ) &
734  + 0.1250_rp * zsfc(i+1,j ) &
735  + 0.1250_rp * zsfc(i ,j-1) &
736  + 0.1250_rp * zsfc(i ,j+1) &
737  + 0.0625_rp * zsfc(i-1,j-1) &
738  + 0.0625_rp * zsfc(i+1,j-1) &
739  + 0.0625_rp * zsfc(i-1,j+1) &
740  + 0.0625_rp * zsfc(i+1,j+1) )
741  enddo
742  enddo
743  !$acc end kernels
744 
745  case( 'LAPLACIAN' )
746 
747  !$omp parallel do
748  !$acc kernels
749  do j = js , je
750  do i = is-1, ie
751  flx_x(i,j) = zsfc(i+1,j) - zsfc(i,j)
752 ! FLX_TMP(i,j) = Zsfc(i+1,j) - Zsfc(i,j)
753  enddo
754  enddo
755  !$acc end kernels
756 !!$ call TOPOGRAPHY_fillhalo( FLX_TMP )
757 !!$ do j = JS , JE
758 !!$ do i = IS-1, IE
759 !!$ FLX_X(i,j) = - ( FLX_TMP(i+1,j) - FLX_TMP(i,j) )
760 !!$ enddo
761 !!$ enddo
762 
763  !$omp parallel do
764  !$acc kernels
765  do j = js-1, je
766  do i = is , ie
767  flx_y(i,j) = zsfc(i,j+1) - zsfc(i,j)
768 ! FLX_TMP(i,j) = Zsfc(i,j+1) - Zsfc(i,j)
769  enddo
770  enddo
771  !$acc end kernels
772 !!$ call TOPOGRAPHY_fillhalo( FLX_TMP )
773 !!$ do j = JS-1, JE
774 !!$ do i = IS , IE
775 !!$ FLX_Y(i,j) = - ( FLX_TMP(i,j+1) - FLX_TMP(i,j) )
776 !!$ enddo
777 !!$ enddo
778 
779 
780  if ( cnvtopo_smooth_local ) then
781  !$omp parallel do &
782  !$omp private(flag)
783  !$acc kernels
784  do j = js , je
785  do i = is-1, ie
786  flag = 0.5_rp &
787  + sign(0.5_rp, max( abs(dzsfc_dxy(i+1,j ,1)), &
788  abs(dzsfc_dxy(i ,j ,1)), &
789  abs(dzsfc_dxy(i-1,j ,1)), &
790  abs(dzsfc_dxy(i+1,j ,2)), &
791  abs(dzsfc_dxy(i+1,j-1,2)), &
792  abs(dzsfc_dxy(i ,j ,2)), &
793  abs(dzsfc_dxy(i ,j-1,2)) &
794  ) - cnvtopo_smooth_maxslope_limit )
795  flx_x(i,j) = flx_x(i,j) * flag
796  enddo
797  enddo
798  !$acc end kernels
799  !$omp parallel do &
800  !$omp private(flag)
801  !$acc kernels
802  do j = js-1, je
803  do i = is , ie
804  flag = 0.5_rp &
805  + sign(0.5_rp, max( abs(dzsfc_dxy(i ,j+1,2)), &
806  abs(dzsfc_dxy(i ,j ,2)), &
807  abs(dzsfc_dxy(i ,j-1,2)), &
808  abs(dzsfc_dxy(i ,j+1,1)), &
809  abs(dzsfc_dxy(i-1,j+1,1)), &
810  abs(dzsfc_dxy(i ,j ,1)), &
811  abs(dzsfc_dxy(i-1,j ,1)) &
812  ) - cnvtopo_smooth_maxslope_limit )
813  flx_y(i,j) = flx_y(i,j) * flag
814  enddo
815  enddo
816  !$acc end kernels
817  endif
818 
819  !$omp parallel do
820  !$acc kernels
821  do j = js, je
822  do i = is, ie
823  zsfc(i,j) = zsfc(i,j) &
824  + 0.1_rp * ( ( flx_x(i,j) - flx_x(i-1,j) ) &
825  + ( flx_y(i,j) - flx_y(i,j-1) ) )
826  enddo
827  enddo
828  !$acc end kernels
829 
830  case default
831  log_error("CNVTOPO_smooth",*) 'Invalid smoothing type'
832  call prc_abort
833  end select
834 
835  if ( cnvtopo_smooth_trim_ocean ) then
836  !$omp parallel do
837  !$acc kernels
838  do j = js, je
839  do i = is, ie
840  zsfc(i,j) = sign( max( zsfc(i,j) * topo_sign(i,j), 0.0_rp ), topo_sign(i,j) )
841  end do
842  end do
843  !$acc end kernels
844  end if
845 
846  enddo
847 
848  if ( ite > cnvtopo_smooth_itelim ) then
849  log_error("CNVTOPO_smooth",*) 'Smoothing did not converge until ', cnvtopo_smooth_itelim,' times of iteration.'
850 
851  log_error_cont(*) 'Please try different parameters of PARAM_CNVTOPO.'
852  log_error_cont(*) '- Number limit of iteration (CNVTOPO_smooth_itelim) = ', cnvtopo_smooth_itelim
853  log_error_cont(*) '- Maximum ratio of slope dZ/dX, dZ/dY (CNVTOPO_smooth_maxslope_ratio) = ', cnvtopo_smooth_maxslope_ratio
854  log_error_cont(*) ' Or, Maximum of slope with degree (CNVTOPO_smooth_maxslope) = ', cnvtopo_smooth_maxslope
855  log_error_cont(*) '- Smoothing type LAPLACIAN/GAUSSIAN/OFF (CNVTOPO_smooth_type) = ', trim(cnvtopo_smooth_type)
856  call prc_abort
857  else
858  log_newline
859  log_info("CNVTOPO_smooth",*) 'smoothing complete.'
860  endif
861 
862 
863 
864  if ( cnvtopo_smooth_hypdiff_niter > 0 ) then
865 
866  log_newline
867  log_info("CNVTOPO_smooth",*) 'Apply hyperdiffusion.'
868 
869  call filter_hyperdiff( ia, is, ie, ja, js, je, &
870  zsfc(:,:), & ! (inout)
871  cnvtopo_smooth_hypdiff_order, cnvtopo_smooth_hypdiff_niter, & ! (in)
872  limiter_sign = topo_sign(:,:) ) ! (in)
873 
874  call topography_fillhalo( zsfc=zsfc(:,:), fill_bnd=.true. )
875 
876  !$omp parallel do
877  !$acc kernels
878  do j = 1, ja
879  do i = 1, ia-1
880  dzsfc_dxy(i,j,1) = atan2( ( zsfc(i+1,j)-zsfc(i,j) ), dxl(i) ) / d2r
881  enddo
882  enddo
883  !$acc end kernels
884  !$acc kernels
885  dzsfc_dxy(ia,:,1) = 0.0_rp
886  !$acc end kernels
887  !$omp parallel do
888  !$acc kernels
889  do j = 1, ja-1
890  do i = 1, ia
891  dzsfc_dxy(i,j,2) = atan2( ( zsfc(i,j+1)-zsfc(i,j) ), dyl(j) ) / d2r
892  enddo
893  enddo
894  !$acc end kernels
895  !$acc kernels
896  dzsfc_dxy(:,ja,2) = 0.0_rp
897  !$acc end kernels
898 
899  !$acc kernels
900  slope(:,:) = max( abs(dzsfc_dxy(:,:,1)), abs(dzsfc_dxy(:,:,2)) )
901  !$acc end kernels
902  call statistics_horizontal_max( ia, is, ie, ja, js, je, &
903  slope(:,:), maxslope )
904 
905  log_info("CNVTOPO_smooth",*) 'maximum slope [deg] : ', maxslope
906 
907  end if
908 
909  call topography_fillhalo( zsfc=zsfc(:,:), fill_bnd=.true. )
910 
911  call statistics_detail( ia, is, ie, ja, js, je, 2, &
912  varname(:), dzsfc_dxy(:,:,:) )
913 
914 
915  if ( cnvtopo_smooth_trim_ocean ) then
916  !$acc exit data delete(TOPO_sign)
917  end if
918  !$acc end data
919 
920  log_newline
921 
922  return
923  end subroutine cnvtopo_smooth
924 
925 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:228
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:45
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:350
mod_cnvtopo::cnvtopo
subroutine, public cnvtopo
Driver.
Definition: mod_cnvtopo.F90:273
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:43
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:35
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:64
mod_copytopo::copytopo
subroutine, public copytopo(TOPO_child)
Setup and Main.
Definition: mod_copytopo.F90:79
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:45
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:37
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:334
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:63
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:55
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:39
mod_cnv2d::cnv2d_tile_init
subroutine, public cnv2d_tile_init(dtype, dlat, dlon, dir, catalogue, interp_type, interp_level, nmax)
Definition: mod_cnv2d.F90:124
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:33
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:43
scale_topography::topography_fillhalo
subroutine, public topography_fillhalo(Zsfc, FILL_BND)
HALO Communication.
Definition: scale_topography.F90:119
scale_io::io_fid_conf
integer, public io_fid_conf
Config file ID.
Definition: scale_io.F90:57
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:50
Generated by   doxygen 1.8.17