SCALE-RM
mod_cnvlanduse.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
9 !-------------------------------------------------------------------------------
11  !-----------------------------------------------------------------------------
12  !
13  !++ used modules
14  !
15  use scale_precision
16  use scale_stdio
17  use scale_prof
19  use scale_tracer
20  !-----------------------------------------------------------------------------
21  implicit none
22  private
23  !-----------------------------------------------------------------------------
24  !
25  !++ Public procedure
26  !
27  public :: cnvlanduse_setup
28  public :: cnvlanduse
29 
30  !-----------------------------------------------------------------------------
31  !
32  !++ Public parameters & variables
33  !
34  logical, public :: cnvlanduse_donothing
35  logical, public :: cnvlanduse_useglccv2 = .false.
36  logical, public :: cnvlanduse_uselu100m = .false.
37  logical, public :: cnvlanduse_usejibis = .false.
38 
39  !-----------------------------------------------------------------------------
40  !
41  !++ Private procedure
42  !
43  private :: cnvlanduse_glccv2
44  private :: cnvlanduse_lu100m
45  private :: cnvlanduse_jibis
46 
47  !-----------------------------------------------------------------------------
48  !
49  !++ Private parameters & variables
50  !
51  real(RP), private :: cnvlanduse_unittile_ddeg ! dx for unit tile [deg]
52  real(RP), private :: cnvlanduse_oversampling_factor = 2.0_rp ! factor of min. dx against the unit tile
53 
54  !-----------------------------------------------------------------------------
55 contains
56  !-----------------------------------------------------------------------------
58  subroutine cnvlanduse_setup
59  use scale_process, only: &
61  use scale_const, only: &
62  d2r => const_d2r
63  use scale_comm, only: &
64  comm_horizontal_min
65  use scale_grid_real, only: &
66  real_dlat, &
67  real_dlon
68  implicit none
69 
70  character(len=H_SHORT) :: CNVLANDUSE_name = 'NONE' ! keep backward compatibility
71 
72  namelist / param_cnvlanduse / &
73  cnvlanduse_name, &
77  cnvlanduse_unittile_ddeg, &
78  cnvlanduse_oversampling_factor
79 
80  real(RP) :: drad(ia,ja)
81  real(RP) :: drad_min
82 
83  integer :: ierr
84  !---------------------------------------------------------------------------
85 
86  if( io_l ) write(io_fid_log,*)
87  if( io_l ) write(io_fid_log,*) '++++++ Module[convert landuseindex] / Categ[preprocess] / Origin[SCALE-RM]'
88  !--- read namelist
89  rewind(io_fid_conf)
90  read(io_fid_conf,nml=param_cnvlanduse,iostat=ierr)
91  if( ierr < 0 ) then !--- missing
92  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
93  elseif( ierr > 0 ) then !--- fatal error
94  write(*,*) 'xxx Not appropriate names in namelist PARAM_CNVLANDUSE. Check!'
95  call prc_mpistop
96  endif
97  if( io_lnml ) write(io_fid_log,nml=param_cnvlanduse)
98 
99  select case(cnvlanduse_name)
100  case('NONE')
101  ! do nothing
102  case('GLCCv2')
103  cnvlanduse_useglccv2 = .true.
104  cnvlanduse_uselu100m = .false.
105  cnvlanduse_usejibis = .false.
106  case('LU100M')
107  cnvlanduse_useglccv2 = .false.
108  cnvlanduse_uselu100m = .true.
109  cnvlanduse_usejibis = .false.
110  case('COMBINE')
111  cnvlanduse_useglccv2 = .true.
112  cnvlanduse_uselu100m = .true.
113  cnvlanduse_usejibis = .true.
114  case default
115  write(*,*) ' xxx Unsupported TYPE:', trim(cnvlanduse_name)
116  call prc_mpistop
117  endselect
118 
119  cnvlanduse_donothing = .true.
120 
121  if ( cnvlanduse_useglccv2 ) then
122  cnvlanduse_donothing = .false.
123  if( io_l ) write(io_fid_log,*) '*** Use GLCC ver.2, global 30 arcsec. data'
124  if ( cnvlanduse_uselu100m ) then
125  if( io_l ) write(io_fid_log,*) '*** Use KSJ landuse 100m data for Japan region'
126  if( io_l ) write(io_fid_log,*) '*** Overwrite Japan region'
127  if ( cnvlanduse_usejibis ) then
128  if( io_l ) write(io_fid_log,*) '*** Use J-IBIS map 100m data for Japan region'
129  if( io_l ) write(io_fid_log,*) '*** Overwrite Japan region (PFT only)'
130  endif
131  endif
132  elseif ( cnvlanduse_uselu100m ) then
133  cnvlanduse_donothing = .false.
134  if( io_l ) write(io_fid_log,*) '*** Use KSJ landuse 100m data, Japan region only'
135  if ( cnvlanduse_usejibis ) then
136  if( io_l ) write(io_fid_log,*) '*** Use J-IBIS map 100m data for Japan region'
137  if( io_l ) write(io_fid_log,*) '*** Overwrite Japan region (PFT only)'
138  endif
139  endif
140 
141  if ( cnvlanduse_donothing ) then
142  if( io_l ) write(io_fid_log,*) '*** Do nothing for landuse index'
143  else
144  drad(:,:) = min( real_dlat(:,:), real_dlon(:,:) )
145  call comm_horizontal_min( drad_min, drad(:,:) )
146 
147  if ( cnvlanduse_unittile_ddeg > 0.0_rp ) then
148  cnvlanduse_oversampling_factor = ( drad_min / d2r ) / cnvlanduse_unittile_ddeg
149  endif
150  cnvlanduse_oversampling_factor = max( 1.0_rp, cnvlanduse_oversampling_factor )
151  cnvlanduse_unittile_ddeg = ( drad_min / d2r ) / cnvlanduse_oversampling_factor
152 
153  if( io_l ) write(io_fid_log,*) '*** The size of tile [deg] = ', cnvlanduse_unittile_ddeg
154  if( io_l ) write(io_fid_log,*) '*** oversampling factor = ', cnvlanduse_oversampling_factor
155  endif
156 
157  return
158  end subroutine cnvlanduse_setup
159 
160  !-----------------------------------------------------------------------------
162  subroutine cnvlanduse
163  use scale_process, only: &
165  use scale_landuse, only: &
168  implicit none
169  !---------------------------------------------------------------------------
170 
171  if ( cnvlanduse_donothing ) then
172  if( io_l ) write(io_fid_log,*)
173  if( io_l ) write(io_fid_log,*) '++++++ SKIP CONVERT LANDUSE DATA ++++++'
174  else
175  if( io_l ) write(io_fid_log,*)
176  if( io_l ) write(io_fid_log,*) '++++++ START CONVERT LANDUSE DATA ++++++'
177 
178  if ( cnvlanduse_useglccv2 ) then
179  call cnvlanduse_glccv2
180  endif
181 
182  if ( cnvlanduse_uselu100m ) then
183  call cnvlanduse_lu100m
184  endif
185 
186  if ( cnvlanduse_usejibis ) then
187  call cnvlanduse_jibis
188  endif
189 
190  ! calculate landuse factors
191  call landuse_calc_fact
192 
193  if( io_l ) write(io_fid_log,*) '++++++ END CONVERT LANDUSE DATA ++++++'
194 
195  ! output landuse file
196  call landuse_write
197  endif
198 
199  return
200  end subroutine cnvlanduse
201 
202  !-----------------------------------------------------------------------------
204  subroutine cnvlanduse_glccv2
205  use scale_process, only: &
207  use scale_const, only: &
208  radius => const_radius, &
209  eps => const_eps, &
210  d2r => const_d2r
211  use scale_landuse, only: &
219  use scale_grid_real, only: &
220  real_laty, &
221  real_lonx
222  implicit none
223 
224  character(len=H_LONG) :: GLCCv2_IN_CATALOGUE = ''
225  character(len=H_LONG) :: GLCCv2_IN_DIR = ''
226  real(RP) :: limit_urban_fraction = 1.0_rp
227 
228  namelist / param_cnvlanduse_glccv2 / &
229  glccv2_in_catalogue, &
230  glccv2_in_dir, &
231  limit_urban_fraction
232 
233  ! data catalogue list
234  integer, parameter :: TILE_nlim = 100
235  integer :: TILE_nmax
236  real(RP) :: TILE_LATS (tile_nlim)
237  real(RP) :: TILE_LATE (tile_nlim)
238  real(RP) :: TILE_LONS (tile_nlim)
239  real(RP) :: TILE_LONE (tile_nlim)
240  character(len=H_LONG) :: TILE_fname(tile_nlim)
241 
242  ! GLCCv2 data
243  integer, parameter :: isize_orig = 3600
244  integer(1) :: TILE_LANDUSE_orig(isize_orig,isize_orig)
245  real(RP) :: TILE_DLAT_orig, TILE_DLON_orig
246 
247  ! GLCCv2 data (oversampling)
248  integer :: ios
249  integer :: isize
250  integer(1), allocatable :: TILE_LANDUSE(:,:)
251  real(RP), allocatable :: TILE_LATH (:)
252  real(RP), allocatable :: TILE_LONH (:)
253  real(RP) :: TILE_DLAT, TILE_DLON
254  real(RP) :: area, area_fraction
255 
256  integer :: iloc
257  integer :: jloc
258  real(RP) :: ifrac_l ! fraction for iloc
259  real(RP) :: jfrac_b ! fraction for jloc
260 
261  real(RP) :: DOMAIN_LATS, DOMAIN_LATE
262  real(RP) :: DOMAIN_LONS, DOMAIN_LONE
263 
264  ! tentative: LANDUSE_PFT_nmax is assumed to be 4
265  real(RP) :: categ_sum(ia,ja,-2:landuse_pft_nmax)
266 
267  integer :: lookuptable(1:25)
268  data lookuptable / 0, & ! 1 Urban and Built-Up Land -> 0 urban
269  7, & ! 2 Dryland Cropland and Pasture -> 7 Dryland Cropland and Pasture
270  8, & ! 3 Irrigated Cropland and Pasture -> 8 Irrigated Cropland and Pasture
271  9, & ! 4 Mixed Cropland and Pasture -> 9 Mixed Cropland and Pasture
272  5, & ! 5 Cropland/Grassland Mosaic -> 5 Cropland/Grassland Mosaic
273  6, & ! 6 Cropland/Woodland Mosaic -> 6 Cropland/Woodland Mosaic
274  2, & ! 7 Grassland -> 2 Grassland
275  3, & ! 8 Shrubland -> 3 Shrubland
276  4, & ! 9 Mixed Shrubland/Grassland -> 4 Mixed Shrubland/Grassland
277  4, & ! 10 Savanna -> 4 Mixed Shrubland/Grassland
278  11, & ! 11 Deciduous Broadleaf Forest -> 11 Deciduous Broadleaf Forest
279  12, & ! 12 Deciduous Needleleaf Forest -> 12 Deciduous Needleleaf Forest
280  13, & ! 13 Evergreen Broadleaf Forest -> 13 Deciduous Broadleaf Forest
281  14, & ! 14 Evergreen Needleleaf Forest -> 14 Deciduous Needleleaf Forest
282  15, & ! 15 Mixed Forest -> 15 Mixed Forest
283  -2, & ! 16 Water Bodies -> -2 Lake/River
284  10, & ! 17 Herbaceous Wetland -> 10 Paddy
285  10, & ! 18 Wooded Wetland -> 10 Paddy
286  1, & ! 19 Barren or Sparsely Vegetated -> 1 Bare Ground
287  1, & ! 20 Herbaceous Tundra -> 1 Bare Ground
288  1, & ! 21 Wooded Tundra -> 1 Bare Ground
289  1, & ! 22 Mixed Tundra -> 1 Bare Ground
290  1, & ! 23 Bare Ground Tundra -> 1 Bare Ground
291  1, & ! 24 Snow or Ice -> 1 Bare Ground
292  -1 / ! 25+ Sea Surface -> -1 Sea Surface
293 
294  real(RP) :: categ_pftsum, allsum
295  real(RP) :: PFT (landuse_pft_nmax)
296  integer :: PFT_idx(landuse_pft_nmax)
297  real(RP) :: temp
298  integer :: temp_idx
299 
300  character(len=H_LONG) :: fname
301 
302  real(RP) :: zerosw
303  logical :: hit_lat, hit_lon
304  integer :: index
305  integer :: fid, ierr
306  integer :: i, j, ii, jj, iii, jjj, t, p, pp
307  !---------------------------------------------------------------------------
308 
309  if( io_l ) write(io_fid_log,*)
310  if( io_l ) write(io_fid_log,*) '+++ Module[GLCCv2]/Categ[CNVLANDUSE]'
311 
312  !--- read namelist
313  rewind(io_fid_conf)
314  read(io_fid_conf,nml=param_cnvlanduse_glccv2,iostat=ierr)
315  if( ierr < 0 ) then !--- missing
316  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
317  elseif( ierr > 0 ) then !--- fatal error
318  write(*,*) 'xxx Not appropriate names in namelist PARAM_CNVLANDUSE_GLCCv2. Check!'
319  call prc_mpistop
320  endif
321  if( io_lnml ) write(io_fid_log,nml=param_cnvlanduse_glccv2)
322 
323  do p = -2, landuse_pft_nmax
324  do j = 1, ja
325  do i = 1, ia
326  categ_sum(i,j,p) = 0.0_rp
327  enddo
328  enddo
329  enddo
330 
331  domain_lats = minval(real_laty(:,:))
332  domain_late = maxval(real_laty(:,:))
333  domain_lons = minval(real_lonx(:,:))
334  domain_lone = maxval(real_lonx(:,:))
335 
336  ios = nint( 30.0_rp / 60.0_rp / 60.0_rp / cnvlanduse_unittile_ddeg - 0.5_rp ) + 1
337  isize = isize_orig * ios
338 
339  allocate( tile_landuse(isize,isize) )
340  allocate( tile_lath(0:isize) )
341  allocate( tile_lonh(0:isize) )
342 
343  if( io_l ) write(io_fid_log,*) '*** Oversampling orig = ', isize_orig, ', use = ', isize
344 
345  tile_dlat_orig = 30.0_rp / 60.0_rp / 60.0_rp * d2r
346  tile_dlon_orig = 30.0_rp / 60.0_rp / 60.0_rp * d2r
347  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT :', tile_dlat_orig/d2r
348  if( io_l ) write(io_fid_log,*) '*** TILE_DLON :', tile_dlon_orig/d2r
349 
350  tile_dlat = tile_dlat_orig / ios
351  tile_dlon = tile_dlon_orig / ios
352  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT (OS) :', tile_dlat/d2r
353  if( io_l ) write(io_fid_log,*) '*** TILE_DLON (OS) :', tile_dlon/d2r
354 
355  !---< READ from external files >---
356 
357  ! catalogue file
358  fname = trim(glccv2_in_dir)//'/'//trim(glccv2_in_catalogue)
359 
360  if( io_l ) write(io_fid_log,*)
361  if( io_l ) write(io_fid_log,*) '+++ Input catalogue file:', trim(fname)
362 
363  fid = io_get_available_fid()
364  open( fid, &
365  file = trim(fname), &
366  form = 'formatted', &
367  status = 'old', &
368  iostat = ierr )
369 
370  if ( ierr /= 0 ) then
371  write(*,*) 'xxx catalogue file not found!', trim(fname)
372  call prc_mpistop
373  endif
374 
375  do t = 1, tile_nlim
376  read(fid,*,iostat=ierr) index, tile_lats(t), tile_late(t), & ! South->North
377  tile_lons(t), tile_lone(t), & ! WEST->EAST
378  tile_fname(t)
379  if ( ierr /= 0 ) exit
380  enddo
381 
382  tile_nmax = t - 1
383  close(fid)
384 
385  ! data file
386  do t = 1, tile_nmax
387  hit_lat = .false.
388  hit_lon = .false.
389 
390  if ( ( tile_lats(t)*d2r >= domain_lats .AND. tile_lats(t)*d2r < domain_late ) &
391  .OR. ( tile_late(t)*d2r >= domain_lats .AND. tile_late(t)*d2r < domain_late ) ) then
392  hit_lat = .true.
393  endif
394 
395  if ( ( domain_lats >= tile_lats(t)*d2r .AND. domain_lats < tile_late(t)*d2r ) &
396  .OR. ( domain_late >= tile_lats(t)*d2r .AND. domain_late < tile_late(t)*d2r ) ) then
397  hit_lat = .true.
398  endif
399 
400  if ( ( tile_lons(t)*d2r >= domain_lons .AND. tile_lons(t)*d2r < domain_lone ) &
401  .OR. ( tile_lone(t)*d2r >= domain_lons .AND. tile_lone(t)*d2r < domain_lone ) ) then
402  hit_lon = .true.
403  endif
404 
405  if ( ( domain_lons >= tile_lons(t)*d2r .AND. domain_lons < tile_lone(t)*d2r ) &
406  .OR. ( domain_lone >= tile_lons(t)*d2r .AND. domain_lone < tile_lone(t)*d2r ) ) then
407  hit_lon = .true.
408  endif
409 
410  if ( hit_lat .AND. hit_lon ) then
411  fname = trim(glccv2_in_dir)//'/'//trim(tile_fname(t))
412 
413  if( io_l ) write(io_fid_log,*)
414  if( io_l ) write(io_fid_log,*) '+++ Input data file :', trim(fname)
415  if( io_l ) write(io_fid_log,*) '*** Domain (LAT) :', domain_lats/d2r, domain_late/d2r
416  if( io_l ) write(io_fid_log,*) '*** (LON) :', domain_lons/d2r, domain_lone/d2r
417  if( io_l ) write(io_fid_log,*) '*** Tile (LAT) :', tile_lats(t), tile_late(t)
418  if( io_l ) write(io_fid_log,*) '*** (LON) :', tile_lons(t), tile_lone(t)
419 
420  fid = io_get_available_fid()
421  open( fid, &
422  file = trim(fname), &
423  form = 'unformatted', &
424  access = 'direct', &
425  status = 'old', &
426  recl = isize_orig*isize_orig*1, &
427  iostat = ierr )
428 
429  if ( ierr /= 0 ) then
430  write(*,*) 'xxx data file not found!'
431  call prc_mpistop
432  endif
433 
434  read(fid,rec=1) tile_landuse_orig(:,:)
435  close(fid)
436 
437  ! oversampling
438  do jj = 1, isize_orig
439  do ii = 1, isize_orig
440  do j = 1, ios
441  do i = 1, ios
442  jjj = (jj-1) * ios + j
443  iii = (ii-1) * ios + i
444 
445  tile_landuse(iii,jjj) = tile_landuse_orig(ii,jj)
446  enddo
447  enddo
448  enddo
449  enddo
450 
451  tile_lath(0) = tile_lats(t) * d2r
452  do jj = 1, isize
453  tile_lath(jj) = tile_lath(jj-1) + tile_dlat
454 ! if( IO_L ) write(IO_FID_LOG,*) jj, TILE_LATH(jj)
455  enddo
456 
457  tile_lonh(0) = tile_lons(t) * d2r
458  do ii = 1, isize
459  tile_lonh(ii) = tile_lonh(ii-1) + tile_dlon
460 ! if( IO_L ) write(IO_FID_LOG,*) ii, TILE_LONH(ii)
461  enddo
462 
463  ! match and calc fraction
464  do jj = 1, isize
465  do ii = 1, isize
466 
467  iloc = 1 ! Z_sfc(1,1) is used for dummy grid
468  ifrac_l = 1.0_rp
469 
470  jloc = 1 ! Z_sfc(1,1) is used for dummy grid
471  jfrac_b = 1.0_rp
472 
473  if ( tile_lonh(ii-1) < domain_lons &
474  .OR. tile_lonh(ii-1) >= domain_lone &
475  .OR. tile_lath(jj-1) < domain_lats &
476  .OR. tile_lath(jj-1) >= domain_late ) then
477  cycle
478  endif
479 
480  jloop: do j = js-1, je+1
481  iloop: do i = is-1, ie+1
482  if ( tile_lonh(ii-1) >= real_lonx(i-1,j ) &
483  .AND. tile_lonh(ii-1) < real_lonx(i ,j ) &
484  .AND. tile_lath(jj-1) >= real_laty(i ,j-1) &
485  .AND. tile_lath(jj-1) < real_laty(i ,j ) ) then
486 
487  iloc = i
488  ifrac_l = min( real_lonx(i,j)-tile_lonh(ii-1), tile_dlon ) / tile_dlon
489 
490  jloc = j
491  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
492  exit jloop
493 
494  endif
495  enddo iloop
496  enddo jloop
497 
498  if( iloc == 1 .AND. jloc == 1 ) cycle
499 
500  area = radius * radius * tile_dlon * ( sin(tile_lath(jj))-sin(tile_lath(jj-1)) )
501 
502  pp = min( max( int(tile_landuse(ii,jj),kind=4), 0 ), 25 )
503  p = lookuptable(pp)
504 
505  !if( IO_L ) write(IO_FID_LOG,*) ii, jj, iloc, jloc, p, pp
506 
507  area_fraction = ( ifrac_l) * ( jfrac_b) * area
508  categ_sum(iloc ,jloc ,p) = categ_sum(iloc ,jloc ,p) + area_fraction
509 
510  area_fraction = (1.0_rp-ifrac_l) * ( jfrac_b) * area
511  categ_sum(iloc+1,jloc ,p) = categ_sum(iloc+1,jloc ,p) + area_fraction
512 
513  area_fraction = ( ifrac_l) * (1.0_rp-jfrac_b) * area
514  categ_sum(iloc ,jloc+1,p) = categ_sum(iloc ,jloc+1,p) + area_fraction
515 
516  area_fraction = (1.0_rp-ifrac_l) * (1.0_rp-jfrac_b) * area
517  categ_sum(iloc+1,jloc+1,p) = categ_sum(iloc+1,jloc+1,p) + area_fraction
518  enddo
519  enddo
520 
521  endif
522  enddo ! tile loop
523 
524  do j = js, je
525  do i = is, ie
526 ! area = RADIUS * RADIUS * (REAL_LONX(i,j)-REAL_LONX(i-1,j)) * ( sin(REAL_LATY(i,j))-sin(REAL_LATY(i,j-1)) )
527 ! allsum = categ_sum(i,j,-2) + categ_sum(i,j,-1) + categ_sum(i,j,0) + categ_pftsum
528 ! if ( abs(allsum/area-1.0_RP) > EPS ) then
529 ! if( IO_L ) write(IO_FID_LOG,*) i,j,allsum/area
530 ! endif
531 
532  do p = 1, landuse_pft_nmax
533  pft(p) = categ_sum(i,j,p)
534  pft_idx(p) = p
535  enddo
536 
537  ! bubble sort
538  do p = 1, landuse_pft_nmax-1
539  do pp = p+1, landuse_pft_nmax
540  if ( pft(pp) > pft(p) ) then
541  temp = pft(p)
542  pft(p) = pft(pp)
543  pft(pp) = temp
544 
545  temp_idx = pft_idx(p)
546  pft_idx(p) = pft_idx(pp)
547  pft_idx(pp) = temp_idx
548  endif
549  enddo
550  enddo
551 
552  categ_pftsum = sum( pft(:) )
553 
554  ! land fraction : 1 - ocean / total
555  allsum = categ_sum(i,j,-2) + categ_sum(i,j,-1) + categ_sum(i,j,0) + categ_pftsum
556  zerosw = 0.5_rp - sign( 0.5_rp, allsum-eps )
557  landuse_frac_land(i,j) = 1.0_rp-zerosw - categ_sum(i,j,-1) * ( 1.0_rp-zerosw ) / ( allsum-zerosw )
558 
559  ! lake fraction : lake / ( total - ocean )
560  allsum = categ_sum(i,j,-2) + categ_sum(i,j,0) + categ_pftsum
561  zerosw = 0.5_rp - sign( 0.5_rp, allsum-eps )
562  landuse_frac_lake(i,j) = categ_sum(i,j,-2) * ( 1.0_rp-zerosw ) / ( allsum-zerosw )
563 
564  ! urban fraction : urban / ( total - ocean - lake )
565  allsum = categ_sum(i,j,0) + categ_pftsum
566  zerosw = 0.5_rp - sign( 0.5_rp, allsum-eps )
567  landuse_frac_urban(i,j) = categ_sum(i,j,0) * ( 1.0_rp-zerosw ) / ( allsum-zerosw )
568 
569  ! PFT fraction : PFT / sum( PFT(1:mosaic) )
570  allsum = sum( pft(1:landuse_pft_mosaic) )
571  if ( abs(allsum) > eps ) then
572  do p = 1, landuse_pft_mosaic
573  landuse_frac_pft(i,j,p) = pft(p) / allsum
574  landuse_index_pft(i,j,p) = pft_idx(p)
575  enddo
576  ! if no second PFT, set to same as PFT1
577  if ( abs(landuse_frac_pft(i,j,1)-1.0_rp) <= eps ) then
578  landuse_frac_pft(i,j,:) = 0.0_rp
579  landuse_frac_pft(i,j,1) = 1.0_rp
580  landuse_index_pft(i,j,:) = pft_idx(1)
581  endif
582  else ! if no PFT, set to bare ground
583  landuse_frac_pft(i,j,:) = 0.0_rp
584  landuse_frac_pft(i,j,1) = 1.0_rp
585  landuse_index_pft(i,j,:) = 1
586  endif
587 
588  enddo
589  enddo
590 
591  if ( limit_urban_fraction < 1.0_rp ) then
592  do j = js, je
593  do i = is, ie
594  if ( landuse_frac_urban(i,j) == 1.0_rp ) then ! if no PFT, set to grassland
595  landuse_frac_pft(i,j,:) = 0.0_rp
596  landuse_frac_pft(i,j,1) = 1.0_rp
597  landuse_index_pft(i,j,:) = 2
598  endif
599  landuse_frac_urban(i,j) = min( landuse_frac_urban(i,j), limit_urban_fraction )
600  enddo
601  enddo
602  endif
603 
604  return
605  end subroutine cnvlanduse_glccv2
606 
607  !-----------------------------------------------------------------------------
609  subroutine cnvlanduse_lu100m
610  use scale_process, only: &
612  use scale_const, only: &
613  radius => const_radius, &
614  eps => const_eps, &
615  d2r => const_d2r
616  use scale_landuse, only: &
624  use scale_grid_real, only: &
625  real_laty, &
626  real_lonx
627  implicit none
628 
629  character(len=H_LONG) :: LU100M_IN_CATALOGUE = ''
630  character(len=H_LONG) :: LU100M_IN_DIR = ''
631  real(RP) :: limit_urban_fraction = 1.0_rp
632 
633  namelist / param_cnvlanduse_lu100m / &
634  lu100m_in_catalogue, &
635  lu100m_in_dir, &
636  limit_urban_fraction
637 
638  ! data catalogue list
639  integer, parameter :: TILE_nlim = 1000
640  integer :: TILE_nmax
641  real(RP) :: TILE_LATS (tile_nlim)
642  real(RP) :: TILE_LATE (tile_nlim)
643  real(RP) :: TILE_LONS (tile_nlim)
644  real(RP) :: TILE_LONE (tile_nlim)
645  character(len=H_LONG) :: TILE_fname(tile_nlim)
646 
647  ! LU100M data
648  integer, parameter :: isize_orig = 800
649  real(SP) :: TILE_LANDUSE_orig(isize_orig,isize_orig)
650  real(RP) :: TILE_DLAT_orig, TILE_DLON_orig
651 
652  ! LU100M data (oversampling)
653  integer :: ios
654  integer :: isize
655  real(RP), allocatable :: TILE_LANDUSE(:,:)
656  real(RP), allocatable :: TILE_LATH (:)
657  real(RP), allocatable :: TILE_LONH (:)
658  real(RP) :: TILE_DLAT, TILE_DLON
659  real(RP) :: area, area_fraction
660 
661  integer :: iloc
662  integer :: jloc
663  real(RP) :: ifrac_l ! fraction for iloc
664  real(RP) :: jfrac_b ! fraction for jloc
665 
666  real(RP) :: DOMAIN_LATS, DOMAIN_LATE
667  real(RP) :: DOMAIN_LONS, DOMAIN_LONE
668 
669  ! tentative: LANDUSE_PFT_nmax is assumed to be 4
670  real(RP) :: categ_sum(ia,ja,-2:landuse_pft_nmax)
671 
672  integer :: lookuptable(0:16)
673  data lookuptable / -1, & ! 0 missing -> -1 Sea Surface
674  10, & ! 1 paddy -> 10 Paddy
675  9, & ! 2 cropland -> 9 Mixed Cropland and Pasture
676  1, & ! 3 UNDEF -> 1 Bare Ground
677  1, & ! 4 UNDEF -> 1 Bare Ground
678  11, & ! 5 forest -> 11 Deciduous Broadleaf Forest
679  1, & ! 6 bareground -> 1 Bare Ground
680  0, & ! 7 urban building -> 0 Urban and Built-up Land
681  1, & ! 8 UNDEF -> 1 Bare Ground
682  0, & ! 9 motorway -> 0 Urban and Built-up Land
683  0, & ! 10 urban ground -> 0 Urban and Built-up Land
684  -2, & ! 11 lake,river -> -2 Lake/River
685  1, & ! 12 UNDEF -> 1 Bare Ground
686  1, & ! 13 UNDEF -> 1 Bare Ground
687  1, & ! 14 seashore -> 1 Bare Ground
688  -1, & ! 15 ocean -> -1 Sea Surface
689  2 / ! 16 golf course -> 2 Grassland
690 
691  real(RP) :: categ_pftsum, allsum
692  real(RP) :: PFT (landuse_pft_nmax)
693  integer :: PFT_idx(landuse_pft_nmax)
694  real(RP) :: temp
695  integer :: temp_idx
696  real(RP) :: frac, mask
697 
698  character(len=H_LONG) :: fname
699 
700  real(RP) :: zerosw
701  logical :: hit_lat, hit_lon
702  integer :: index
703  integer :: fid, ierr
704  integer :: i, j, ii, jj, iii, jjj, t, p, pp
705  !---------------------------------------------------------------------------
706 
707  if( io_l ) write(io_fid_log,*)
708  if( io_l ) write(io_fid_log,*) '+++ Module[LU100M]/Categ[CNVLANDUSE]'
709 
710  !--- read namelist
711  rewind(io_fid_conf)
712  read(io_fid_conf,nml=param_cnvlanduse_lu100m,iostat=ierr)
713  if( ierr < 0 ) then !--- missing
714  if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
715  elseif( ierr > 0 ) then !--- fatal error
716  write(*,*) 'xxx Not appropriate names in namelist PARAM_CNVLANDUSE_LU100M. Check!'
717  call prc_mpistop
718  endif
719  if( io_lnml ) write(io_fid_log,nml=param_cnvlanduse_lu100m)
720 
721  do p = -2, landuse_pft_nmax
722  do j = 1, ja
723  do i = 1, ia
724  categ_sum(i,j,p) = 0.0_rp
725  enddo
726  enddo
727  enddo
728 
729  domain_lats = minval(real_laty(:,:))
730  domain_late = maxval(real_laty(:,:))
731  domain_lons = minval(real_lonx(:,:))
732  domain_lone = maxval(real_lonx(:,:))
733 
734  ios = nint( 5.0_rp / 60.0_rp / 100.0_rp / cnvlanduse_unittile_ddeg - 0.5_rp ) + 1
735  isize = isize_orig * ios
736 
737  allocate( tile_landuse(isize,isize) )
738  allocate( tile_lath(0:isize) )
739  allocate( tile_lonh(0:isize) )
740 
741  if( io_l ) write(io_fid_log,*) '*** Oversampling orig = ', isize_orig, ', use = ', isize
742 
743  tile_dlat_orig = 5.0_rp / 60.0_rp / 100.0_rp * d2r
744  tile_dlon_orig = 7.5_rp / 60.0_rp / 100.0_rp * d2r
745  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT :', tile_dlat_orig/d2r
746  if( io_l ) write(io_fid_log,*) '*** TILE_DLON :', tile_dlon_orig/d2r
747 
748  tile_dlat = tile_dlat_orig / ios
749  tile_dlon = tile_dlon_orig / ios
750  if( io_l ) write(io_fid_log,*) '*** TILE_DLAT (OS) :', tile_dlat/d2r
751  if( io_l ) write(io_fid_log,*) '*** TILE_DLON (OS) :', tile_dlon/d2r
752 
753  !---< READ from external files >---
754 
755  ! catalogue file
756  fname = trim(lu100m_in_dir)//'/'//trim(lu100m_in_catalogue)
757 
758  if( io_l ) write(io_fid_log,*)
759  if( io_l ) write(io_fid_log,*) '+++ Input catalogue file:', trim(fname)
760 
761  fid = io_get_available_fid()
762  open( fid, &
763  file = trim(fname), &
764  form = 'formatted', &
765  status = 'old', &
766  iostat = ierr )
767 
768  if ( ierr /= 0 ) then
769  write(*,*) 'xxx catalogue file not found!', trim(fname)
770  call prc_mpistop
771  endif
772 
773  do t = 1, tile_nlim
774  read(fid,*,iostat=ierr) index, tile_lats(t), tile_late(t), & ! South->North
775  tile_lons(t), tile_lone(t), & ! WEST->EAST
776  tile_fname(t)
777  if ( ierr /= 0 ) exit
778  enddo
779 
780  tile_nmax = t - 1
781  close(fid)
782 
783  ! data file
784  do t = 1, tile_nmax
785  hit_lat = .false.
786  hit_lon = .false.
787 
788  if ( ( tile_lats(t)*d2r >= domain_lats .AND. tile_lats(t)*d2r < domain_late ) &
789  .OR. ( tile_late(t)*d2r >= domain_lats .AND. tile_late(t)*d2r < domain_late ) ) then
790  hit_lat = .true.
791  endif
792 
793  if ( ( domain_lats >= tile_lats(t)*d2r .AND. domain_lats < tile_late(t)*d2r ) &
794  .OR. ( domain_late >= tile_lats(t)*d2r .AND. domain_late < tile_late(t)*d2r ) ) then
795  hit_lat = .true.
796  endif
797 
798  if ( ( tile_lons(t)*d2r >= domain_lons .AND. tile_lons(t)*d2r < domain_lone ) &
799  .OR. ( tile_lone(t)*d2r >= domain_lons .AND. tile_lone(t)*d2r < domain_lone ) ) then
800  hit_lon = .true.
801  endif
802 
803  if ( ( domain_lons >= tile_lons(t)*d2r .AND. domain_lons < tile_lone(t)*d2r ) &
804  .OR. ( domain_lone >= tile_lons(t)*d2r .AND. domain_lone < tile_lone(t)*d2r ) ) then
805  hit_lon = .true.
806  endif
807 
808  if ( hit_lat .AND. hit_lon ) then
809  fname = trim(lu100m_in_dir)//'/'//trim(tile_fname(t))
810 
811  if( io_l ) write(io_fid_log,*)
812  if( io_l ) write(io_fid_log,*) '+++ Input data file :', trim(fname)
813  if( io_l ) write(io_fid_log,*) '*** Domain (LAT) :', domain_lats/d2r, domain_late/d2r
814  if( io_l ) write(io_fid_log,*) '*** (LON) :', domain_lons/d2r, domain_lone/d2r
815  if( io_l ) write(io_fid_log,*) '*** Tile (LAT) :', tile_lats(t), tile_late(t)
816  if( io_l ) write(io_fid_log,*) '*** (LON) :', tile_lons(t), tile_lone(t)
817 
818  fid = io_get_available_fid()
819  open( fid, &
820  file = trim(fname), &
821  form = 'unformatted', &
822  access = 'direct', &
823  status = 'old', &
824  recl = isize_orig*isize_orig*4, &
825  iostat = ierr )
826 
827  if ( ierr /= 0 ) then
828  write(*,*) 'xxx data file not found!'
829  call prc_mpistop
830  endif
831 
832  read(fid,rec=1) tile_landuse_orig(:,:)
833  close(fid)
834 
835  ! oversampling
836  do jj = 1, isize_orig
837  do ii = 1, isize_orig
838  do j = 1, ios
839  do i = 1, ios
840  jjj = (jj-1) * ios + j
841  iii = (ii-1) * ios + i
842 
843  tile_landuse(iii,jjj) = real( TILE_LANDUSE_orig(ii,jj), kind=rp )
844  enddo
845  enddo
846  enddo
847  enddo
848 
849  tile_lath(0) = tile_lats(t) * d2r
850  do jj = 1, isize
851  tile_lath(jj) = tile_lath(jj-1) + tile_dlat
852 ! if( IO_L ) write(IO_FID_LOG,*) jj, TILE_LATH(jj)
853  enddo
854 
855  tile_lonh(0) = tile_lons(t) * d2r
856  do ii = 1, isize
857  tile_lonh(ii) = tile_lonh(ii-1) + tile_dlon
858 ! if( IO_L ) write(IO_FID_LOG,*) ii, TILE_LONH(ii)
859  enddo
860 
861  ! match and calc fraction
862  do jj = 1, isize
863  do ii = 1, isize
864 
865  iloc = 1 ! Z_sfc(1,1) is used for dummy grid
866  ifrac_l = 1.0_rp
867 
868  jloc = 1 ! Z_sfc(1,1) is used for dummy grid
869  jfrac_b = 1.0_rp
870 
871  if ( tile_lonh(ii-1) < domain_lons &
872  .OR. tile_lonh(ii-1) >= domain_lone &
873  .OR. tile_lath(jj-1) < domain_lats &
874  .OR. tile_lath(jj-1) >= domain_late ) then
875  cycle
876  endif
877 
878  jloop: do j = js-1, je+1
879  iloop: do i = is-1, ie+1
880  if ( tile_lonh(ii-1) >= real_lonx(i-1,j ) &
881  .AND. tile_lonh(ii-1) < real_lonx(i ,j ) &
882  .AND. tile_lath(jj-1) >= real_laty(i ,j-1) &
883  .AND. tile_lath(jj-1) < real_laty(i ,j ) ) then
884 
885  iloc = i
886  ifrac_l = min( real_lonx(i,j)-tile_lonh(ii-1), tile_dlon ) / tile_dlon
887 
888  jloc = j
889  jfrac_b = min( real_laty(i,j)-tile_lath(jj-1), tile_dlat ) / tile_dlat
890  exit jloop
891 
892  endif
893  enddo iloop
894  enddo jloop
895 
896  if( iloc == 1 .AND. jloc == 1 ) cycle
897 
898  area = radius * radius * tile_dlon * ( sin(tile_lath(jj))-sin(tile_lath(jj-1)) )
899 
900  pp = int( max( tile_landuse(ii,jj), 0.0_rp ) )
901  p = lookuptable(pp)
902 
903  !if( IO_L ) write(IO_FID_LOG,*) ii, jj, iloc, jloc, p, pp
904 
905  area_fraction = ( ifrac_l) * ( jfrac_b) * area
906  categ_sum(iloc ,jloc ,p) = categ_sum(iloc ,jloc ,p) + area_fraction
907 
908  area_fraction = (1.0_rp-ifrac_l) * ( jfrac_b) * area
909  categ_sum(iloc+1,jloc ,p) = categ_sum(iloc+1,jloc ,p) + area_fraction
910 
911  area_fraction = ( ifrac_l) * (1.0_rp-jfrac_b) * area
912  categ_sum(iloc ,jloc+1,p) = categ_sum(iloc ,jloc+1,p) + area_fraction
913 
914  area_fraction = (1.0_rp-ifrac_l) * (1.0_rp-jfrac_b) * area
915  categ_sum(iloc+1,jloc+1,p) = categ_sum(iloc+1,jloc+1,p) + area_fraction
916  enddo
917  enddo
918 
919  endif
920  enddo ! tile loop
921 
922  do j = js, je
923  do i = is, ie
924 ! area = RADIUS * RADIUS * (REAL_LONX(i,j)-REAL_LONX(i-1,j)) * ( sin(REAL_LATY(i,j))-sin(REAL_LATY(i,j-1)) )
925 ! allsum = categ_sum(i,j,-2) + categ_sum(i,j,-1) + categ_sum(i,j,0) + categ_pftsum
926 ! if ( abs(allsum/area-1.0_RP) > EPS ) then
927 ! if( IO_L ) write(IO_FID_LOG,*) i,j,allsum/area
928 ! endif
929 
930  do p = 1, landuse_pft_nmax
931  pft(p) = categ_sum(i,j,p)
932  pft_idx(p) = p
933  enddo
934 
935  ! bubble sort
936  do p = 1, landuse_pft_nmax-1
937  do pp = p+1, landuse_pft_nmax
938  if ( pft(pp) > pft(p) ) then
939  temp = pft(p)
940  pft(p) = pft(pp)
941  pft(pp) = temp
942 
943  temp_idx = pft_idx(p)
944  pft_idx(p) = pft_idx(pp)
945  pft_idx(pp) = temp_idx
946  endif
947  enddo
948  enddo
949 
950  categ_pftsum = sum( pft(:) )
951 
952  ! overwrite existing map (GLCCv2) by LU100
953  allsum = categ_sum(i,j,-2) + categ_sum(i,j,0) + categ_pftsum
954  mask = 0.5_rp + sign( 0.5_rp, allsum-eps ) ! if any land data is found, overwrite
955 
956  ! land fraction : 1 - ocean / total
957  allsum = categ_sum(i,j,-2) + categ_sum(i,j,-1) + categ_sum(i,j,0) + categ_pftsum
958  zerosw = 0.5_rp - sign( 0.5_rp, allsum-eps )
959  frac = ( 1.0_rp-zerosw ) - categ_sum(i,j,-1) * ( 1.0_rp-zerosw ) / ( allsum-zerosw )
960  landuse_frac_land(i,j) = ( mask ) * frac & ! overwrite
961  + ( 1.0_rp-mask ) * landuse_frac_land(i,j) ! keep existing value
962 
963  ! lake fraction : lake / ( total - ocean )
964  allsum = categ_sum(i,j,-2) + categ_sum(i,j,0) + categ_pftsum
965  zerosw = 0.5_rp - sign( 0.5_rp, allsum-eps )
966  frac = categ_sum(i,j,-2) * ( 1.0_rp-zerosw ) / ( allsum-zerosw )
967  landuse_frac_lake(i,j) = ( mask ) * frac & ! overwrite
968  + ( 1.0_rp-mask ) * landuse_frac_lake(i,j) ! keep existing value
969 
970  ! urban fraction : urban / ( total - ocean - lake )
971  allsum = categ_sum(i,j,0) + categ_pftsum
972  zerosw = 0.5_rp - sign( 0.5_rp, allsum-eps )
973  frac = categ_sum(i,j,0) * ( 1.0_rp-zerosw ) / ( allsum-zerosw )
974  landuse_frac_urban(i,j) = ( mask ) * frac & ! overwrite
975  + ( 1.0_rp-mask ) * landuse_frac_urban(i,j) ! keep existing value
976 
977  ! PFT fraction : PFT / sum( PFT(1:mosaic) )
978  allsum = sum( pft(1:landuse_pft_mosaic) )
979  if ( abs(allsum) > eps ) then
980  do p = 1, landuse_pft_mosaic
981  landuse_frac_pft(i,j,p) = ( mask ) * pft(p) / allsum & ! overwrite
982  + ( 1.0_rp-mask ) * landuse_frac_pft(i,j,p) ! keep existing value
983  landuse_index_pft(i,j,p) = ( mask ) * pft_idx(p) & ! overwrite
984  + ( 1.0_rp-mask ) * landuse_index_pft(i,j,p) ! keep existing value
985  enddo
986  ! if no second PFT, set to same as PFT1
987  if ( abs(landuse_frac_pft(i,j,1)-1.0_rp) <= eps ) then
988  landuse_frac_pft(i,j,:) = 0.0_rp
989  landuse_frac_pft(i,j,1) = 1.0_rp
990  landuse_index_pft(i,j,:) = pft_idx(1)
991  endif
992  else ! if no PFT, set to bare ground
993  do p = 1, landuse_pft_mosaic
994  landuse_frac_pft(i,j,p) = ( mask ) * 0.0_rp & ! overwrite
995  + ( 1.0_rp-mask ) * landuse_frac_pft(i,j,p) ! keep existing value
996  landuse_index_pft(i,j,p) = ( mask ) * 1 & ! overwrite
997  + ( 1.0_rp-mask ) * landuse_index_pft(i,j,p) ! keep existing value
998  enddo
999  landuse_frac_pft(i,j,1) = ( mask ) * 1.0_rp & ! overwrite
1000  + ( 1.0_rp-mask ) * landuse_frac_pft(i,j,1) ! keep existing value
1001  endif
1002 
1003  enddo
1004  enddo
1005 
1006  if ( limit_urban_fraction < 1.0_rp ) then
1007  do j = js, je
1008  do i = is, ie
1009  if ( landuse_frac_urban(i,j) == 1.0_rp ) then ! if no PFT, set to grassland
1010  landuse_frac_pft(i,j,:) = 0.0_rp
1011  landuse_frac_pft(i,j,1) = 1.0_rp
1012  landuse_index_pft(i,j,:) = 2
1013  endif
1014  landuse_frac_urban(i,j) = min( landuse_frac_urban(i,j), limit_urban_fraction )
1015  enddo
1016  enddo
1017  endif
1018 
1019  return
1020  end subroutine cnvlanduse_lu100m
1021 
1022  !-----------------------------------------------------------------------------
1024  subroutine cnvlanduse_jibis
1025  implicit none
1026  !---------------------------------------------------------------------------
1027 
1028  return
1029  end subroutine cnvlanduse_jibis
1030 
1031 end module mod_cnvlanduse
integer, public is
start point of inner domain: x, local
integer, public je
end point of inner domain: y, local
real(rp), dimension(:,:,:), allocatable, public landuse_frac_pft
fraction of PFT for each mosaic
subroutine, public prc_mpistop
Abort MPI.
logical, public io_l
output log or not? (this process)
Definition: scale_stdio.F90:59
real(rp), public const_radius
radius of the planet [m]
Definition: scale_const.F90:46
module STDIO
Definition: scale_stdio.F90:12
subroutine, public landuse_calc_fact
subroutine, public cnvlanduse
Driver.
real(rp), public const_d2r
degree to radian
Definition: scale_const.F90:35
real(rp), dimension(:,:), allocatable, public landuse_frac_urban
urban fraction
logical, public cnvlanduse_useglccv2
module grid index
module TRACER
integer, public ia
of x whole cells (local, with HALO)
integer function, public io_get_available_fid()
search & get available file ID
module GRID (real space)
module LANDUSE
real(rp), dimension(:,:), allocatable, public real_dlon
delta longitude
subroutine, public landuse_write
Write landuse data.
module COMMUNICATION
Definition: scale_comm.F90:23
integer, public js
start point of inner domain: y, local
module PROCESS
integer, public landuse_pft_mosaic
number of PFT mosaic
real(rp), dimension(:,:), allocatable, public real_dlat
delta latitude
integer, dimension(:,:,:), allocatable, public landuse_index_pft
index of PFT for each mosaic
module CONSTANT
Definition: scale_const.F90:14
module Convert LandUseIndex
module profiler
Definition: scale_prof.F90:10
integer, public ie
end point of inner domain: x, local
subroutine, public cnvlanduse_setup
Setup.
real(rp), public const_eps
small number
Definition: scale_const.F90:36
real(rp), dimension(:,:), allocatable, public landuse_frac_lake
lake fraction
logical, public io_lnml
output log or not? (for namelist, this process)
Definition: scale_stdio.F90:60
logical, public cnvlanduse_uselu100m
module PRECISION
logical, public cnvlanduse_usejibis
integer, public landuse_pft_nmax
number of plant functional type(PFT)
integer, public io_fid_conf
Config file ID.
Definition: scale_stdio.F90:55
integer, public io_fid_log
Log file ID.
Definition: scale_stdio.F90:56
integer, parameter, public rp
logical, public cnvlanduse_donothing
real(rp), dimension(:,:), allocatable, public real_lonx
longitude at staggered point (uy) [rad,0-2pi]
real(rp), dimension(:,:), allocatable, public landuse_frac_land
land fraction
real(rp), dimension(:,:), allocatable, public real_laty
latitude at staggered point (xv) [rad,-pi,pi]
integer, public ja
of y whole cells (local, with HALO)