From 934e977fa31b048ae4d00277fa0eb4adb707ee9b Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Tue, 30 May 2023 12:48:43 -0400 Subject: [PATCH 01/16] clean up make_bcs --- .../GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs | 2 +- .../Utils/Raster/makebcs/make_bcs_questionary.py | 4 ++-- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 10 +++++++++- 3 files changed, 12 insertions(+), 4 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs index 709bb037d..662206d74 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs @@ -1403,7 +1403,7 @@ if( $grid == ease ) then set BCNAME = ${EASEVERSION}_${HRCODE} set BCDIR = $EXPDIR/$OUTDIR/$BCNAME.scratch set BCJOB = $BCDIR/$BCNAME.j - set THISGRID = ${EASEVERSION}-${HRCODE} +# set THISGRID = ${EASEVERSION}-${HRCODE} set nfiles = `find $EXPDIR -maxdepth 5 -name ${BCNAME}".j" | wc -l` if( $nfiles >= 1 ) then diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py index 10b9b6f3d..c629a1be8 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py @@ -231,8 +231,8 @@ def ask_questions(default_grid="Cubed-Sphere"): "name": "EASEv1", "message": "Select EASEv1 grid resolution: \n ", "choices": [ \ - "M01 -- 1km $34668x14688$", \ - "M03 -- 3km $11556x4896$", \ +# "M01 -- 1km $34668x14688$", \ +# "M03 -- 3km $11556x4896$", \ "M09 -- 9km $3852x1632$", \ "M25 -- 25km $1383x586$", \ "M36 -- 36km $963x408$"], diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 6f6f4efa6..cee3368de 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -136,12 +136,20 @@ PROGRAM mkEASETilesParam gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid + if (index(EASELabel, 'v1') /=0) then + print*, " EASEv1_M03 is not implemented" + stop + endif regrid = .true. NC = 21600 NR = 10800 NT = 500000000 endif - if (index(EASELabel,'M03') /=0) then ! EASE 1 km grid + if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid + if (index(EASELabel, 'v1') /=0) then + print*, " EASEv1_M01 is not implemented" + stop + endif regrid = .true. NC = 43200 NR = 21600 From de99f877555febb585c7de3dfe965d8f95b332dc Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Wed, 7 Jun 2023 13:21:07 -0400 Subject: [PATCH 02/16] remove unnecessary check for M03 and M02 of EASEv1 --- .../Utils/Raster/makebcs/make_bcs_questionary.py | 4 ++-- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 8 -------- 2 files changed, 2 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py index c629a1be8..10b9b6f3d 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py @@ -231,8 +231,8 @@ def ask_questions(default_grid="Cubed-Sphere"): "name": "EASEv1", "message": "Select EASEv1 grid resolution: \n ", "choices": [ \ -# "M01 -- 1km $34668x14688$", \ -# "M03 -- 3km $11556x4896$", \ + "M01 -- 1km $34668x14688$", \ + "M03 -- 3km $11556x4896$", \ "M09 -- 9km $3852x1632$", \ "M25 -- 25km $1383x586$", \ "M36 -- 36km $963x408$"], diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index cee3368de..588be2d0c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -136,20 +136,12 @@ PROGRAM mkEASETilesParam gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid - if (index(EASELabel, 'v1') /=0) then - print*, " EASEv1_M03 is not implemented" - stop - endif regrid = .true. NC = 21600 NR = 10800 NT = 500000000 endif if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid - if (index(EASELabel, 'v1') /=0) then - print*, " EASEv1_M01 is not implemented" - stop - endif regrid = .true. NC = 43200 NR = 21600 From 9be02790caf070b3b7d403bb57ea61d4e684ce97 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Wed, 14 Jun 2023 18:12:03 -0400 Subject: [PATCH 03/16] clean ease grid options --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 331 +++++++++--------- 1 file changed, 171 insertions(+), 160 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 588be2d0c..840a04d2c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -28,7 +28,7 @@ PROGRAM mkEASETilesParam implicit none integer i,j,ig,jg,i0,iop,n,d1,d2,j1,j2,i1,i2,ix, jx,icount,pcount - integer :: NC = i_raster, NR = j_raster, NT = 16330000, ND = 10000, ND_raster = 10000 + integer :: NC, NR, NT, ND = 10000, ND_raster = 10000 integer, parameter :: nc_esa = 129600, nr_esa = 64800 @@ -62,7 +62,7 @@ PROGRAM mkEASETilesParam character(3) :: easegrid real :: clat, clon, r_ease, s_ease, da real :: fr_gcm - integer :: ind_col, ind_row, status, ncid, varid, nciv,nland_cells, DOM_INDX + integer :: ind_col, ind_row, status, ncid, varid, nciv,nland_cells, DOM_INDX, error_code, MPI_COMM_WORLD !REAL (kind=8), PARAMETER :: RADIUS=6378137.0,pi=3.14159265358979323846 character*100 :: veg_class (12) character*200 :: gfile,gtopo30 @@ -135,18 +135,26 @@ PROGRAM mkEASETilesParam write(nr_string, '(i0)') nr_ease gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) - if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid + if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid regrid = .true. - NC = 21600 - NR = 10800 - NT = 500000000 + NT = 1500000000 endif - if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid + if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid regrid = .true. - NC = 43200 - NR = 21600 NT = 1500000000 endif + if (index(EASELabel,'M09') /=0) then ! EASE 9 km grid + regrid = .true. + NT = 16330000 + endif + if (index(EASELabel,'M25') /=0) then ! EASE 25 km grid + regrid = .true. + NT = 16330000 + endif + if (index(EASELabel,'M36') /=0) then ! EASE 36 km grid + regrid = .true. + NT = 16330000 + endif allocate(land_id (1:NT)) allocate(water_id (1:NT)) @@ -178,7 +186,8 @@ PROGRAM mkEASETilesParam !else ! if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile !endif - + + ! Only MASK file we use with EASE grid in bcs package is GEOS5_10arcsec_mask, all other are not supported. if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then ! New ESA (Veg) + SRTM (catchments) based mask file ! is overlaid on the EASE grid @@ -315,155 +324,156 @@ PROGRAM mkEASETilesParam print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) else - - ! Old IGBP (Veg) + HYDRO1k (catchments) based mask will - ! Overlaid on EASE mask - ! ----------------------------------------------------- - - allocate(iaster (i_raster,j_raster)) - allocate(i2aster (i_raster,j_raster)) - allocate(veg (1:nc,1:nr)) - allocate(catid (1:nc,1:nr)) - allocate(catid_index (1:nc,1:nr)) - allocate(tileid_index(1:nc,1:nr)) - - dx = 360._8/nc - dy = 180._8/nr - d2r = MAPL_PI_R8/180._8 - !da = MAPL_radius*MAPL_radius*pi*pi*dx*dy/180./180./1000000. - - tileid_index = 0 - - ! Simple Biosphere 2 Model Legend - ! Value Class Name - ! (ftp://edcftp.cr.usgs.gov/pub/data/glcc/globe/latlon/sib22_0.leg) - ! the types vary 0-11 (array index minus 1) - - veg_class(1) = 'Ocean' - veg_class(2) = 'Broadleaf Evergreen Trees' - veg_class(3) = 'Broadleaf Deciduous Trees' - veg_class(4) = 'Broadleaf and Needleleaf Trees' - veg_class(5) = 'Needleleaf Evergreen Trees' - veg_class(6) = 'Needleleaf Deciduous Trees' - veg_class(7) = 'Short Vegetation/C4 Grassland' - veg_class(8) = 'Shrubs with Bare Soil' - veg_class(9) = 'Dwarf Trees and Shrubs' - veg_class(10) = 'Agriculture or C3 Grassland' - veg_class(11) = 'Water, Wetlands' - veg_class(12) = 'Ice/Snow' - - ! reading SiB2 land cover classification data - the origin of the - ! 2.5'x2.5' vegetation raster file is global 1min IGBP data - ! (ftp://edcftp.cr.usgs.gov/pub/data/glcc/globe/latlon/sib22_0.leg) - - open (10,file=trim(MAKE_BCS_INPUT_DIR)//'/land/veg/pft/v1/sib22.5_v2.0.dat', & - form='unformatted', & - action='read', convert='big_endian',status='old') - - READ(10)i2aster - - close (10,status='keep') - - if(regrid) then - call RegridRaster1 (i2aster,veg) - else - veg = i2aster - endif - - deallocate (i2aster) - - ! reading 2.5'x2.5' global raster file of Pfafstetter Catchment IDs - ! In this version, the dateline has been overlaid over the catchments those straddle - ! across. The numbers contain for - ! 1 global ocean catchment : Pfafstetter ID 0 - ! 36716 global land catchments : Pfafstetter IDs 1000-5999900 - ! 1 global inland water (lakes) catchment : Pfafstetter ID 6190000 - ! 1 global ice catchment : Pfafstetter ID 6200000 - - open (10,file= trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/global.cat_id.catch.DL', form='formatted', & - action='read', status='old')! - - do j=1,j_raster - read(10,*)(iaster(i,j),i=1,i_raster) - end do - - close (10,status='keep') - - if(regrid) then - call RegridRaster(iaster,catid) - else - catid = iaster - endif - - print *,'Read global.cat_id.catch.DL' - print *,'Min and Max of Pfafstetter IDs:', minval(catid),maxval(catid) - - ! reading the 2.5'x2.5' global raster file of tile indices for the - ! above Pfafstetter Catchments - ! 1 global ocean catchment : tile_index 36719 - ! 36716 global land catchments : tile_index 1-36716 - ! 1 global inland water (lakes) catchment : tile_index 36717 - ! 1 global ice catchment : tile_index 36718 - ! ------------------------------------------------------------ - - open (10,file=trim(MAKE_BCS_INPUT_DIR)//'/land/topo/' & - //'PfafstatterDL.rst', form='unformatted', & - action='read',convert='little_endian', status='old') - - do j=1,j_raster - read(10)(iaster(i,j),i=1,i_raster) - end do - - close (10,status='keep') - - if(regrid) then - call RegridRaster(iaster,catid_index) - else - catid_index = iaster - endif - - deallocate (iaster) - - print *,'Read PfafstatterDL.rst' - print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) - - ! While looping through the nc x nr grid (tile raster), this section counts # of - ! EASE grid cells that contain land, ice or water, seperately. - ! Each EASE grid cell is assigned with an ID = ind_row*ND + ind_col. - ! This is just the prelimiminery assessment in the process of assigning separate - ! tiles for land, water and ice fractions within the EASE Grid cell - ! The program checks each nc x nr pixels whether there is a EASE grid cell underneath, and counts - ! number of water, land and ice pixels as seen on veg raster. - ! ----------------------------------------------------------------------------------------------- - - - do i = 1 ,nc - - clon = -180. + float(i-1)*dx + dx/2. - - do j =nr ,1 ,-1 - - clat = -90. + float(j-1)*dy + dy/2. - call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) - - ind_col = nint(r_ease) + 1 - ind_row = nint(s_ease) + 1 - - if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_ease)) then - l= ind_row*ND + ind_col - - if(veg(i,j)==LakeType) then - water_id(l) = 1 - else if(veg(i,j)==IceType) then - ice_id (l) = 1 - else - land_id (l) = 1 - endif - endif - end do - end do - - endif + print *,'Only mask file we support with EASE grid is GEOS5_10arcsec_mask' + call MPI_Abort(MPI_COMM_WORLD,error_code,status) +! ! Old IGBP (Veg) + HYDRO1k (catchments) based mask will +! ! Overlaid on EASE mask +! ! ----------------------------------------------------- +! +! allocate(iaster (i_raster,j_raster)) +! allocate(i2aster (i_raster,j_raster)) +! allocate(veg (1:nc,1:nr)) +! allocate(catid (1:nc,1:nr)) +! allocate(catid_index (1:nc,1:nr)) +! allocate(tileid_index(1:nc,1:nr)) +! +! dx = 360._8/nc +! dy = 180._8/nr +! d2r = MAPL_PI_R8/180._8 +! !da = MAPL_radius*MAPL_radius*pi*pi*dx*dy/180./180./1000000. +! +! tileid_index = 0 +! +! ! Simple Biosphere 2 Model Legend +! ! Value Class Name +! ! (ftp://edcftp.cr.usgs.gov/pub/data/glcc/globe/latlon/sib22_0.leg) +! ! the types vary 0-11 (array index minus 1) +! +! veg_class(1) = 'Ocean' +! veg_class(2) = 'Broadleaf Evergreen Trees' +! veg_class(3) = 'Broadleaf Deciduous Trees' +! veg_class(4) = 'Broadleaf and Needleleaf Trees' +! veg_class(5) = 'Needleleaf Evergreen Trees' +! veg_class(6) = 'Needleleaf Deciduous Trees' +! veg_class(7) = 'Short Vegetation/C4 Grassland' +! veg_class(8) = 'Shrubs with Bare Soil' +! veg_class(9) = 'Dwarf Trees and Shrubs' +! veg_class(10) = 'Agriculture or C3 Grassland' +! veg_class(11) = 'Water, Wetlands' +! veg_class(12) = 'Ice/Snow' +! +! ! reading SiB2 land cover classification data - the origin of the +! ! 2.5'x2.5' vegetation raster file is global 1min IGBP data +! ! (ftp://edcftp.cr.usgs.gov/pub/data/glcc/globe/latlon/sib22_0.leg) +! +! open (10,file=trim(MAKE_BCS_INPUT_DIR)//'/land/veg/pft/v1/sib22.5_v2.0.dat', & +! form='unformatted', & +! action='read', convert='big_endian',status='old') +! +! READ(10)i2aster +! +! close (10,status='keep') +! +! if(regrid) then +! call RegridRaster1 (i2aster,veg) +! else +! veg = i2aster +! endif +! +! deallocate (i2aster) +! +! ! reading 2.5'x2.5' global raster file of Pfafstetter Catchment IDs +! ! In this version, the dateline has been overlaid over the catchments those straddle +! ! across. The numbers contain for +! ! 1 global ocean catchment : Pfafstetter ID 0 +! ! 36716 global land catchments : Pfafstetter IDs 1000-5999900 +! ! 1 global inland water (lakes) catchment : Pfafstetter ID 6190000 +! ! 1 global ice catchment : Pfafstetter ID 6200000 +! +! open (10,file= trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/global.cat_id.catch.DL', form='formatted', & +! action='read', status='old')! +! +! do j=1,j_raster +! read(10,*)(iaster(i,j),i=1,i_raster) +! end do +! +! close (10,status='keep') +! +! if(regrid) then +! call RegridRaster(iaster,catid) +! else +! catid = iaster +! endif +! +! print *,'Read global.cat_id.catch.DL' +! print *,'Min and Max of Pfafstetter IDs:', minval(catid),maxval(catid) +! +! ! reading the 2.5'x2.5' global raster file of tile indices for the +! ! above Pfafstetter Catchments +! ! 1 global ocean catchment : tile_index 36719 +! ! 36716 global land catchments : tile_index 1-36716 +! ! 1 global inland water (lakes) catchment : tile_index 36717 +! ! 1 global ice catchment : tile_index 36718 +! ! ------------------------------------------------------------ +! +! open (10,file=trim(MAKE_BCS_INPUT_DIR)//'/land/topo/' & +! //'PfafstatterDL.rst', form='unformatted', & +! action='read',convert='little_endian', status='old') +! +! do j=1,j_raster +! read(10)(iaster(i,j),i=1,i_raster) +! end do +! +! close (10,status='keep') +! +! if(regrid) then +! call RegridRaster(iaster,catid_index) +! else +! catid_index = iaster +! endif +! +! deallocate (iaster) +! +! print *,'Read PfafstatterDL.rst' +! print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) +! +! ! While looping through the nc x nr grid (tile raster), this section counts # of +! ! EASE grid cells that contain land, ice or water, seperately. +! ! Each EASE grid cell is assigned with an ID = ind_row*ND + ind_col. +! ! This is just the prelimiminery assessment in the process of assigning separate +! ! tiles for land, water and ice fractions within the EASE Grid cell +! ! The program checks each nc x nr pixels whether there is a EASE grid cell underneath, and counts +! ! number of water, land and ice pixels as seen on veg raster. +! ! ----------------------------------------------------------------------------------------------- +! +! +! do i = 1 ,nc +! +! clon = -180. + float(i-1)*dx + dx/2. +! +! do j =nr ,1 ,-1 +! +! clat = -90. + float(j-1)*dy + dy/2. +! call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) +! +! ind_col = nint(r_ease) + 1 +! ind_row = nint(s_ease) + 1 +! +! if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_ease)) then +! l= ind_row*ND + ind_col +! +! if(veg(i,j)==LakeType) then +! water_id(l) = 1 +! else if(veg(i,j)==IceType) then +! ice_id (l) = 1 +! else +! land_id (l) = 1 +! endif +! endif +! end do +! end do +! + endif ! (GEOS5_10arcsec_mask) ! Reading SRTM elevation data - to be consistent with AGCM ! -------------------------------------------------------- @@ -668,7 +678,8 @@ PROGRAM mkEASETilesParam if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then pfaf = cindex else - pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) + print *,'Only mask file we support with EASE grid is GEOS5_10arcsec_mask' + ! pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) endif if ((l > l_index).and.(l <= w_index)) typ =19 From de78331a34ad8f7f8b9410b865543e46bfe63b71 Mon Sep 17 00:00:00 2001 From: biljanaorescanin Date: Thu, 15 Jun 2023 09:11:31 -0400 Subject: [PATCH 04/16] missing stop + remove M01 until ready for users --- .../Utils/Raster/makebcs/make_bcs_questionary.py | 4 ++-- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py index 10b9b6f3d..6b6954274 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py @@ -231,7 +231,7 @@ def ask_questions(default_grid="Cubed-Sphere"): "name": "EASEv1", "message": "Select EASEv1 grid resolution: \n ", "choices": [ \ - "M01 -- 1km $34668x14688$", \ + #"M01 -- 1km $34668x14688$", \ "M03 -- 3km $11556x4896$", \ "M09 -- 9km $3852x1632$", \ "M25 -- 25km $1383x586$", \ @@ -243,7 +243,7 @@ def ask_questions(default_grid="Cubed-Sphere"): "name": "EASEv2", "message": "Select EASEv2 grid resolution: \n ", "choices": [ \ - "M01 -- 1km $34704x14616$", \ + #"M01 -- 1km $34704x14616$", \ "M03 -- 3km $11568x4872$", \ "M09 -- 9km $3856x1624$", \ "M25 -- 25km $1388x584$", \ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 840a04d2c..6d5c52b6f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -679,6 +679,7 @@ PROGRAM mkEASETilesParam pfaf = cindex else print *,'Only mask file we support with EASE grid is GEOS5_10arcsec_mask' + call MPI_Abort(MPI_COMM_WORLD,error_code,status) ! pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) endif From 2e5cbb4157dfbf0e149ef1d8a397006c40ce451f Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 17 Jun 2023 14:12:39 -0400 Subject: [PATCH 05/16] major cleanup of mkEASETilesParam.F90: - rewrote some code blocks for clarity and efficiency - removed obsolete variable declarations - removed obsolete use statements - added "only:" qualifier to some use statement - commented out obsolete code blocks - renamed some variables for clarity - removed repetition of identical operations - added comments - white-space changes for improved readability --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 796 +++++++++++------- 1 file changed, 492 insertions(+), 304 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 6d5c52b6f..c985952fb 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -3,7 +3,7 @@ PROGRAM mkEASETilesParam - ! This program constructs land and lake tiles for EASE grid tile spaces such as + ! This program constructs land, landice, and lake tiles for EASE grid tile spaces such as ! those used for the SMAP Level-4 products and other offline projects ! This program resulted from the merger and cleanup of mkSMAPTilesPara.F90 @@ -14,77 +14,100 @@ PROGRAM mkEASETilesParam ! ! - wjiang + reichle, 21 Sep 2022 - - use EASE_conv + ! additional major cleanup (reichle, 15 Jun 2023) + ! - rewrote some code blocks for clarity and efficiency + ! - removed obsolete variable declarations + ! - removed obsolete use statements + ! - added "only:" qualifier to some use statement + ! - commented out obsolete code blocks + ! - renamed some variables for clarity + ! - removed repetition of identical operations + ! - added comments + ! - white-space changes for improved readability + + use EASE_conv, only : EASE_extent, EASE_convert, EASE_inverse use rmTinyCatchParaMod, only : i_raster, j_raster, SRTM_maxcat - use rmTinyCatchParaMod, only : RegridRaster, RegridRaster1, RegridRasterReal + use rmTinyCatchParaMod, only : RegridRasterReal use rmTinyCatchParaMod, only : MAKE_BCS_INPUT_DIR - use process_hres_data + use process_hres_data, only : histogram use MAPL_SortMod use MAPL_ConstantsMod use MAPL_ExceptionHandling - use LogRectRasterizeMod use netcdf implicit none - integer i,j,ig,jg,i0,iop,n,d1,d2,j1,j2,i1,i2,ix, jx,icount,pcount - integer :: NC, NR, NT, ND = 10000, ND_raster = 10000 + + integer :: i, j, ig, jg, n + integer :: NC, NR, N_ease_grid_cells, NDND, ND_raster integer, parameter :: nc_esa = 129600, nr_esa = 64800 ! For regridding - integer, allocatable, target, dimension (:,:) & - :: geos_msk - REAL, allocatable, DIMENSION (:) :: loc_val - INTEGER, ALLOCATABLE, DIMENSION (:) :: density, loc_int - logical, dimension (:), allocatable :: unq_mask - integer, pointer , dimension (:,:) :: subset - integer, pointer , dimension (:) :: subset1 - real, pointer , dimension (:) :: subset2 - integer :: dx_esa, dy_esa, NBINS, NPLUS - - integer*8, allocatable, dimension (:) :: SRTM_catid - real(kind=8),allocatable, dimension(:):: SRTM_catid_r8 - - integer,allocatable, dimension (:,:), target :: tileid_index,catid_index - integer,allocatable, dimension (:,:) :: catid, iaster - integer,allocatable, dimension (:) :: land_id,water_id,ice_id - integer,allocatable, dimension (:) :: my_land, all_id - real, allocatable, dimension (:) :: ease_grid_area,tile_area,SRTM_CatchArea - integer*1,allocatable, dimension (:,:) :: veg, i2aster - real*4, dimension (:,:), allocatable :: q0,raster - REAL, dimension (:), allocatable :: tile_ele, tile_area_land - - INTEGER*8 :: PFAF_CODE - integer l,imn,imx,jmn,jmx,mval,l_index,i_index,w_index,typ,pfaf,cindex - integer :: LakeType, IceType, OceanType - character(3) :: easegrid - real :: clat, clon, r_ease, s_ease, da - real :: fr_gcm - integer :: ind_col, ind_row, status, ncid, varid, nciv,nland_cells, DOM_INDX, error_code, MPI_COMM_WORLD - !REAL (kind=8), PARAMETER :: RADIUS=6378137.0,pi=3.14159265358979323846 - character*100 :: veg_class (12) - character*200 :: gfile,gtopo30 - integer :: nc_ease,nr_ease, N_args, command_argument_count - REAL :: dx,dy,d2r,lats,mnx,mxx,mny,mxy,sum1,sum2,jgv, VDUM,pix_area - character(40) :: arg, EASElabel_ + integer, allocatable, dimension(:,:), target :: geos_msk + + REAL, allocatable, DIMENSION(:) :: loc_val + INTEGER, ALLOCATABLE, DIMENSION(:) :: density, loc_int + logical, allocatable, dimension(:) :: unq_mask + integer, dimension(:,:), pointer :: subset + + integer :: dx_esa, dy_esa, NBINS, NPLUS + + integer*8, allocatable, dimension(:) :: SRTM_catid + real(kind=8),allocatable, dimension(:) :: SRTM_catid_r8 + + integer, allocatable, dimension(:,:), target :: tileid_index, catid_index + + ! integer, allocatable, dimension(:,:) :: catid, iaster + + integer, allocatable, dimension(:) :: land_id, water_id, ice_id + integer, allocatable, dimension(:) :: my_land, all_id + real, allocatable, dimension(:) :: ease_grid_area, tile_area !, SRTM_CatchArea + integer*1, allocatable, dimension(:,:) :: veg !, i2aster + real*4, allocatable, dimension(:,:) :: q0, raster + REAL, allocatable, dimension(:) :: tile_ele + + !INTEGER*8 :: PFAF_CODE + + integer :: l, l_index, i_index, w_index, typ, pfaf, cindex + integer :: LandType, LakeType, IceType, OceanType + real :: clat, clon, r_ease, s_ease + real :: fr_gcm + integer :: ind_col, ind_row, status, ncid, varid, nciv !, DOM_INDX + integer :: n_land, n_lake, n_landice, n_landlakelandice + + ! REAL (kind=8), PARAMETER :: RADIUS=6378137.0,pi=3.14159265358979323846 + + !character*100 :: veg_class(12) + + character*200 :: gfile, gtopo30 + integer :: nc_ease, nr_ease, N_args, command_argument_count + REAL :: dx, dy, d2r, lats, mnx, mxx, mny, mxy, jgv, pix_area ! VDUM + real :: dx_ease, mean_land_ele + character(40) :: arg, EASElabel_ + character(len=:), allocatable :: EASElabel - character*200 :: tmpstring, tmpstring1, tmpstring2 - logical :: regrid = .false. + + !character*200 :: tmpstring, tmpstring1, tmpstring2 + + logical :: regrid = .false. character*128 :: MaskFile - logical :: pfaf_til = .false. + + !logical :: pfaf_til = .false. + character*1 :: PF - character(len=6) :: EASE_Version + + !character(len=6) :: EASE_Version + character(len=10) :: nc_string, nr_string character(len=128) :: usage1, usage2 character(len=128) :: Iam = "mkEASETilesParam" - call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) - ! -------------------------------------------------------------------------------------- - usage1 = 'USAGE : bin/mkEASETilesParam.x -ease_label EASELabel ' + call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) + + usage1 = 'USAGE : bin/mkEASETilesParam.x -ease_label EASELabel ' usage2 = ' where EASELabel = *EASEv[x]_M[yy]*, x={1,2}, yy={01,03,09,25,36}' N_args = command_argument_count() @@ -120,55 +143,62 @@ PROGRAM mkEASETilesParam end do - ! WY note, remove this verification. There can be all combination - ! if (MGRID /= 'M25' .and. EASE_version == 'EASEv1') then - ! stop ("EASEv1 only supports M25") - ! endif - - ! Setting EASE Grid specifications + ! Get EASE Grid specifications ! -------------------------------- - + EASElabel = trim(EASELabel_) - - call ease_extent(EASELabel, nc_ease, nr_ease ) + + call ease_extent( EASELabel, nc_ease, nr_ease ) + write(nc_string, '(i0)') nc_ease write(nr_string, '(i0)') nr_ease + gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) - - if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid - regrid = .true. - NT = 1500000000 - endif - if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid - regrid = .true. - NT = 1500000000 - endif - if (index(EASELabel,'M09') /=0) then ! EASE 9 km grid - regrid = .true. - NT = 16330000 - endif - if (index(EASELabel,'M25') /=0) then ! EASE 25 km grid - regrid = .true. - NT = 16330000 - endif - if (index(EASELabel,'M36') /=0) then ! EASE 36 km grid - regrid = .true. - NT = 16330000 - endif - - allocate(land_id (1:NT)) - allocate(water_id (1:NT)) - allocate(ice_id (1:NT)) + +! ! get appropriate number for length of land_id, water_id, ice_id vectors +! +! if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid +! N_ease_grid_cells = 1500000000 +! else if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid +! N_ease_grid_cells = 1500000000 +! else if (index(EASELabel,'M09') /=0) then ! EASE 9 km grid +! N_ease_grid_cells = 16330000 +! else if (index(EASELabel,'M25') /=0) then ! EASE 25 km grid +! N_ease_grid_cells = 16330000 +! else if (index(EASELabel,'M36') /=0) then ! EASE 36 km grid +! N_ease_grid_cells = 16330000 +! else +! print *,"ERROR: Unknown EASELabel, stopping." +! stop +! endif + + N_ease_grid_cells = nc_ease * nr_ease + + allocate(land_id (1:N_ease_grid_cells)) + allocate(water_id(1:N_ease_grid_cells)) + allocate(ice_id (1:N_ease_grid_cells)) land_id = 0 water_id = 0 ice_id = 0 + + ! define tile types used for processing here (possibly from ESA mask) + OceanType = 0 - IceType = 11 - LakeType = 10 + IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 + LakeType = 10 ! lake type used for processing here; in GEOS, lake tiles are type= 19 + LandType = 1 ! land type used for processing here; in GEOS, land tiles are type=100 - ND = 10*10**(nint(log10(1.*nr_ease))) + ! Step size for conversion of 2-dim indexing into 1-dim indexing: l = ind_row*NDND + ind_col + ! If conversion was simply to get from 2-dim to 1-dim indexing for EASE grid cells, + ! NDND=nc_ease would suffice + ! Here, NDND is computed as a power of 10, leaving room for "empty" indices that may or may not be needed later. + ! - reichle, 15 Jun 2023 + + !NDND = 10*10**(nint(log10(1.*nr_ease))) + NDND = nc_ease + ! Check for the 10 arc-sec MaskFile ! ----------------------------------- @@ -186,63 +216,87 @@ PROGRAM mkEASETilesParam !else ! if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile !endif - - ! Only MASK file we use with EASE grid in bcs package is GEOS5_10arcsec_mask, all other are not supported. + if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then + ! New ESA (Veg) + SRTM (catchments) based mask file ! is overlaid on the EASE grid ! ------------------------------------------------- - nc = 43200 ! Number of rows in raster file + ! ESA/SRTM mask file is 10-arcsec resolution; coarsen here to 30-arcsec raster grid: + + nc = 43200 ! Number of rows in raster file nr = 21600 ! Number of columns in raster file regrid = .true. - dx_esa = nc_esa / nc ! x-dimension (or # of ESA columns within the raster grid cell) - dy_esa = nr_esa / nr ! y-dimension (or # of ESA rows within the raster grid cell) + + dx_esa = nc_esa / nc ! =3 (# of ESA columns within a raster grid cell in x-dim) + dy_esa = nr_esa / nr ! =3 (# of ESA rows within a raster grid cell in y-dim) allocate(tileid_index(1:nc,1:nr)) - allocate(SRTM_catid (1:SRTM_maxcat+2)) - allocate(SRTM_catid_r8(1:SRTM_maxcat+2), source = 0.d0) + allocate(SRTM_catid (1:SRTM_maxcat+2) ) + allocate(SRTM_catid_r8 (1:SRTM_maxcat+2), source = 0.d0) allocate(catid_index (1:nc,1:nr)) allocate(veg (1:nc,1:nr)) allocate(geos_msk (1:nc_esa,1:dy_esa)) - allocate(SRTM_CatchArea (1:SRTM_maxcat)) - - OPEN (10, FILE = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat', & - FORM = 'FORMATTED',STATUS='OLD',ACTION='READ') - - READ (10,*) I - DO N = 1, I - READ (10, '(i8,i15,4(1x,f9.4),1x,e10.3,4(1x,e9.3),I8,6(1x,f9.4))') & - DOM_INDX,PFAF_CODE,VDUM,VDUM,VDUM,VDUM,VDUM, & - SRTM_CatchArea (N) - END DO - CLOSE (10, STATUS='KEEP') - - dx = 360._8/nc - dy = 180._8/nr - d2r = MAPL_PI_R8/180._8 + + ! the following block is not needed (perhaps for routing) + +! allocate(SRTM_CatchArea(1:SRTM_maxcat)) +! +! OPEN (10, FILE = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat', & +! FORM = 'FORMATTED',STATUS='OLD',ACTION='READ') +! +! READ (10,*) I +! DO N = 1, I +! READ (10, '(i8,i15,4(1x,f9.4),1x,e10.3,4(1x,e9.3),I8,6(1x,f9.4))') & +! DOM_INDX,PFAF_CODE,VDUM,VDUM,VDUM,VDUM,VDUM, & +! SRTM_CatchArea (N) +! END DO +! CLOSE (10, STATUS='KEEP') + + dx = 360._8/nc ! raster grid spacing (dlon) + dy = 180._8/nr ! raster grid spacing (dlat) + + d2r = MAPL_PI_R8/180._8 ! degree-to-radians conversion factor -- d2r declared as REAL ?!?!?! + !da = MAPL_radius*MAPL_radius*pi*pi*dx*dy/180./180./1000000. tileid_index = 0 catid_index = 0 - veg = 0 + veg = OceanType ! initialize to ocean (type=0) + ! read list of Pfafstetter catchment IDs ('PfafID') from 10-arcsec mask file into variable 'SRTM_catid[_r8]' + ! [vector of length SRTM_maxcat] + status = NF90_OPEN (trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/GEOS5_10arcsec_mask.nc', NF90_NOWRITE, ncid) status = nf90_inq_varid(ncid, name='PfafID', varid=varid) - status = nf90_get_var(ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/)) + status = nf90_get_var(ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/)) if(status /=0) then PRINT *, NF90_STRERROR(STATUS) print *, 'Problem with NF90_OPEN',trim(MaskFile) endif - SRTM_catid = int8(SRTM_catid_r8) - SRTM_catid (SRTM_maxcat + 1) = 190000000 - SRTM_catid (SRTM_maxcat + 2) = 200000000 - i1 = 0 ! count # of 30-arcsec pixels + SRTM_catid = int8(SRTM_catid_r8) ! convert data to integer*8 -- contains 12-digit Pfaf code + SRTM_catid (SRTM_maxcat + 1) = 190000000 ! append ID for Lake type + SRTM_catid (SRTM_maxcat + 2) = 200000000 ! append ID for Landice type + + ! ----------------------------------------------- + ! + ! first loop through 30-arcsec raster grid cells + ! - aggregate mask from 10-arcsec resolution to 30-arcsec raster grid by determining dominant tile type + ! - for each EASE grid cell, determine tile types based on aggregated (30-arcsec) mask + + !i1 = 0 ! count # of 30-arcsec pixels - NOT USED + status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) + do j=1,nr + + clat = -90. + float(j-1)*dy + dy/2. ! center lat of raster grid cell (*,j) + + ! read slice of variable 'CatchIndex' from 10-arcsec mask file into variable 'geos_msk' + ! [2d-array: 129600-by-3] - clat = -90. + float(j-1)*dy + dy/2. status = NF90_GET_VAR (ncid, varid, geos_msk, (/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/)) ! Read 10-arcsec rows that lie within the raster row 'j' !status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' @@ -253,79 +307,105 @@ PROGRAM mkEASETilesParam do i = 1,nc - clon = -180. + float(i-1)*dx + dx/2. + clon = -180. + float(i-1)*dx + dx/2. ! center lon of raster grid cell (i,*) + + ! extract [3-by-3] subset of ESA/SRTM 10-arcsec mask file that corresponds to 30-arcsec raster grid cell (i,j) if (associated (subset)) NULLIFY (subset) + subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains ESA pixels that lie within the raster grid cell at i,j + if(maxval (subset) > SRTM_maxcat) then - where (subset == 190000000) subset = SRTM_maxcat + 1 - where (subset == 200000000) subset = SRTM_maxcat + 2 + where (subset == 190000000) subset = SRTM_maxcat + 1 ! Lake type (convert ID from 190000000 to SRTM_maxcat+1) + where (subset == 200000000) subset = SRTM_maxcat + 2 ! Landice type (convert ID from 200000000 to SRTM_maxcat+2) endif - - if (maxval (subset) > 0) then ! check whether there are Non-ocean ESA pixels - ! catID of the middle pixel - veg (i,j) = 1 ! veg is set to land + if (maxval(subset) > 0) then ! check if there are Non-ocean ESA pixels - NPLUS = count(subset >= 1 .and. subset <= SRTM_maxcat + 2) ! Count non-ocean ESA pixels within - allocate (loc_int (1:NPLUS)) - allocate (unq_mask(1:NPLUS)) - loc_int = pack(subset,mask = (subset >= 1 .and. subset <= SRTM_maxcat + 2)) ! loc_int contains catch_indices of non-ocean ESA pixels - call MAPL_Sort (loc_int) + ! raster grid cell has at least one 10-arcsec pixel that is land or landice or lake + ! + ! now find out the *dominant* Pfafstetter ID (could be Lake or Landice) within the raster grid cell + + NPLUS = count( subset>=1 .and. subset<=SRTM_maxcat+2 ) ! Count non-ocean ESA pixels within raster grid cell + allocate(loc_int (1:NPLUS)) + allocate(unq_mask(1:NPLUS)) + loc_int = pack( subset, mask = ( subset>= 1 .and. subset<=SRTM_maxcat+2) ) ! loc_int contains catch_indices of non-ocean ESA pixels + call MAPL_Sort(loc_int) unq_mask = .true. - do n = 2,NPLUS - unq_mask(n) = .not.(loc_int(n) == loc_int(n-1)) ! count number of unique numbers in loc_int for binning + do n=2,NPLUS + unq_mask(n) = .not. (loc_int(n)==loc_int(n-1)) ! count number of unique numbers in loc_int for binning end do NBINS = count(unq_mask) - + if (NBINS > 1) then - allocate(loc_val (1:NBINS)) - allocate(density (1:NBINS)) + allocate(loc_val(1:NBINS)) + allocate(density(1:NBINS)) loc_val = 1.*pack(loc_int,mask =unq_mask) ! loc_val contains available non-ocean catch_indices within the i,j grid cell, ! Those numbers will be used as bin values - call histogram (dx_esa*dy_esa, NBINS, density, loc_val, real(subset)) ! density is the pixel count for each bin value - catid_index (i,j) = loc_val (maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j - deallocate (loc_val, density) + call histogram( dx_esa*dy_esa, NBINS, density, loc_val, real(subset) ) ! density is the pixel count for each bin value + catid_index(i,j) = loc_val(maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j + deallocate(loc_val, density) else - catid_index (i,j) = loc_int (1) + catid_index(i,j) = loc_int (1) endif - deallocate (loc_int, unq_mask) + deallocate(loc_int, unq_mask) - if(catid_index (i,j) == SRTM_maxcat + 1) veg (i,j) = LakeType - if(catid_index (i,j) == SRTM_maxcat + 2) veg (i,j) = IceType - if((catid_index(i,j) >= 1).and.(catid_index (i,j) <= SRTM_maxcat)) i1 = i1 + 1 + ! now catid_index(i,j) = dominant PfafID (or LakeID or LandiceID) in raster grid cell (i,j) + + if (catid_index (i,j) <= SRTM_maxcat ) veg(i,j) = LandType + if (catid_index (i,j) == SRTM_maxcat + 1) veg(i,j) = LakeType + if (catid_index (i,j) == SRTM_maxcat + 2) veg(i,j) = IceType + !if( (catid_index(i,j)>= 1) .and. (catid_index(i,j)<=SRTM_maxcat) ) i1 = i1 + 1 ! i1 = counter for raster grid cells with land as dominant type -- NOT USED? + + ! set land_id or water_id or ice_id of EASE grid cell to 1 accordingly + ! --> EASE grid cell may contain more than one tile type + ! count in if this is i,j pixel is a land, lake or ice within ind_col,ind_row EASE grid cell + ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) ind_col = nint(r_ease) + 1 - ind_row = nint(s_ease) + 1 + ind_row = nint(s_ease) + 1 ! can be negative or greater than nr_ease (lat near N/S pole) - if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_ease)) then - l= ind_row*ND + ind_col + if( (veg(i,j).ne.OceanType) .and. (ind_row.ge.1) .and. (ind_row.le.nr_ease) ) then - if(veg(i,j)==LakeType) then + ! raster grid cell has tile type other than ocean and is located within the lat band covered + ! by the EASE grid (approx 85S-85N): 1 <= ind_row <= nr_ease + + ! 1-dim indexing used here appears wasteful; why not ll=(ind_row-1)*nc_ease+ind_col, with NT=nc_ease*nr_ease ?? + + !l = ind_row*NDND + ind_col + + l = (ind_row-1)*NDND + ind_col ! 1-dim index for all EASE grid cells + + if( veg(i,j)==LakeType) then water_id(l) = 1 - else if(veg(i,j)==IceType) then + else if(veg(i,j)==IceType ) then ice_id (l) = 1 else land_id (l) = 1 endif endif endif - end do - enddo - + + end do ! i=1,nc + enddo ! j=1,nr + status = NF90_CLOSE (ncid) deallocate (geos_msk) - print *,'Read ', trim (MaskFile) - print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) + print *,'Done reading ', trim(MaskFile) + print *,'Min and Max of tile indices:', minval(catid_index), maxval(catid_index) else - print *,'Only mask file we support with EASE grid is GEOS5_10arcsec_mask' - call MPI_Abort(MPI_COMM_WORLD,error_code,status) + + print *,'MaskFile = ', trim(MaskFile) + print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' + stop + ! ! Old IGBP (Veg) + HYDRO1k (catchments) based mask will ! ! Overlaid on EASE mask ! ! ----------------------------------------------------- @@ -439,7 +519,7 @@ PROGRAM mkEASETilesParam ! ! ! While looping through the nc x nr grid (tile raster), this section counts # of ! ! EASE grid cells that contain land, ice or water, seperately. -! ! Each EASE grid cell is assigned with an ID = ind_row*ND + ind_col. +! ! Each EASE grid cell is assigned with an ID = ind_row*NDND + ind_col. ! ! This is just the prelimiminery assessment in the process of assigning separate ! ! tiles for land, water and ice fractions within the EASE Grid cell ! ! The program checks each nc x nr pixels whether there is a EASE grid cell underneath, and counts @@ -460,7 +540,7 @@ PROGRAM mkEASETilesParam ! ind_row = nint(s_ease) + 1 ! ! if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_ease)) then -! l= ind_row*ND + ind_col +! l= ind_row*NDND + ind_col ! ! if(veg(i,j)==LakeType) then ! water_id(l) = 1 @@ -472,14 +552,16 @@ PROGRAM mkEASETilesParam ! endif ! end do ! end do -! - endif ! (GEOS5_10arcsec_mask) + + + endif ! (GEOS5_10arcsec_mask) - ! Reading SRTM elevation data - to be consistent with AGCM ! -------------------------------------------------------- + ! + ! Read SRTM elevation data - to be consistent with AGCM - allocate(raster (i_raster,j_raster)) - allocate(q0(nc,nr)) + allocate(raster(i_raster,j_raster)) ! 2.5-min raster grid ( 8640-by-4320 ) + allocate(q0 (nc, nr) ) ! 30-arcsec raster grid (43200-by-21600) gtopo30 = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/srtm30_withKMS_2.5x2.5min.data' @@ -487,39 +569,54 @@ PROGRAM mkEASETilesParam read (10) raster close (10,status='keep') + ! remap SRTM elevation data from 2.5-min to 30-arcsec raster grid + if(regrid) then call RegridRasterReal(raster,q0) else - q0 = raster + q0 = raster endif deallocate (raster) - print *,'# of Land pixels in EASE: ',sum (land_id) - print *,'# of water pixels in EASE: ',sum (water_id) - print *,'# of ice pixels in EASE: ',sum (ice_id) + ! --------------------------------------------------------- + ! + ! determine number of land, lake, and landice tiles + + n_land = sum(land_id) ! number of land tiles + n_lake = sum(water_id) ! number of lake tiles + n_landice = sum(ice_id) ! number of landice tiles + + n_landlakelandice = n_land + n_lake + n_landice + + print *,'# of Land tiles: ', n_land + print *,'# of Lake tiles: ', n_lake + print *,'# of Landice tiles: ', n_landice + + l_index = 0 ! start index for land EASE cells + w_index = n_land ! start index for lake EASE cells + i_index = w_index + n_lake ! start index for landice EASE cells - l_index=0 - w_index=sum (land_id) - i_index=sum (land_id) + sum (water_id) - nland_cells = w_index + allocate(ease_grid_area(1:N_ease_grid_cells)) + allocate(tile_area (1:n_landlakelandice)) + allocate(my_land (1:n_landlakelandice)) + allocate(all_id (1:n_landlakelandice)) - allocate(tile_area (1:i_index + sum (ice_id))) - allocate(ease_grid_area (1:NT)) - allocate(tile_ele (1:w_index)) - allocate(tile_area_land(1:w_index)) - allocate(my_land (1:i_index + sum (ice_id))) - allocate(all_id (1:i_index + sum (ice_id))) + allocate(tile_ele (1:n_land)) - land_id = 0 - water_id= 0 - ice_id = 0 + ! =========================================================================== + ! + ! prepare for second loop through raster grid cells + + land_id = 0 + water_id = 0 + ice_id = 0 + + my_land = 0 + all_id = 0 - my_land = 0 - all_id = 0 ease_grid_area = 0. - tile_area_land = 0. tile_ele = 0. tile_area = 0. @@ -537,75 +634,107 @@ PROGRAM mkEASETilesParam ! is derived in the below loop ND_raster = 10*10**(nint(log10(1.*NR))) - i2 = 1 + + !i2 = 1 -- NOT USED do i = 1 ,nc - clon = -180. + float(i-1)*dx + dx/2. + clon = -180. + float(i-1)*dx + dx/2. ! center lon of raster grid cell (*,j) do j =nr ,1 ,-1 - lats = -90._8 + (j - 0.5_8)*dy - clat = -90. + float(j-1)*dy + dy/2. - call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) + + lats = -90._8 + (j - 0.5_8)*dy ! center lat of raster grid cell (*,j) -- lats declared REAL ?!?!?! same as clat ?!?!?! + clat = -90. + float(j-1)*dy + dy/2. ! center lat of raster grid cell (*,j) + + ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) + + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) ind_col = nint(r_ease) + 1 - ind_row = nint(s_ease) + 1 - - l= ind_row*ND + ind_col - pix_area =(sin(d2r*(lats+0.5*dy)) -sin(d2r*(lats-0.5*dy)))*(dx*d2r) - - if((ind_row.ge.1).and.(veg(i,j).ge.1).and.(ind_row.le.nr_ease)) then - - if(veg(i,j)==LakeType) then - if(water_id(l)==0) then - w_index = w_index + 1 - water_id(l) = w_index - tileid_index(i,j)= water_id(l) - else - tileid_index(i,j)= water_id(l) - endif - endif + ind_row = nint(s_ease) + 1 ! can be negative or greater than nr_ease (lat near N/S pole) + + if( (ind_row.ge.1) .and. (ind_row.le.nr_ease) ) then - if(veg(i,j)==IceType) then - if(ice_id(l)==0) then - i_index = i_index + 1 - ice_id (l) = i_index - tileid_index(i,j)= ice_id (l) !i_index - else - tileid_index(i,j)= ice_id (l) !i_index - endif - endif + ! raster grid cell (i,j) is located within the lat band covered by the (cylindrical) EASE grid (approx 85S-85N) - if(veg(i,j).lt.LakeType) then - if(land_id(l)==0) then - l_index = l_index + 1 - land_id (l) = l_index - tileid_index(i,j)= land_id (l) !1-l_index - else - tileid_index(i,j)= land_id (l) !1-l_index - endif - - endif - tile_area(tileid_index(i,j))= tile_area(tileid_index(i,j)) + & - pix_area - my_land(tileid_index(i,j)) = l - all_id (tileid_index(i,j)) = j*ND_raster + i - endif - - if((ind_row.ge.1).and.(ind_row.le.nr_ease)) then - ease_grid_area(l) = ease_grid_area(l) + & - pix_area - endif + !l= ind_row*NDND + ind_col ! 1-dim index for EASE grid cells + + l = (ind_row-1)*NDND + ind_col ! 1-dim index for all EASE grid cells + + pix_area = ( sin(d2r*(lats+0.5*dy)) - sin(d2r*(lats-0.5*dy)) )*(dx*d2r) ! area of 30-arcsec raster grid cell + + if (veg(i,j).ne.OceanType) then + + ! set tile ID and compute tile area + + select case (veg(i,j)) + + case (LakeType) ! raster grid cell (i,j) is Lake + + if(water_id(l)==0) then + ! raster grid cell is first Lake type seen for this EASE grid cell, tile ID has not yet been set + ! recall: tile IDs for Lake: [(n_land+1):(n_land+n_lake)], and w_index was initalized to n_land above + w_index = w_index + 1 + tileid_index(i,j)= w_index + else + ! no action needed (tile ID does not depend on number of contributing Lake raster grid cells) + endif + + case (IceType) ! raster grid cell (i,j) is Landice + + if(ice_id(l)==0) then + ! raster grid cell is first Landice type seen for this EASE grid cell, tile ID has not yet been set + ! recall: tile IDs for Landice: [(n_land+n_lake+1):n_landlakelandice], and i_index was initalized to n_land+n_lake above + i_index = i_index + 1 + tileid_index(i,j)= i_index + else + ! no action needed (tile ID does not depend on number of contributing Landice raster grid cells) + endif + + case (LandType) ! raster grid cell (i,j) is Land + + if(land_id(l)==0) then + ! raster grid cell is first Land type seen for this EASE grid cell, tile ID has not yet been set + ! recall: tile IDs for Land: [1:n_land], and l_index was initalized to 0 above + l_index = l_index + 1 + tileid_index(i,j)= l_index + else + ! no action needed (tile ID does not depend on number of contributing Land raster grid cells) + endif + + ! sum up area and (area-weighted) elevation (only over raster grid cells of type land!) + + tile_ele( tileid_index(i,j)) = tile_ele( tileid_index(i,j)) + q0(i,j) * pix_area ! q0 = elevation + + ! tile_area_land should be obsolete because identical to tile_area(1:n_land) + !tile_area_land(tileid_index(i,j)) = tile_area_land(tileid_index(i,j)) + pix_area + + case default + + print *,'ERROR: unknown tile type value in veg(i,j): ', veg(i,j), ' STOPPING.' + stop + + end select + + ! sum up area of raster grid cells contributing to each tile (land or lake or landice) + + tile_area(tileid_index(i,j)) = tile_area(tileid_index(i,j)) + pix_area + + my_land(tileid_index(i,j)) = l ! store 1-dim index of all EASE grid cells + all_id( tileid_index(i,j)) = j*ND_raster + i ! store 1-dim index of all 30-arcsec raster grid cells; BUG???? - ! computing tile area/elevation - ! ----------------------------- + ! BUG??? all_id stores only the 1-dim index of the 30-arcsec raster grid cells that is the + ! last one to contribute to the EASE tile specified by tileid_index(i,j) + ! This does not seem to be the desired dominant ID across the EASE grid cell + + endif ! raster grid cell has tile type other than ocean + + ! compute total area of raster grid cells of any tile type that contribute to EASE grid cell + + ease_grid_area(l) = ease_grid_area(l) + pix_area - if((tileid_index(i,j) > 0).and.(tileid_index(i,j) <= nland_cells))then - tile_ele(tileid_index(i,j)) = tile_ele(tileid_index(i,j)) + q0(i,j) * & - pix_area - tile_area_land(tileid_index(i,j)) = tile_area_land(tileid_index(i,j)) + & - pix_area - endif + endif ! raster grid cell is located within lat band of EASE grid + end do end do @@ -613,33 +742,36 @@ PROGRAM mkEASETilesParam deallocate(water_id) deallocate(ice_id ) - tile_ele = tile_ele/tile_area_land + tile_ele = tile_ele/tile_area(1:n_land) ! finalize tile elevation ! adjustment Global Mean Topography to 614.649 (615.662 GTOPO 30) m ! ----------------------------------------------------------------- - sum1=0. - sum2=0. + mean_land_ele=0. do j=1,l_index - sum1 = sum1 + tile_ele(j)*tile_area(j) + mean_land_ele = mean_land_ele + tile_ele(j)*tile_area(j) enddo - if(sum1/sum(tile_area(1:l_index)).ne. 614.649D0 ) then - print *,'Global Mean Elevation (over land): ', sum1/sum(tile_area(1:l_index)) - tile_ele =tile_ele*(614.649D0 / (sum1/sum(tile_area(1:l_index)))) - sum1=0. - sum2=0. + mean_land_ele = mean_land_ele/sum(tile_area(1:l_index)) + + if(mean_land_ele .ne. 614.649D0 ) then + print *,'Global Mean Elevation (over land): ', mean_land_ele + tile_ele = tile_ele*(614.649D0/mean_land_ele) + ! verify that adjust elevation is correct + mean_land_ele=0. do j=1,l_index - sum1 = sum1 + tile_ele(j)*tile_area(j) + mean_land_ele = mean_land_ele + tile_ele(j)*tile_area(j) enddo - print *,'Global Mean Elevation after scaling to SRTM : ',sum1/sum(tile_area(1:l_index)) + print *,'Global Mean Elevation after scaling to SRTM : ',mean_land_ele/sum(tile_area(1:l_index)) endif - print *,'Total Land Area :', sum(tile_area(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000., & - sum(tile_area_land(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. - - print *,'Creating ... ', trim(gfile)//'.rst' + print *,'Total Land Area :', sum(tile_area(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. ! , & + !sum(tile_area_land(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. - !------------------------------------------- + !------------------------------------------- + ! + ! write *.rst + + print *,'Writing ... ', trim(gfile)//'.rst' open (10, file ='rst/'//trim(gfile)//'.rst',form='unformatted',status='unknown', & action='write') @@ -650,94 +782,150 @@ PROGRAM mkEASETilesParam close (10,status='keep') - print *,'Creating ... ', trim(gfile)//'.til ,catchment.def' - - !----------------------------------------------------------- - + !----------------------------------------------------------- + + ! write catchment.def and *.til files + + print *,'Writing ... ', trim(gfile)//'.til and catchment.def' + open (11,file='clsm/catchment.def', & form='formatted',status='unknown') write(11,*)l_index open (10, file ='til/'//trim(gfile)//'.til',form='formatted',status='unknown',action='write') - write (10,*)i_index,SRTM_maxcat, nc, nr - write (10,*)1 - write (10,*)EASELabel - write (10,*)nc_ease - write (10,*)nr_ease - ! write (10,*)'NO-OCEAN' - ! write (10,*) -9999 - ! write (10,*) -9999 + write (10,*) i_index,SRTM_maxcat, nc, nr + write (10,*) 1 + write (10,*) EASELabel + write (10,*) nc_ease + write (10,*) nr_ease + + ! write (10,*)'NO-OCEAN' + ! write (10,*) -9999 + ! write (10,*) -9999 + + dx_ease = 180./real(nc_ease) do l=1,i_index - ig = my_land(l)-ND*(my_land(l)/ND) - jg = my_land(l)/ND + !ig = my_land(l)-NDND*(my_land(l)/NDND) + !jg = my_land(l)/NDND + ! get row and column indices from 1-dim index (that is, invert l=(ind_row-1)*NDND+ind_col) + + jg = (my_land(l)-1)/NDND + 1 ! = ind_row (1-based) [note integer division] + + ig = my_land(l) - NDND*(jg-1) ! = ind_col (1-based) + + ! extract original PfafID from catid_index + ! + ! BUG??? catid_index is in 30-arcsec raster space and contains the dominant PfafID + ! (or LakeID, or LandiceID) within the 30-arcsec raster grid cell, + ! where dominant is w.r.t. the 10-arcsec mask with PfafIDs + ! It seems that here we just pick the PfafID associated with a single + ! 30-arcsec raster grid cell within the EASE tile, which happens to be the + ! last of the 30-arcsec raster grid cells that contribute to the EASE tile + ! + + + + cindex= catid_index(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then pfaf = cindex else - print *,'Only mask file we support with EASE grid is GEOS5_10arcsec_mask' - call MPI_Abort(MPI_COMM_WORLD,error_code,status) - ! pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) + print *,'MaskFile = ', trim(MaskFile) + print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' + stop + ! pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) endif + + if ((l>l_index) .and. (l<=w_index)) typ = 19 ! Lake tile + + if ( l>w_index ) typ = 20 ! Landice tile + + if (l<=l_index) then - if ((l > l_index).and.(l <= w_index)) typ =19 - if (l > w_index) typ = 20 + typ = 100 ! Land tile - if (l <= l_index) then - typ = 100 - call EASE_inverse (EASELabel, real(ig-1),real(jg-1), clat, clon) + ! get min/max lat/lon of EASE grid cell + ! BUG: This is *not* the desired min/max lat/lon of the land tile!!! + + call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) - mnx = clon - 180./real(nc_ease) - mxx = clon + 180./real(nc_ease) + mnx = clon - dx_ease + mxx = clon + dx_ease jgv = real(jg-1) + 0.5 - call EASE_inverse (EASELabel, real(ig-1),jgv, clat, clon) + call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) mny = clat jgv = real(jg-1) - 0.5 - call EASE_inverse (EASELabel, real(ig-1),jgv, clat, clon) + call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) mxy = clat + + ! write tile properties into catchment.def file - write (11,'(i10,i8,5(2x,f9.4), i4)')l,pfaf,mnx,mxx,mny,mxy,tile_ele(l) + write (11,'(i10,i8,5(2x,f9.4), i4)') l, pfaf, mnx, mxx, mny, mxy, tile_ele(l) endif - call EASE_inverse (EASELabel, real(ig-1), real(jg-1), clat, clon) + ! get area fraction of tile within EASE grid cell + ! NOTE: the area of the EASE grid cell here has to be the sum of the areas of the + ! contributing raster grid cells, which is *not* the same for all EASE grid cells + + call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) + + fr_gcm = tile_area(l) / ease_grid_area((jg-1)*NDND+ig) + + ! write tile properties into *.til file - fr_gcm= tile_area(l)/ease_grid_area(jg*ND + ig) - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - write(10,'(i10,i9,2f10.4,2i6,f19.12,i10,i15,e13.4)') & - typ,pfaf,clon,clat,ig-1,jg-1,fr_gcm ,pfaf,SRTM_catid(cindex) + + ! Note: If-condition uses MaskFile as a proxy for format & specs of *.til file (columns written) + + ! pfaf = running *index* of Pfafstetter catchment (1:SRTM_maxcat) + ! SRTM_catid = 12-digit Pfafstetter code (encoding network/routing info) + + write(10,'(i10,i9,2f10.4,2i6,f19.12,i10,i15,e13.4)') & + typ, pfaf, clon, clat, ig-1, jg-1, fr_gcm, pfaf, SRTM_catid(cindex) else - write(10,'(i10,i9,2f10.4,2i5,f19.12,i10,e13.4,i8)') & - typ,pfaf,clon,clat,ig-1,jg-1,fr_gcm ,cindex + + ! write statement below was used with MaskFile prior to availability of SRTM-based Pfafstetter catchments + ! obsolete with EASE grid + ! - reichle, 15 Jun 2023 + + print *,'MaskFile = ', trim(MaskFile) + print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' + stop +! +! write(10,'(i10,i9,2f10.4,2i5,f19.12,i10,e13.4,i8)') & +! typ,pfaf,clon,clat,ig-1,jg-1,fr_gcm ,cindex endif end do - close (10,status='keep') - close (11,status='keep') - - deallocate (tileid_index,catid_index,veg) - deallocate (tile_area, ease_grid_area, tile_ele, tile_area_land, my_land, all_id) - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - - print *,'Creating SMAP-Catch_TransferData.nc files.' - - !--------------------------------------------------- - - deallocate (SRTM_CatchArea, SRTM_catid, SRTM_catid_r8) - - endif + close(10,status='keep') + close(11,status='keep') + + deallocate( tileid_index, catid_index,veg ) + deallocate( tile_area, ease_grid_area, tile_ele, tile_area_land, my_land, all_id ) + + ! Commented out "empty" if-block. -rreichle, 15 Jun 2023 +! +! if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then +! +! print *,'Creating SMAP-Catch_TransferData.nc files.' +! +! !--------------------------------------------------- +! +! deallocate (SRTM_CatchArea, SRTM_catid, SRTM_catid_r8) +! +! endif ! create Grid2Catch transfer file ! ------------------------------- From bd375d547f0432d96fb41bff9e36e934eff7dbe2 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 17 Jun 2023 14:37:20 -0400 Subject: [PATCH 06/16] fixed minor errors in previous commit --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 23 ++++++++++--------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index c985952fb..1fd234e27 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -36,12 +36,21 @@ PROGRAM mkEASETilesParam use netcdf implicit none + + integer, parameter :: nc_esa = 129600 ! number of cols in 10-arcsec ESA mask file + integer, parameter :: nr_esa = 64800 ! number of rows in 10-arcsec ESA mask file + + ! define tile types used for processing here (values may be from ESA mask?) + + integer, parameter :: OceanType = 0 + integer, parameter :: LandType = 1 ! land type used for processing here; in GEOS, land tiles are type=100 + integer, parameter :: LakeType = 10 ! lake type used for processing here; in GEOS, lake tiles are type= 19 + integer, parameter :: IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 + integer :: i, j, ig, jg, n integer :: NC, NR, N_ease_grid_cells, NDND, ND_raster - integer, parameter :: nc_esa = 129600, nr_esa = 64800 - ! For regridding integer, allocatable, dimension(:,:), target :: geos_msk @@ -70,7 +79,6 @@ PROGRAM mkEASETilesParam !INTEGER*8 :: PFAF_CODE integer :: l, l_index, i_index, w_index, typ, pfaf, cindex - integer :: LandType, LakeType, IceType, OceanType real :: clat, clon, r_ease, s_ease real :: fr_gcm integer :: ind_col, ind_row, status, ncid, varid, nciv !, DOM_INDX @@ -182,13 +190,6 @@ PROGRAM mkEASETilesParam water_id = 0 ice_id = 0 - ! define tile types used for processing here (possibly from ESA mask) - - OceanType = 0 - IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 - LakeType = 10 ! lake type used for processing here; in GEOS, lake tiles are type= 19 - LandType = 1 ! land type used for processing here; in GEOS, land tiles are type=100 - ! Step size for conversion of 2-dim indexing into 1-dim indexing: l = ind_row*NDND + ind_col ! If conversion was simply to get from 2-dim to 1-dim indexing for EASE grid cells, ! NDND=nc_ease would suffice @@ -913,7 +914,7 @@ PROGRAM mkEASETilesParam close(11,status='keep') deallocate( tileid_index, catid_index,veg ) - deallocate( tile_area, ease_grid_area, tile_ele, tile_area_land, my_land, all_id ) + deallocate( tile_area, ease_grid_area, tile_ele, my_land, all_id ) ! Commented out "empty" if-block. -rreichle, 15 Jun 2023 ! From fa5fce3561db2472373d039e636ebe6931f593d3 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 17 Jun 2023 15:20:34 -0400 Subject: [PATCH 07/16] fixed netcdf error print statements (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 1fd234e27..b542d7fbe 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -270,12 +270,12 @@ PROGRAM mkEASETilesParam ! read list of Pfafstetter catchment IDs ('PfafID') from 10-arcsec mask file into variable 'SRTM_catid[_r8]' ! [vector of length SRTM_maxcat] - status = NF90_OPEN (trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/GEOS5_10arcsec_mask.nc', NF90_NOWRITE, ncid) - status = nf90_inq_varid(ncid, name='PfafID', varid=varid) - status = nf90_get_var(ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/)) + status = NF90_OPEN( trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/GEOS5_10arcsec_mask.nc', NF90_NOWRITE, ncid ) + status = nf90_inq_varid( ncid, name='PfafID', varid=varid ) + status = nf90_get_var( ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/) ) if(status /=0) then - PRINT *, NF90_STRERROR(STATUS) - print *, 'Problem with NF90_OPEN',trim(MaskFile) + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'Problem with NF90_OPEN(): ', trim(MaskFile) endif SRTM_catid = int8(SRTM_catid_r8) ! convert data to integer*8 -- contains 12-digit Pfaf code SRTM_catid (SRTM_maxcat + 1) = 190000000 ! append ID for Lake type @@ -289,7 +289,7 @@ PROGRAM mkEASETilesParam !i1 = 0 ! count # of 30-arcsec pixels - NOT USED - status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) + status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) do j=1,nr @@ -298,12 +298,12 @@ PROGRAM mkEASETilesParam ! read slice of variable 'CatchIndex' from 10-arcsec mask file into variable 'geos_msk' ! [2d-array: 129600-by-3] - status = NF90_GET_VAR (ncid, varid, geos_msk, (/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/)) ! Read 10-arcsec rows that lie within the raster row 'j' + status = NF90_GET_VAR( ncid, varid, geos_msk, (/1,(j-1)*dy_esa +1/), (/nc_esa,dy_esa/) ) ! Read 10-arcsec rows that lie within the raster row 'j' !status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' if(status /=0) then - PRINT *, NF90_STRERROR(STATUS) - print *, 'Problem with NF_GET_VARA_INT',trim(MaskFile),status + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'Problem with NF_GET_VAR(): ', trim(MaskFile), status endif do i = 1,nc @@ -317,8 +317,8 @@ PROGRAM mkEASETilesParam subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains ESA pixels that lie within the raster grid cell at i,j if(maxval (subset) > SRTM_maxcat) then - where (subset == 190000000) subset = SRTM_maxcat + 1 ! Lake type (convert ID from 190000000 to SRTM_maxcat+1) - where (subset == 200000000) subset = SRTM_maxcat + 2 ! Landice type (convert ID from 200000000 to SRTM_maxcat+2) + where (subset == 190000000) subset = SRTM_maxcat + 1 ! Lake type (convert ID from 190000000 to SRTM_maxcat+1) + where (subset == 200000000) subset = SRTM_maxcat + 2 ! Landice type (convert ID from 200000000 to SRTM_maxcat+2) endif if (maxval(subset) > 0) then ! check if there are Non-ocean ESA pixels From dc438cdbfe24f47ea436239ab3d9f3d3a3acd123 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 17 Jun 2023 16:29:02 -0400 Subject: [PATCH 08/16] clarified print statements for NF90 status checks; clarified resolution of data in GEOS5_10arcsec_mask.nc (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 35 ++++++++++++++++--- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index b542d7fbe..8cd5c844c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -47,7 +47,6 @@ PROGRAM mkEASETilesParam integer, parameter :: LakeType = 10 ! lake type used for processing here; in GEOS, lake tiles are type= 19 integer, parameter :: IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 - integer :: i, j, ig, jg, n integer :: NC, NR, N_ease_grid_cells, NDND, ND_raster @@ -111,6 +110,8 @@ PROGRAM mkEASETilesParam character(len=128) :: usage1, usage2 character(len=128) :: Iam = "mkEASETilesParam" + character(len=512) :: fname_mask + ! -------------------------------------------------------------------------------------- call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) @@ -226,6 +227,9 @@ PROGRAM mkEASETilesParam ! ESA/SRTM mask file is 10-arcsec resolution; coarsen here to 30-arcsec raster grid: + ! NOTE: Only coastlines (and lake shorelines?) are at 10-arcsec resolution. + ! The original Pfaf catchments are at 1-arcmin resolution (~2 km). + nc = 43200 ! Number of rows in raster file nr = 21600 ! Number of columns in raster file @@ -270,13 +274,30 @@ PROGRAM mkEASETilesParam ! read list of Pfafstetter catchment IDs ('PfafID') from 10-arcsec mask file into variable 'SRTM_catid[_r8]' ! [vector of length SRTM_maxcat] - status = NF90_OPEN( trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/GEOS5_10arcsec_mask.nc', NF90_NOWRITE, ncid ) + fname_mask = trim(MAKE_BCS_INPUT_DIR) // '/shared/mask/GEOS5_10arcsec_mask.nc' + + print *, 'Opening ', trim(fname_mask) + + status = NF90_OPEN( fname_mask, NF90_NOWRITE, ncid ) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'Problem with NF90_OPEN(): ', trim(fname_mask), ' ', ncid, status + else + print *, 'ncid=', ncid + endif + status = nf90_inq_varid( ncid, name='PfafID', varid=varid ) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'Problem with NF90_INQ_VARID(): PfafID ', ncid, varid, status + endif + status = nf90_get_var( ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/) ) if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF90_OPEN(): ', trim(MaskFile) + print *, 'Problem with NF90_GET_VAR(): PfafID ', ncid, varid, SRTM_maxcat, status endif + SRTM_catid = int8(SRTM_catid_r8) ! convert data to integer*8 -- contains 12-digit Pfaf code SRTM_catid (SRTM_maxcat + 1) = 190000000 ! append ID for Lake type SRTM_catid (SRTM_maxcat + 2) = 200000000 ! append ID for Landice type @@ -290,6 +311,10 @@ PROGRAM mkEASETilesParam !i1 = 0 ! count # of 30-arcsec pixels - NOT USED status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'Problem with NF90_INQ_VARID(): CatchIndex ', ncid, varid, status + endif do j=1,nr @@ -298,12 +323,12 @@ PROGRAM mkEASETilesParam ! read slice of variable 'CatchIndex' from 10-arcsec mask file into variable 'geos_msk' ! [2d-array: 129600-by-3] - status = NF90_GET_VAR( ncid, varid, geos_msk, (/1,(j-1)*dy_esa +1/), (/nc_esa,dy_esa/) ) ! Read 10-arcsec rows that lie within the raster row 'j' + status = NF90_GET_VAR( ncid, varid, geos_msk, (/1,(j-1)*dy_esa+1/), (/nc_esa,dy_esa/) ) ! Read 10-arcsec rows that lie within the raster row 'j' !status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF_GET_VAR(): ', trim(MaskFile), status + print *, 'Problem with NF_GET_VAR(): CatchIndex ', ncid, varid, j, dy_esa, nc_esa, status endif do i = 1,nc From 6c3753ce371298fdd5371e89db1134cef35d9afc Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 17 Jun 2023 17:25:58 -0400 Subject: [PATCH 09/16] fixed bug introduced in earlier clean-up commit (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 24 +++++++++++-------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 8cd5c844c..4650e3bd3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -114,7 +114,7 @@ PROGRAM mkEASETilesParam ! -------------------------------------------------------------------------------------- - call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) + call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR ) usage1 = 'USAGE : bin/mkEASETilesParam.x -ease_label EASELabel ' usage2 = ' where EASELabel = *EASEv[x]_M[yy]*, x={1,2}, yy={01,03,09,25,36}' @@ -615,9 +615,10 @@ PROGRAM mkEASETilesParam n_landlakelandice = n_land + n_lake + n_landice - print *,'# of Land tiles: ', n_land - print *,'# of Lake tiles: ', n_lake - print *,'# of Landice tiles: ', n_landice + print *,'# of Land tiles: ', n_land + print *,'# of Lake tiles: ', n_lake + print *,'# of Landice tiles: ', n_landice + print *,'# of Land+Lake+Landice tiles: ', n_landlakelandice l_index = 0 ! start index for land EASE cells w_index = n_land ! start index for lake EASE cells @@ -700,8 +701,9 @@ PROGRAM mkEASETilesParam if(water_id(l)==0) then ! raster grid cell is first Lake type seen for this EASE grid cell, tile ID has not yet been set ! recall: tile IDs for Lake: [(n_land+1):(n_land+n_lake)], and w_index was initalized to n_land above - w_index = w_index + 1 - tileid_index(i,j)= w_index + w_index = w_index + 1 + water_id(l) = w_index ! needed so if condition will be false for next raster grid cell of same type + tileid_index(i,j) = w_index else ! no action needed (tile ID does not depend on number of contributing Lake raster grid cells) endif @@ -711,8 +713,9 @@ PROGRAM mkEASETilesParam if(ice_id(l)==0) then ! raster grid cell is first Landice type seen for this EASE grid cell, tile ID has not yet been set ! recall: tile IDs for Landice: [(n_land+n_lake+1):n_landlakelandice], and i_index was initalized to n_land+n_lake above - i_index = i_index + 1 - tileid_index(i,j)= i_index + i_index = i_index + 1 + ice_id(l) = i_index ! needed so if condition will be false for next raster grid cell of same type + tileid_index(i,j) = i_index else ! no action needed (tile ID does not depend on number of contributing Landice raster grid cells) endif @@ -722,8 +725,9 @@ PROGRAM mkEASETilesParam if(land_id(l)==0) then ! raster grid cell is first Land type seen for this EASE grid cell, tile ID has not yet been set ! recall: tile IDs for Land: [1:n_land], and l_index was initalized to 0 above - l_index = l_index + 1 - tileid_index(i,j)= l_index + l_index = l_index + 1 + land_id(l) = l_index ! needed so if condition will be false for next raster grid cell of same type + tileid_index(i,j) = l_index else ! no action needed (tile ID does not depend on number of contributing Land raster grid cells) endif From 8492415ad4db2eff01831182c88c141044af676a Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sat, 17 Jun 2023 20:15:18 -0400 Subject: [PATCH 10/16] removed obsolete $THISGRID from make_bcs --- .../Utils/Raster/makebcs/make_bcs | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs index 662206d74..4bfe47cfa 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs @@ -1403,7 +1403,6 @@ if( $grid == ease ) then set BCNAME = ${EASEVERSION}_${HRCODE} set BCDIR = $EXPDIR/$OUTDIR/$BCNAME.scratch set BCJOB = $BCDIR/$BCNAME.j -# set THISGRID = ${EASEVERSION}-${HRCODE} set nfiles = `find $EXPDIR -maxdepth 5 -name ${BCNAME}".j" | wc -l` if( $nfiles >= 1 ) then @@ -1444,17 +1443,6 @@ cd data cd ../ limit stacksize unlimited -#if ( $EASEVERSION == EASEv2 ) then - ## This section was used to make Irrigated Tiles - ##if(${HRCODE} == M09 | ${HRCODE} == M36) then - ## bin/mkLandRaster.x -x ${NX} -y ${NY} -v -t ${NT} - ## bin/mkEASETilesParam.x -ease_label ${HRCODE} -pfaf_til T - ## bin/CombineRasters.x -f 0 -t 232000000 ${THISGRID} Pfafstetter > /dev/null - ## bin/CombineRasters.x -t 232000000 ${THISGRID} ${THISGRID}-Pfafstetter - ## /bin/mv til/${THISGRID}_${THISGRID}-Pfafstetter.til til/${THISGRID}_${THISGRID}-Pfafstetter.ind - ##endif -#endif - setenv MASKFILE ${MASKFILE} setenv OMP_NUM_THREADS 1 bin/mkEASETilesParam.x -ease_label ${BCNAME} @@ -1468,7 +1456,6 @@ bin/create_README.csh /bin/mv clsm clsm.${IM}x${JM} /bin/cp til/${EASEVERSION}_${HRCODE}_${RS}.til clsm.${IM}x${JM} -##/bin/cp til//${THISGRID}_${THISGRID}-Pfafstetter.TIL clsm.${IM}x${JM} cd clsm.${IM}x${JM} /bin/mv irrig.dat irrigation_${RS}_DE.dat @@ -1549,7 +1536,6 @@ mkdir -p $BCNAME/land/${EASEVERSION}_${HRCODE} #mkdir -p IRRIG/$BCNAME/clsm IRRIG/$BCNAME/rst #bin/mkIrrigTiles.x -x 43200 -y 21600 -b $BCNAME -t ${EASEVERSION}_${HRCODE}_${RS}.til -r _${RS}_DE -p $IRRIGTHRES -#/bin/cp $BCNAME/${THISGRID}_${THISGRID}-Pfafstetter.TIL IRRIG/$BCNAME/. _EOF_ From a097bf4518ed59603e91662e6efb358fd0683ad7 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Sun, 18 Jun 2023 15:43:46 -0400 Subject: [PATCH 11/16] fixed bug introduced during earlier clean-up commit; improved error print statements (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 41 +++++++++++-------- 1 file changed, 24 insertions(+), 17 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 4650e3bd3..11e6d9d14 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -281,7 +281,9 @@ PROGRAM mkEASETilesParam status = NF90_OPEN( fname_mask, NF90_NOWRITE, ncid ) if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF90_OPEN(): ', trim(fname_mask), ' ', ncid, status + print *, 'ERROR -- NF90_OPEN(): ', trim(fname_mask), ' ', ncid, status + print *, trim(Iam), '.x STOPPING.' + stop else print *, 'ncid=', ncid endif @@ -289,13 +291,17 @@ PROGRAM mkEASETilesParam status = nf90_inq_varid( ncid, name='PfafID', varid=varid ) if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF90_INQ_VARID(): PfafID ', ncid, varid, status + print *, 'ERROR -- NF90_INQ_VARID(): PfafID ', ncid, varid, status + print *, trim(Iam), '.x STOPPING.' + stop endif status = nf90_get_var( ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/) ) if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF90_GET_VAR(): PfafID ', ncid, varid, SRTM_maxcat, status + print *, 'ERROR -- NF90_GET_VAR(): PfafID ', ncid, varid, SRTM_maxcat, status + print *, trim(Iam), '.x STOPPING.' + stop endif SRTM_catid = int8(SRTM_catid_r8) ! convert data to integer*8 -- contains 12-digit Pfaf code @@ -313,7 +319,9 @@ PROGRAM mkEASETilesParam status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF90_INQ_VARID(): CatchIndex ', ncid, varid, status + print *, 'ERROR -- NF90_INQ_VARID(): CatchIndex ', ncid, varid, status + print *, trim(Iam), '.x STOPPING.' + stop endif do j=1,nr @@ -328,7 +336,9 @@ PROGRAM mkEASETilesParam if(status /=0) then PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'Problem with NF_GET_VAR(): CatchIndex ', ncid, varid, j, dy_esa, nc_esa, status + print *, 'ERROR -- NF_GET_VAR(): CatchIndex ', ncid, varid, j, dy_esa, nc_esa, status + print *, trim(Iam), '.x STOPPING.' + stop endif do i = 1,nc @@ -703,10 +713,9 @@ PROGRAM mkEASETilesParam ! recall: tile IDs for Lake: [(n_land+1):(n_land+n_lake)], and w_index was initalized to n_land above w_index = w_index + 1 water_id(l) = w_index ! needed so if condition will be false for next raster grid cell of same type - tileid_index(i,j) = w_index - else - ! no action needed (tile ID does not depend on number of contributing Lake raster grid cells) - endif + end if + + tileid_index(i,j) = water_id(l) case (IceType) ! raster grid cell (i,j) is Landice @@ -715,10 +724,9 @@ PROGRAM mkEASETilesParam ! recall: tile IDs for Landice: [(n_land+n_lake+1):n_landlakelandice], and i_index was initalized to n_land+n_lake above i_index = i_index + 1 ice_id(l) = i_index ! needed so if condition will be false for next raster grid cell of same type - tileid_index(i,j) = i_index - else - ! no action needed (tile ID does not depend on number of contributing Landice raster grid cells) - endif + end if + + tileid_index(i,j) = ice_id(l) case (LandType) ! raster grid cell (i,j) is Land @@ -727,10 +735,9 @@ PROGRAM mkEASETilesParam ! recall: tile IDs for Land: [1:n_land], and l_index was initalized to 0 above l_index = l_index + 1 land_id(l) = l_index ! needed so if condition will be false for next raster grid cell of same type - tileid_index(i,j) = l_index - else - ! no action needed (tile ID does not depend on number of contributing Land raster grid cells) - endif + end if + + tileid_index(i,j) = land_id(l) ! sum up area and (area-weighted) elevation (only over raster grid cells of type land!) From 6ce84112c82cba934b665b2a4895d10d03d395a9 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 19 Jun 2023 09:59:17 -0400 Subject: [PATCH 12/16] more efficient 1-dim indexing for raster grid; added checks for number of tiles (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 57 ++++++++++++++----- 1 file changed, 42 insertions(+), 15 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 11e6d9d14..58de108ca 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -48,7 +48,7 @@ PROGRAM mkEASETilesParam integer, parameter :: IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 integer :: i, j, ig, jg, n - integer :: NC, NR, N_ease_grid_cells, NDND, ND_raster + integer :: nc, nr, N_ease_grid_cells, NDND, ND_raster ! For regridding @@ -191,14 +191,10 @@ PROGRAM mkEASETilesParam water_id = 0 ice_id = 0 - ! Step size for conversion of 2-dim indexing into 1-dim indexing: l = ind_row*NDND + ind_col - ! If conversion was simply to get from 2-dim to 1-dim indexing for EASE grid cells, - ! NDND=nc_ease would suffice - ! Here, NDND is computed as a power of 10, leaving room for "empty" indices that may or may not be needed later. - ! - reichle, 15 Jun 2023 + ! conversion from 2-dim to 1-dim index for EASE grid: l = (ind_row-1)*NDND + ind_col !NDND = 10*10**(nint(log10(1.*nr_ease))) - + NDND = nc_ease ! Check for the 10 arc-sec MaskFile @@ -670,8 +666,12 @@ PROGRAM mkEASETilesParam ! global nc x nr array of tileid_index(nc,nr) contains corresponding tile_index values which ! is derived in the below loop - ND_raster = 10*10**(nint(log10(1.*NR))) + ! conversion from 2-dim index to 1-dim index for raster grid: all_ind(l) = (j-1)*ND_raster + i + + !ND_raster = 10*10**(nint(log10(1.*NR))) + ND_raster = nc + !i2 = 1 -- NOT USED do i = 1 ,nc @@ -757,8 +757,8 @@ PROGRAM mkEASETilesParam tile_area(tileid_index(i,j)) = tile_area(tileid_index(i,j)) + pix_area - my_land(tileid_index(i,j)) = l ! store 1-dim index of all EASE grid cells - all_id( tileid_index(i,j)) = j*ND_raster + i ! store 1-dim index of all 30-arcsec raster grid cells; BUG???? + my_land(tileid_index(i,j)) = l ! for this tile, store 1-dim index for EASE grid cells + all_id( tileid_index(i,j)) = (j-1)*ND_raster + i ! for this tile, store 1-dim index for raster grid cells - last in prevails! ! BUG??? all_id stores only the 1-dim index of the 30-arcsec raster grid cells that is the ! last one to contribute to the EASE tile specified by tileid_index(i,j) @@ -774,7 +774,29 @@ PROGRAM mkEASETilesParam end do end do + + ! verify l_index, w_index, i_index against tile counts from first loop + if (l_index/= n_land ) then + print *, 'ERROR: l_index = ', l_index, ' must now match # of Land tiles = ', n_land + print *, 'STOPPING.' + stop + end if + + if (w_index/=(n_land+n_lake) ) then + print *, 'ERROR: w_index = ', w_index, ' must now match # of Land+Lake tiles = ', n_land+n_lake + print *, 'STOPPING.' + stop + + end if + + if (i_index/= n_landlakelandice ) then + print *, 'ERROR: i_index = ', i_index, ' must now match # of Land+Lake+Landice tiles = ', n_landlakelandice + print *, 'STOPPING.' + stop + + end if + deallocate(land_id, q0) deallocate(water_id) deallocate(ice_id ) @@ -847,7 +869,7 @@ PROGRAM mkEASETilesParam !ig = my_land(l)-NDND*(my_land(l)/NDND) !jg = my_land(l)/NDND - ! get row and column indices from 1-dim index (that is, invert l=(ind_row-1)*NDND+ind_col) + ! for EASE grid: get row and column indices from 1-dim index (that is, invert l=(ind_row-1)*NDND+ind_col) jg = (my_land(l)-1)/NDND + 1 ! = ind_row (1-based) [note integer division] @@ -863,11 +885,16 @@ PROGRAM mkEASETilesParam ! last of the 30-arcsec raster grid cells that contribute to the EASE tile ! + ! for raster grid: get row and column indices from 1-dim index (that is, invert all_id(l)=(j-1)*ND_raster + i - - - cindex= catid_index(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) - + j = (all_id(l) -1)/ND_raster + 1 + + i = all_id(l) - ND_raster*(j-1) + + !cindex= catid_index(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) + + cindex = catid_index(i,j) + if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then pfaf = cindex else From dae09d97f7b8b5a39e463a0c449b91175e8da1be Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Mon, 19 Jun 2023 12:20:32 -0400 Subject: [PATCH 13/16] additional cleanup: - cleaned-up 1-dim indexing for EASE and raster grid cells - renamed loop counter variables for clarity - cleaned up hardwired target for global mean elevation - edited comments for clarity - removed obsolete code --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 1588 ++++++++--------- 1 file changed, 782 insertions(+), 806 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 58de108ca..2fcff34f0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -5,7 +5,12 @@ PROGRAM mkEASETilesParam ! This program constructs land, landice, and lake tiles for EASE grid tile spaces such as ! those used for the SMAP Level-4 products and other offline projects - + ! + ! Files generated by this program: + ! [EASElabel]_[EASEncol]x[EASEnrow].til + ! [EASElabel]_[EASEncol]x[EASEnrow].rst + ! catchment.def + ! This program resulted from the merger and cleanup of mkSMAPTilesPara.F90 ! and mkSMAPTilesPara_v2.F90 in September 2022. ! Before the merger and cleanup, the EASE grid parameters were hard-coded here. @@ -14,7 +19,7 @@ PROGRAM mkEASETilesParam ! ! - wjiang + reichle, 21 Sep 2022 - ! additional major cleanup (reichle, 15 Jun 2023) + ! Additional *major* cleanup (reichle, 15 Jun 2023) ! - rewrote some code blocks for clarity and efficiency ! - removed obsolete variable declarations ! - removed obsolete use statements @@ -25,419 +30,388 @@ PROGRAM mkEASETilesParam ! - added comments ! - white-space changes for improved readability - use EASE_conv, only : EASE_extent, EASE_convert, EASE_inverse - use rmTinyCatchParaMod, only : i_raster, j_raster, SRTM_maxcat - use rmTinyCatchParaMod, only : RegridRasterReal - use rmTinyCatchParaMod, only : MAKE_BCS_INPUT_DIR - use process_hres_data, only : histogram - use MAPL_SortMod - use MAPL_ConstantsMod - use MAPL_ExceptionHandling - use netcdf - - implicit none - - integer, parameter :: nc_esa = 129600 ! number of cols in 10-arcsec ESA mask file - integer, parameter :: nr_esa = 64800 ! number of rows in 10-arcsec ESA mask file - - ! define tile types used for processing here (values may be from ESA mask?) - - integer, parameter :: OceanType = 0 - integer, parameter :: LandType = 1 ! land type used for processing here; in GEOS, land tiles are type=100 - integer, parameter :: LakeType = 10 ! lake type used for processing here; in GEOS, lake tiles are type= 19 - integer, parameter :: IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 - - integer :: i, j, ig, jg, n - integer :: nc, nr, N_ease_grid_cells, NDND, ND_raster - - ! For regridding - - integer, allocatable, dimension(:,:), target :: geos_msk - - REAL, allocatable, DIMENSION(:) :: loc_val - INTEGER, ALLOCATABLE, DIMENSION(:) :: density, loc_int - logical, allocatable, dimension(:) :: unq_mask - integer, dimension(:,:), pointer :: subset - - integer :: dx_esa, dy_esa, NBINS, NPLUS - - integer*8, allocatable, dimension(:) :: SRTM_catid - real(kind=8),allocatable, dimension(:) :: SRTM_catid_r8 - - integer, allocatable, dimension(:,:), target :: tileid_index, catid_index - - ! integer, allocatable, dimension(:,:) :: catid, iaster - - integer, allocatable, dimension(:) :: land_id, water_id, ice_id - integer, allocatable, dimension(:) :: my_land, all_id - real, allocatable, dimension(:) :: ease_grid_area, tile_area !, SRTM_CatchArea - integer*1, allocatable, dimension(:,:) :: veg !, i2aster - real*4, allocatable, dimension(:,:) :: q0, raster - REAL, allocatable, dimension(:) :: tile_ele - - !INTEGER*8 :: PFAF_CODE - - integer :: l, l_index, i_index, w_index, typ, pfaf, cindex - real :: clat, clon, r_ease, s_ease - real :: fr_gcm - integer :: ind_col, ind_row, status, ncid, varid, nciv !, DOM_INDX - integer :: n_land, n_lake, n_landice, n_landlakelandice - - ! REAL (kind=8), PARAMETER :: RADIUS=6378137.0,pi=3.14159265358979323846 - - !character*100 :: veg_class(12) - - character*200 :: gfile, gtopo30 - integer :: nc_ease, nr_ease, N_args, command_argument_count - REAL :: dx, dy, d2r, lats, mnx, mxx, mny, mxy, jgv, pix_area ! VDUM - real :: dx_ease, mean_land_ele - character(40) :: arg, EASElabel_ - - character(len=:), allocatable :: EASElabel - - !character*200 :: tmpstring, tmpstring1, tmpstring2 - - logical :: regrid = .false. - character*128 :: MaskFile - - !logical :: pfaf_til = .false. - - character*1 :: PF - - !character(len=6) :: EASE_Version - - character(len=10) :: nc_string, nr_string - character(len=128) :: usage1, usage2 - character(len=128) :: Iam = "mkEASETilesParam" - - character(len=512) :: fname_mask - - ! -------------------------------------------------------------------------------------- - - call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR ) - - usage1 = 'USAGE : bin/mkEASETilesParam.x -ease_label EASELabel ' - usage2 = ' where EASELabel = *EASEv[x]_M[yy]*, x={1,2}, yy={01,03,09,25,36}' - - N_args = command_argument_count() - - if(N_args < 1) then + use EASE_conv, only : EASE_extent, EASE_convert, EASE_inverse + use rmTinyCatchParaMod, only : i_raster, j_raster, SRTM_maxcat + use rmTinyCatchParaMod, only : RegridRasterReal + use rmTinyCatchParaMod, only : MAKE_BCS_INPUT_DIR + use process_hres_data, only : histogram + use MAPL_SortMod + use MAPL_ConstantsMod + use MAPL_ExceptionHandling + use netcdf + + implicit none + + integer, parameter :: nc_esa = 129600 ! number of cols in 10-arcsec ESA mask file + integer, parameter :: nr_esa = 64800 ! number of rows in 10-arcsec ESA mask file + + ! define tile types used for processing here (values may be from ESA mask?) + + integer, parameter :: OceanType = 0 + integer, parameter :: LandType = 1 ! land type used for processing here; in GEOS, land tiles are type=100 + integer, parameter :: LakeType = 10 ! lake type used for processing here; in GEOS, lake tiles are type= 19 + integer, parameter :: IceType = 11 ! landice type used for processing here; in GEOS, landice tiles are type= 20 + + real*8, parameter :: Target_mean_land_elev = 614.649D0 + + ! ------------------------------------------------------------ + + integer :: i, j, ig, jg, nn, kkE, kkR, mm + integer :: nc, nr, N_ease_grid_cells + + ! For regridding + + integer, allocatable, dimension(:,:), target :: geos_msk + + REAL, allocatable, DIMENSION(:) :: loc_val + INTEGER, ALLOCATABLE, DIMENSION(:) :: density, loc_int + logical, allocatable, dimension(:) :: unq_mask + integer, dimension(:,:), pointer :: subset + + integer :: dx_esa, dy_esa, NBINS, NPLUS + + integer*8, allocatable, dimension(:) :: SRTM_catid + real(kind=8),allocatable, dimension(:) :: SRTM_catid_r8 + + integer, allocatable, dimension(:,:), target :: tileid_index, catid_index + + ! integer, allocatable, dimension(:,:) :: catid, iaster + + integer, allocatable, dimension(:) :: land_id, water_id, ice_id + integer, allocatable, dimension(:) :: my_land, all_id + real, allocatable, dimension(:) :: ease_grid_area, tile_area + integer*1, allocatable, dimension(:,:) :: veg + real*4, allocatable, dimension(:,:) :: q0, raster + REAL, allocatable, dimension(:) :: tile_elev + + !INTEGER*8 :: PFAF_CODE + + integer :: l, l_index, i_index, w_index, typ, pfaf_index + real :: clat, clon, r_ease, s_ease + real :: fr_gcm + integer :: ind_col, ind_row, status, ncid, varid, nciv + integer :: n_land, n_lake, n_landice, n_landlake, n_landlakelandice + + character*200 :: gfile, gtopo30 + integer :: nc_ease, nr_ease, N_args, command_argument_count + REAL :: dx, dy, d2r, lats, mnx, mxx, mny, mxy, jgv, pix_area + real :: dx_ease, mean_land_elev, total_land_area + character(40) :: arg, EASElabel_ + + character(len=:), allocatable :: EASElabel + + logical :: regrid = .false. + character*128 :: MaskFile + + !logical :: pfaf_til = .false. + + character(len=10) :: nc_string, nr_string + character(len=128) :: usage1, usage2 + character(len=128) :: Iam = "mkEASETilesParam" + character(len=512) :: fname_mask + + ! -------------------------------------------------------------------------------------- + + call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR ) + + usage1 = 'USAGE : bin/mkEASETilesParam.x -ease_label EASELabel ' + usage2 = ' where EASELabel = *EASEv[x]_M[yy]*, x={1,2}, yy={01,03,09,25,36}' + + N_args = command_argument_count() + + if(N_args < 1) then + print *,trim(usage1) + print *,trim(usage2) + stop + end if + + i=0 + do while ( i < N_args ) + + i = i+1 + + call get_command_argument(i,arg) + + if ( trim(arg) == '-ease_label' ) then + i = i+1 + call get_command_argument(i,EASELabel_) + + ! WY noted: this may be used in the future for irrigation tiles + !elseif ( trim(arg) == '-pfaf_til' ) then + ! i = i+1 + ! call get_command_argument(i,PF) + ! if (PF == 'T') pfaf_til = .true. + + else ! stop for any other arguments print *,trim(usage1) print *,trim(usage2) stop - end if - - i=0 - do while ( i < N_args ) - - i = i+1 - - call get_command_argument(i,arg) - - if ( trim(arg) == '-ease_label' ) then - i = i+1 - call get_command_argument(i,EASELabel_) - - ! WY noted: this may be used in the future for irrigation tiles - !elseif ( trim(arg) == '-pfaf_til' ) then - ! i = i+1 - ! call get_command_argument(i,PF) - ! if (PF == 'T') pfaf_til = .true. - - else ! stop for any other arguments - print *,trim(usage1) - print *,trim(usage2) - stop - endif - - end do - - ! Get EASE Grid specifications - ! -------------------------------- - - EASElabel = trim(EASELabel_) - - call ease_extent( EASELabel, nc_ease, nr_ease ) - - write(nc_string, '(i0)') nc_ease - write(nr_string, '(i0)') nr_ease - - gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) - -! ! get appropriate number for length of land_id, water_id, ice_id vectors -! -! if (index(EASELabel,'M01') /=0) then ! EASE 1 km grid -! N_ease_grid_cells = 1500000000 -! else if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid -! N_ease_grid_cells = 1500000000 -! else if (index(EASELabel,'M09') /=0) then ! EASE 9 km grid -! N_ease_grid_cells = 16330000 -! else if (index(EASELabel,'M25') /=0) then ! EASE 25 km grid -! N_ease_grid_cells = 16330000 -! else if (index(EASELabel,'M36') /=0) then ! EASE 36 km grid -! N_ease_grid_cells = 16330000 -! else -! print *,"ERROR: Unknown EASELabel, stopping." -! stop -! endif - - N_ease_grid_cells = nc_ease * nr_ease - - allocate(land_id (1:N_ease_grid_cells)) - allocate(water_id(1:N_ease_grid_cells)) - allocate(ice_id (1:N_ease_grid_cells)) - - land_id = 0 - water_id = 0 - ice_id = 0 - - ! conversion from 2-dim to 1-dim index for EASE grid: l = (ind_row-1)*NDND + ind_col - - !NDND = 10*10**(nint(log10(1.*nr_ease))) - - NDND = nc_ease - - ! Check for the 10 arc-sec MaskFile - ! ----------------------------------- - - call get_environment_variable ("MASKFILE" ,MaskFile ) - - print *, 'Using MaskFile ', trim(MaskFile) - - ! This section was used to make Irrigated Tiles - !if(pfaf_til) then - - ! nc = 43200 ! Number of rows in raster file - ! nr = 21600 - ! call mkEASEv2Raster - - !else - ! if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile - !endif - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - - ! New ESA (Veg) + SRTM (catchments) based mask file - ! is overlaid on the EASE grid - ! ------------------------------------------------- - - ! ESA/SRTM mask file is 10-arcsec resolution; coarsen here to 30-arcsec raster grid: - - ! NOTE: Only coastlines (and lake shorelines?) are at 10-arcsec resolution. - ! The original Pfaf catchments are at 1-arcmin resolution (~2 km). - - nc = 43200 ! Number of rows in raster file - nr = 21600 ! Number of columns in raster file - - regrid = .true. - - dx_esa = nc_esa / nc ! =3 (# of ESA columns within a raster grid cell in x-dim) - dy_esa = nr_esa / nr ! =3 (# of ESA rows within a raster grid cell in y-dim) + endif + + end do + + ! Get EASE Grid specifications + ! -------------------------------- + + EASElabel = trim(EASELabel_) + + call ease_extent( EASELabel, nc_ease, nr_ease ) + + write(nc_string, '(i0)') nc_ease + write(nr_string, '(i0)') nr_ease + + ! assemble output file name (minus file name extension) + + gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) + + N_ease_grid_cells = nc_ease * nr_ease + + allocate(land_id (1:N_ease_grid_cells)) + allocate(water_id(1:N_ease_grid_cells)) + allocate(ice_id (1:N_ease_grid_cells)) + + land_id = 0 + water_id = 0 + ice_id = 0 - allocate(tileid_index(1:nc,1:nr)) - allocate(SRTM_catid (1:SRTM_maxcat+2) ) - allocate(SRTM_catid_r8 (1:SRTM_maxcat+2), source = 0.d0) - allocate(catid_index (1:nc,1:nr)) - allocate(veg (1:nc,1:nr)) - allocate(geos_msk (1:nc_esa,1:dy_esa)) + ! ----------------------------------- + ! + ! Get MaskFile + + call get_environment_variable( "MASKFILE", MaskFile ) + + print *, 'Using MaskFile ', trim(MaskFile) + + ! This section was used to make Irrigated Tiles + !if(pfaf_til) then + + ! nc = 43200 ! Number of rows in raster file + ! nr = 21600 + ! call mkEASEv2Raster + + !else + ! if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile + !endif + + if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then + + ! New ESA (Veg) + SRTM (catchments) based mask file + ! is overlaid on the EASE grid + ! ------------------------------------------------- + + ! ESA/SRTM mask is on 10-arcsec grid + + ! NOTE: Only coastlines and lake shorelines are at 10-arcsec resolution. + ! The original Pfaf catchments are at 1-arcmin resolution (~2 km). + + ! Coarsen 10-arcsec mask to 30-arcsec "raster" grid - ! the following block is not needed (perhaps for routing) - -! allocate(SRTM_CatchArea(1:SRTM_maxcat)) -! -! OPEN (10, FILE = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat', & -! FORM = 'FORMATTED',STATUS='OLD',ACTION='READ') + nc = 43200 ! Number of rows of raster grid + nr = 21600 ! Number of columns of raster grid + + regrid = .true. + + dx_esa = nc_esa / nc ! =3 (# of ESA columns within a raster grid cell in x-dim) + dy_esa = nr_esa / nr ! =3 (# of ESA rows within a raster grid cell in y-dim) + + allocate(tileid_index(1:nc,1:nr)) + allocate(SRTM_catid (1:SRTM_maxcat+2) ) + allocate(SRTM_catid_r8 (1:SRTM_maxcat+2), source = 0.d0) + allocate(catid_index (1:nc,1:nr)) + allocate(veg (1:nc,1:nr)) + allocate(geos_msk (1:nc_esa,1:dy_esa)) + + ! the following block is not needed (perhaps for routing?) + + ! allocate(SRTM_CatchArea(1:SRTM_maxcat)) + ! + ! OPEN (10, FILE = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat', & + ! FORM = 'FORMATTED',STATUS='OLD',ACTION='READ') + ! + ! READ (10,*) I + ! DO N = 1, I + ! READ (10, '(i8,i15,4(1x,f9.4),1x,e10.3,4(1x,e9.3),I8,6(1x,f9.4))') & + ! DOM_INDX,PFAF_CODE,VDUM,VDUM,VDUM,VDUM,VDUM, & + ! SRTM_CatchArea (N) + ! END DO + ! CLOSE (10, STATUS='KEEP') + + dx = 360._8/nc ! raster grid spacing (dlon) + dy = 180._8/nr ! raster grid spacing (dlat) + + d2r = MAPL_PI_R8/180._8 ! degree-to-radians conversion factor -- d2r declared as REAL ?!?!?! + + tileid_index = 0 + catid_index = 0 + veg = OceanType ! initialize to ocean (type=0) + + ! read list of Pfafstetter catchment IDs ('PfafID') from 10-arcsec mask file into variable 'SRTM_catid[_r8]' + ! [vector of length SRTM_maxcat] + + fname_mask = trim(MAKE_BCS_INPUT_DIR) // '/shared/mask/GEOS5_10arcsec_mask.nc' + + print *, 'Opening ', trim(fname_mask) + + status = NF90_OPEN( fname_mask, NF90_NOWRITE, ncid ) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'ERROR -- NF90_OPEN(): ', trim(fname_mask), ' ', ncid, status + print *, trim(Iam), '.x STOPPING.' + stop + else + print *, 'ncid=', ncid + endif + + status = nf90_inq_varid( ncid, name='PfafID', varid=varid ) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'ERROR -- NF90_INQ_VARID(): PfafID ', ncid, varid, status + print *, trim(Iam), '.x STOPPING.' + stop + endif + + status = nf90_get_var( ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/) ) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'ERROR -- NF90_GET_VAR(): PfafID ', ncid, varid, SRTM_maxcat, status + print *, trim(Iam), '.x STOPPING.' + stop + endif + + SRTM_catid = int8(SRTM_catid_r8) ! convert data to integer*8 -- contains 12-digit Pfaf code + SRTM_catid (SRTM_maxcat + 1) = 190000000 ! append ID for Lake type + SRTM_catid (SRTM_maxcat + 2) = 200000000 ! append ID for Landice type + + ! ----------------------------------------------- + ! + ! first loop through 30-arcsec raster grid cells + ! - aggregate mask from 10-arcsec resolution to 30-arcsec raster grid by determining dominant tile type + ! - for each EASE grid cell, determine tile types based on aggregated (30-arcsec) mask + + status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'ERROR -- NF90_INQ_VARID(): CatchIndex ', ncid, varid, status + print *, trim(Iam), '.x STOPPING.' + stop + endif + + do j=1,nr + + clat = -90. + float(j-1)*dy + dy/2. ! center lat of raster grid cell (*,j) + + ! read slice of variable 'CatchIndex' from 10-arcsec mask file into variable 'geos_msk' + ! [2d-array: 129600-by-3] + + status = NF90_GET_VAR( ncid, varid, geos_msk, (/1,(j-1)*dy_esa+1/), (/nc_esa,dy_esa/) ) ! Read 10-arcsec rows that lie within the raster row 'j' + + if(status /=0) then + PRINT *, trim(NF90_STRERROR(STATUS)) + print *, 'ERROR -- NF_GET_VAR(): CatchIndex ', ncid, varid, j, dy_esa, nc_esa, status + print *, trim(Iam), '.x STOPPING.' + stop + endif + + do i = 1,nc + + clon = -180. + float(i-1)*dx + dx/2. ! center lon of raster grid cell (i,*) + + ! extract [3-by-3] subset of ESA/SRTM 10-arcsec mask file that corresponds to 30-arcsec raster grid cell (i,j) + + if (associated (subset)) NULLIFY (subset) + + subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains ESA pixels that lie within the raster grid cell at i,j + + if(maxval (subset) > SRTM_maxcat) then + where (subset == 190000000) subset = SRTM_maxcat + 1 ! Lake type (convert ID from 190000000 to SRTM_maxcat+1) + where (subset == 200000000) subset = SRTM_maxcat + 2 ! Landice type (convert ID from 200000000 to SRTM_maxcat+2) + endif + + if (maxval(subset) > 0) then ! check if there are Non-ocean ESA pixels + + ! raster grid cell has at least one 10-arcsec pixel that is land or landice or lake + ! + ! now find out the *dominant* Pfafstetter ID (could be Lake or Landice) within the raster grid cell + + NPLUS = count( subset>=1 .and. subset<=SRTM_maxcat+2 ) ! Count non-ocean ESA pixels within raster grid cell + allocate(loc_int (1:NPLUS)) + allocate(unq_mask(1:NPLUS)) + loc_int = pack( subset, mask = ( subset>= 1 .and. subset<=SRTM_maxcat+2) ) ! loc_int contains catch_indices of non-ocean ESA pixels + call MAPL_Sort(loc_int) + unq_mask = .true. + do mm=2,NPLUS + unq_mask(mm) = .not. (loc_int(mm)==loc_int(mm-1)) ! count number of unique numbers in loc_int for binning + end do + NBINS = count(unq_mask) + + if (NBINS > 1) then + allocate(loc_val(1:NBINS)) + allocate(density(1:NBINS)) + loc_val = 1.*pack(loc_int,mask =unq_mask) ! loc_val contains available non-ocean catch_indices within the i,j grid cell, + ! Those numbers will be used as bin values + call histogram( dx_esa*dy_esa, NBINS, density, loc_val, real(subset) ) ! density is the pixel count for each bin value + catid_index(i,j) = loc_val(maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j + deallocate(loc_val, density) + else + catid_index(i,j) = loc_int (1) + endif + deallocate(loc_int, unq_mask) + + ! now catid_index(i,j) = dominant PfafID (or LakeID or LandiceID) in raster grid cell (i,j) + + if (catid_index (i,j) <= SRTM_maxcat ) veg(i,j) = LandType + if (catid_index (i,j) == SRTM_maxcat + 1) veg(i,j) = LakeType + if (catid_index (i,j) == SRTM_maxcat + 2) veg(i,j) = IceType + + ! set land_id or water_id or ice_id of EASE grid cell to 1 accordingly + ! --> EASE grid cell may contain more than one tile type + + ! count in if this is i,j pixel is a land, lake or ice within ind_col,ind_row EASE grid cell + + ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) + + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) + + ind_col = nint(r_ease) + 1 + ind_row = nint(s_ease) + 1 ! can be negative or greater than nr_ease (lat near N/S pole) + + if( (veg(i,j).ne.OceanType) .and. (ind_row.ge.1) .and. (ind_row.le.nr_ease) ) then + + ! raster grid cell has tile type other than ocean and is located within the lat band covered + ! by the EASE grid (approx 85S-85N): 1 <= ind_row <= nr_ease + + kkE = (ind_row-1)*nc_ease + ind_col ! 1-dim index for all EASE grid cells + + if( veg(i,j)==LakeType) then + water_id(kkE) = 1 + else if(veg(i,j)==IceType ) then + ice_id (kkE) = 1 + else + land_id (kkE) = 1 + endif + endif + endif + + end do ! i=1,nc + enddo ! j=1,nr + + status = NF90_CLOSE (ncid) + deallocate (geos_msk) + + print *,'Done reading ', trim(MaskFile) + print *,'Min and Max of tile indices:', minval(catid_index), maxval(catid_index) + + else + + print *,'MaskFile = ', trim(MaskFile) + print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' + stop + +! This block processes outdated mask file that was never used with EASE tile space. +! It is kept in the file for reference but has not been maintained after it was commented out +! at the start of the Jun 2023 major cleanup. +! - reichle, 19 Jun 2023 ! -! READ (10,*) I -! DO N = 1, I -! READ (10, '(i8,i15,4(1x,f9.4),1x,e10.3,4(1x,e9.3),I8,6(1x,f9.4))') & -! DOM_INDX,PFAF_CODE,VDUM,VDUM,VDUM,VDUM,VDUM, & -! SRTM_CatchArea (N) -! END DO -! CLOSE (10, STATUS='KEEP') - - dx = 360._8/nc ! raster grid spacing (dlon) - dy = 180._8/nr ! raster grid spacing (dlat) - - d2r = MAPL_PI_R8/180._8 ! degree-to-radians conversion factor -- d2r declared as REAL ?!?!?! - - !da = MAPL_radius*MAPL_radius*pi*pi*dx*dy/180./180./1000000. - - tileid_index = 0 - catid_index = 0 - veg = OceanType ! initialize to ocean (type=0) - - ! read list of Pfafstetter catchment IDs ('PfafID') from 10-arcsec mask file into variable 'SRTM_catid[_r8]' - ! [vector of length SRTM_maxcat] - - fname_mask = trim(MAKE_BCS_INPUT_DIR) // '/shared/mask/GEOS5_10arcsec_mask.nc' - - print *, 'Opening ', trim(fname_mask) - - status = NF90_OPEN( fname_mask, NF90_NOWRITE, ncid ) - if(status /=0) then - PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'ERROR -- NF90_OPEN(): ', trim(fname_mask), ' ', ncid, status - print *, trim(Iam), '.x STOPPING.' - stop - else - print *, 'ncid=', ncid - endif - - status = nf90_inq_varid( ncid, name='PfafID', varid=varid ) - if(status /=0) then - PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'ERROR -- NF90_INQ_VARID(): PfafID ', ncid, varid, status - print *, trim(Iam), '.x STOPPING.' - stop - endif - - status = nf90_get_var( ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/) ) - if(status /=0) then - PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'ERROR -- NF90_GET_VAR(): PfafID ', ncid, varid, SRTM_maxcat, status - print *, trim(Iam), '.x STOPPING.' - stop - endif - - SRTM_catid = int8(SRTM_catid_r8) ! convert data to integer*8 -- contains 12-digit Pfaf code - SRTM_catid (SRTM_maxcat + 1) = 190000000 ! append ID for Lake type - SRTM_catid (SRTM_maxcat + 2) = 200000000 ! append ID for Landice type - - ! ----------------------------------------------- - ! - ! first loop through 30-arcsec raster grid cells - ! - aggregate mask from 10-arcsec resolution to 30-arcsec raster grid by determining dominant tile type - ! - for each EASE grid cell, determine tile types based on aggregated (30-arcsec) mask - - !i1 = 0 ! count # of 30-arcsec pixels - NOT USED - - status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) - if(status /=0) then - PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'ERROR -- NF90_INQ_VARID(): CatchIndex ', ncid, varid, status - print *, trim(Iam), '.x STOPPING.' - stop - endif - - do j=1,nr - - clat = -90. + float(j-1)*dy + dy/2. ! center lat of raster grid cell (*,j) - - ! read slice of variable 'CatchIndex' from 10-arcsec mask file into variable 'geos_msk' - ! [2d-array: 129600-by-3] - - status = NF90_GET_VAR( ncid, varid, geos_msk, (/1,(j-1)*dy_esa+1/), (/nc_esa,dy_esa/) ) ! Read 10-arcsec rows that lie within the raster row 'j' - !status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' - - if(status /=0) then - PRINT *, trim(NF90_STRERROR(STATUS)) - print *, 'ERROR -- NF_GET_VAR(): CatchIndex ', ncid, varid, j, dy_esa, nc_esa, status - print *, trim(Iam), '.x STOPPING.' - stop - endif - - do i = 1,nc - - clon = -180. + float(i-1)*dx + dx/2. ! center lon of raster grid cell (i,*) - - ! extract [3-by-3] subset of ESA/SRTM 10-arcsec mask file that corresponds to 30-arcsec raster grid cell (i,j) - - if (associated (subset)) NULLIFY (subset) - - subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains ESA pixels that lie within the raster grid cell at i,j - - if(maxval (subset) > SRTM_maxcat) then - where (subset == 190000000) subset = SRTM_maxcat + 1 ! Lake type (convert ID from 190000000 to SRTM_maxcat+1) - where (subset == 200000000) subset = SRTM_maxcat + 2 ! Landice type (convert ID from 200000000 to SRTM_maxcat+2) - endif - - if (maxval(subset) > 0) then ! check if there are Non-ocean ESA pixels - - ! raster grid cell has at least one 10-arcsec pixel that is land or landice or lake - ! - ! now find out the *dominant* Pfafstetter ID (could be Lake or Landice) within the raster grid cell - - NPLUS = count( subset>=1 .and. subset<=SRTM_maxcat+2 ) ! Count non-ocean ESA pixels within raster grid cell - allocate(loc_int (1:NPLUS)) - allocate(unq_mask(1:NPLUS)) - loc_int = pack( subset, mask = ( subset>= 1 .and. subset<=SRTM_maxcat+2) ) ! loc_int contains catch_indices of non-ocean ESA pixels - call MAPL_Sort(loc_int) - unq_mask = .true. - do n=2,NPLUS - unq_mask(n) = .not. (loc_int(n)==loc_int(n-1)) ! count number of unique numbers in loc_int for binning - end do - NBINS = count(unq_mask) - - if (NBINS > 1) then - allocate(loc_val(1:NBINS)) - allocate(density(1:NBINS)) - loc_val = 1.*pack(loc_int,mask =unq_mask) ! loc_val contains available non-ocean catch_indices within the i,j grid cell, - ! Those numbers will be used as bin values - call histogram( dx_esa*dy_esa, NBINS, density, loc_val, real(subset) ) ! density is the pixel count for each bin value - catid_index(i,j) = loc_val(maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j - deallocate(loc_val, density) - else - catid_index(i,j) = loc_int (1) - endif - deallocate(loc_int, unq_mask) - - ! now catid_index(i,j) = dominant PfafID (or LakeID or LandiceID) in raster grid cell (i,j) - - if (catid_index (i,j) <= SRTM_maxcat ) veg(i,j) = LandType - if (catid_index (i,j) == SRTM_maxcat + 1) veg(i,j) = LakeType - if (catid_index (i,j) == SRTM_maxcat + 2) veg(i,j) = IceType - - !if( (catid_index(i,j)>= 1) .and. (catid_index(i,j)<=SRTM_maxcat) ) i1 = i1 + 1 ! i1 = counter for raster grid cells with land as dominant type -- NOT USED? - - ! set land_id or water_id or ice_id of EASE grid cell to 1 accordingly - ! --> EASE grid cell may contain more than one tile type - - ! count in if this is i,j pixel is a land, lake or ice within ind_col,ind_row EASE grid cell - - ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) - - call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) - - ind_col = nint(r_ease) + 1 - ind_row = nint(s_ease) + 1 ! can be negative or greater than nr_ease (lat near N/S pole) - - if( (veg(i,j).ne.OceanType) .and. (ind_row.ge.1) .and. (ind_row.le.nr_ease) ) then - - ! raster grid cell has tile type other than ocean and is located within the lat band covered - ! by the EASE grid (approx 85S-85N): 1 <= ind_row <= nr_ease - - ! 1-dim indexing used here appears wasteful; why not ll=(ind_row-1)*nc_ease+ind_col, with NT=nc_ease*nr_ease ?? - - !l = ind_row*NDND + ind_col - - l = (ind_row-1)*NDND + ind_col ! 1-dim index for all EASE grid cells - - if( veg(i,j)==LakeType) then - water_id(l) = 1 - else if(veg(i,j)==IceType ) then - ice_id (l) = 1 - else - land_id (l) = 1 - endif - endif - endif - - end do ! i=1,nc - enddo ! j=1,nr - - status = NF90_CLOSE (ncid) - deallocate (geos_msk) - - print *,'Done reading ', trim(MaskFile) - print *,'Min and Max of tile indices:', minval(catid_index), maxval(catid_index) - - else - - print *,'MaskFile = ', trim(MaskFile) - print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' - stop - ! ! Old IGBP (Veg) + HYDRO1k (catchments) based mask will ! ! Overlaid on EASE mask ! ! ----------------------------------------------------- @@ -586,400 +560,400 @@ PROGRAM mkEASETilesParam ! end do - endif ! (GEOS5_10arcsec_mask) - - ! -------------------------------------------------------- - ! - ! Read SRTM elevation data - to be consistent with AGCM - - allocate(raster(i_raster,j_raster)) ! 2.5-min raster grid ( 8640-by-4320 ) - allocate(q0 (nc, nr) ) ! 30-arcsec raster grid (43200-by-21600) - - gtopo30 = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/srtm30_withKMS_2.5x2.5min.data' + endif ! (GEOS5_10arcsec_mask) + + ! -------------------------------------------------------- + ! + ! Read SRTM elevation data - to be consistent with AGCM + + allocate(raster(i_raster,j_raster)) ! 2.5-min raster grid ( 8640-by-4320 ) + allocate(q0 (nc, nr) ) ! 30-arcsec raster grid (43200-by-21600) + + gtopo30 = trim(MAKE_BCS_INPUT_DIR)//'/land/topo/v1/srtm30_withKMS_2.5x2.5min.data' + + open (10,file=trim(gtopo30),form='unformatted',status='old',convert='little_endian') + read (10) raster + close (10,status='keep') + + ! remap SRTM elevation data from 2.5-min to 30-arcsec raster grid + + if(regrid) then + call RegridRasterReal(raster,q0) + else + q0 = raster + endif + + deallocate (raster) + + ! --------------------------------------------------------- + ! + ! determine number of land, lake, and landice tiles + + n_land = sum(land_id) ! number of land tiles + n_lake = sum(water_id) ! number of lake tiles + n_landice = sum(ice_id) ! number of landice tiles + + n_landlake = n_land + n_lake + n_landlakelandice = n_landlake + n_landice + + print *,'# of Land tiles: ', n_land + print *,'# of Lake tiles: ', n_lake + print *,'# of Landice tiles: ', n_landice + print *,'# of Land+Lake tiles: ', n_landlake + print *,'# of Land+Lake+Landice tiles: ', n_landlakelandice + + l_index = 0 ! initialize index for land tiles + w_index = n_land ! initialize index for lake tiles + i_index = n_landlake ! initialize index for landice tiles + + allocate(ease_grid_area(1:N_ease_grid_cells)) + + allocate(tile_area (1:n_landlakelandice)) + allocate(my_land (1:n_landlakelandice)) + allocate(all_id (1:n_landlakelandice)) + + allocate(tile_elev (1:n_land)) + + ! =========================================================================== + ! + ! prepare for second loop through raster grid cells + + land_id = 0 + water_id = 0 + ice_id = 0 + + my_land = 0 + all_id = 0 + + ease_grid_area = 0. + tile_elev = 0. + tile_area = 0. + + ! The second loop through the nc x nr raster grid derives land, ice and water tiles. + ! Each EASE grid cell is assigned with an ID [1-dim index] kkE = (ind_row-1)*nc_ease + ind_col + ! where ind_col, ind_row are overlying EASE grid cell indices + ! Based on the above sums: + ! n_land EASE grid cells have non-zero land fraction (sum(land_id )) + ! n_lake EASE grid cells have non-zero inland water fractions (sum(water_id)) + ! n_landice EASE grid cells have non-zero landice fractions (sum(ice_id )) + ! hence, Tiles with index 1 to n_land represent land tiles + ! Tiles with index n_land+1 to n_landlake represent lake tiles + ! Tiles with index n_landlake+1 to n_landlakelandice represent landice tiles + ! The global nc x nr array of tileid_index(nc,nr) contains corresponding tile_index values. + + ! Conversion from 2-dim index to 1-dim index for *raster* grid: all_id(nn) = (j-1)*nc + i + + do i = 1 ,nc - open (10,file=trim(gtopo30),form='unformatted',status='old',convert='little_endian') - read (10) raster - close (10,status='keep') - - ! remap SRTM elevation data from 2.5-min to 30-arcsec raster grid - - if(regrid) then - call RegridRasterReal(raster,q0) - else - q0 = raster - endif - - deallocate (raster) - - ! --------------------------------------------------------- - ! - ! determine number of land, lake, and landice tiles - - n_land = sum(land_id) ! number of land tiles - n_lake = sum(water_id) ! number of lake tiles - n_landice = sum(ice_id) ! number of landice tiles - - n_landlakelandice = n_land + n_lake + n_landice - - print *,'# of Land tiles: ', n_land - print *,'# of Lake tiles: ', n_lake - print *,'# of Landice tiles: ', n_landice - print *,'# of Land+Lake+Landice tiles: ', n_landlakelandice - - l_index = 0 ! start index for land EASE cells - w_index = n_land ! start index for lake EASE cells - i_index = w_index + n_lake ! start index for landice EASE cells - - allocate(ease_grid_area(1:N_ease_grid_cells)) - - allocate(tile_area (1:n_landlakelandice)) - allocate(my_land (1:n_landlakelandice)) - allocate(all_id (1:n_landlakelandice)) - - allocate(tile_ele (1:n_land)) - - ! =========================================================================== - ! - ! prepare for second loop through raster grid cells - - land_id = 0 - water_id = 0 - ice_id = 0 - - my_land = 0 - all_id = 0 - - ease_grid_area = 0. - tile_ele = 0. - tile_area = 0. - - ! While looping through the nc x nr grid, this section derives land, ice and water tiles. - ! Each EASE grid cell is assigned with an ID = ind_row*ND + ind_col - ! ind_col, ind_row are overlying EASE grid cell indices - ! Based on the above sums: - ! l_index Grid cells have land fractions (sum(land_id)) - ! w_index EASE Grid cells have inland water fractions (sum(water_id)) - ! i_index EASE Grid cells have ice fractions (sum(ice_id)) - ! hence, tile_index 1 to l_index represent land tiles - ! tile_index l_index +1 to l_index + w_index represent water (lakes) tiles - ! tile_index l_index + w_index +1 to l_index + w_index + i_index represent ice tiles - ! global nc x nr array of tileid_index(nc,nr) contains corresponding tile_index values which - ! is derived in the below loop - - ! conversion from 2-dim index to 1-dim index for raster grid: all_ind(l) = (j-1)*ND_raster + i - - !ND_raster = 10*10**(nint(log10(1.*NR))) - - ND_raster = nc - - !i2 = 1 -- NOT USED - - do i = 1 ,nc - - clon = -180. + float(i-1)*dx + dx/2. ! center lon of raster grid cell (*,j) - - do j =nr ,1 ,-1 - - lats = -90._8 + (j - 0.5_8)*dy ! center lat of raster grid cell (*,j) -- lats declared REAL ?!?!?! same as clat ?!?!?! - clat = -90. + float(j-1)*dy + dy/2. ! center lat of raster grid cell (*,j) - - ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) - - call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) - - ind_col = nint(r_ease) + 1 - ind_row = nint(s_ease) + 1 ! can be negative or greater than nr_ease (lat near N/S pole) - - if( (ind_row.ge.1) .and. (ind_row.le.nr_ease) ) then - - ! raster grid cell (i,j) is located within the lat band covered by the (cylindrical) EASE grid (approx 85S-85N) - - !l= ind_row*NDND + ind_col ! 1-dim index for EASE grid cells - - l = (ind_row-1)*NDND + ind_col ! 1-dim index for all EASE grid cells - - pix_area = ( sin(d2r*(lats+0.5*dy)) - sin(d2r*(lats-0.5*dy)) )*(dx*d2r) ! area of 30-arcsec raster grid cell - - if (veg(i,j).ne.OceanType) then - - ! set tile ID and compute tile area - - select case (veg(i,j)) - - case (LakeType) ! raster grid cell (i,j) is Lake - - if(water_id(l)==0) then - ! raster grid cell is first Lake type seen for this EASE grid cell, tile ID has not yet been set - ! recall: tile IDs for Lake: [(n_land+1):(n_land+n_lake)], and w_index was initalized to n_land above - w_index = w_index + 1 - water_id(l) = w_index ! needed so if condition will be false for next raster grid cell of same type - end if - - tileid_index(i,j) = water_id(l) - - case (IceType) ! raster grid cell (i,j) is Landice - - if(ice_id(l)==0) then - ! raster grid cell is first Landice type seen for this EASE grid cell, tile ID has not yet been set - ! recall: tile IDs for Landice: [(n_land+n_lake+1):n_landlakelandice], and i_index was initalized to n_land+n_lake above - i_index = i_index + 1 - ice_id(l) = i_index ! needed so if condition will be false for next raster grid cell of same type - end if - - tileid_index(i,j) = ice_id(l) - - case (LandType) ! raster grid cell (i,j) is Land - - if(land_id(l)==0) then - ! raster grid cell is first Land type seen for this EASE grid cell, tile ID has not yet been set - ! recall: tile IDs for Land: [1:n_land], and l_index was initalized to 0 above - l_index = l_index + 1 - land_id(l) = l_index ! needed so if condition will be false for next raster grid cell of same type - end if - - tileid_index(i,j) = land_id(l) - - ! sum up area and (area-weighted) elevation (only over raster grid cells of type land!) - - tile_ele( tileid_index(i,j)) = tile_ele( tileid_index(i,j)) + q0(i,j) * pix_area ! q0 = elevation - - ! tile_area_land should be obsolete because identical to tile_area(1:n_land) - !tile_area_land(tileid_index(i,j)) = tile_area_land(tileid_index(i,j)) + pix_area - - case default - - print *,'ERROR: unknown tile type value in veg(i,j): ', veg(i,j), ' STOPPING.' - stop - - end select - - ! sum up area of raster grid cells contributing to each tile (land or lake or landice) - - tile_area(tileid_index(i,j)) = tile_area(tileid_index(i,j)) + pix_area - - my_land(tileid_index(i,j)) = l ! for this tile, store 1-dim index for EASE grid cells - all_id( tileid_index(i,j)) = (j-1)*ND_raster + i ! for this tile, store 1-dim index for raster grid cells - last in prevails! - - ! BUG??? all_id stores only the 1-dim index of the 30-arcsec raster grid cells that is the - ! last one to contribute to the EASE tile specified by tileid_index(i,j) - ! This does not seem to be the desired dominant ID across the EASE grid cell - - endif ! raster grid cell has tile type other than ocean - - ! compute total area of raster grid cells of any tile type that contribute to EASE grid cell - - ease_grid_area(l) = ease_grid_area(l) + pix_area - - endif ! raster grid cell is located within lat band of EASE grid - - end do - end do - - ! verify l_index, w_index, i_index against tile counts from first loop - - if (l_index/= n_land ) then - print *, 'ERROR: l_index = ', l_index, ' must now match # of Land tiles = ', n_land - print *, 'STOPPING.' - stop - end if - - if (w_index/=(n_land+n_lake) ) then - print *, 'ERROR: w_index = ', w_index, ' must now match # of Land+Lake tiles = ', n_land+n_lake - print *, 'STOPPING.' - stop - - end if - - if (i_index/= n_landlakelandice ) then - print *, 'ERROR: i_index = ', i_index, ' must now match # of Land+Lake+Landice tiles = ', n_landlakelandice - print *, 'STOPPING.' - stop - - end if - - deallocate(land_id, q0) - deallocate(water_id) - deallocate(ice_id ) - - tile_ele = tile_ele/tile_area(1:n_land) ! finalize tile elevation - - ! adjustment Global Mean Topography to 614.649 (615.662 GTOPO 30) m - ! ----------------------------------------------------------------- - mean_land_ele=0. - - do j=1,l_index - mean_land_ele = mean_land_ele + tile_ele(j)*tile_area(j) - enddo - - mean_land_ele = mean_land_ele/sum(tile_area(1:l_index)) - - if(mean_land_ele .ne. 614.649D0 ) then - print *,'Global Mean Elevation (over land): ', mean_land_ele - tile_ele = tile_ele*(614.649D0/mean_land_ele) - ! verify that adjust elevation is correct - mean_land_ele=0. - do j=1,l_index - mean_land_ele = mean_land_ele + tile_ele(j)*tile_area(j) - enddo - print *,'Global Mean Elevation after scaling to SRTM : ',mean_land_ele/sum(tile_area(1:l_index)) - endif - print *,'Total Land Area :', sum(tile_area(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. ! , & - !sum(tile_area_land(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. - - !------------------------------------------- - ! - ! write *.rst - - print *,'Writing ... ', trim(gfile)//'.rst' - - open (10, file ='rst/'//trim(gfile)//'.rst',form='unformatted',status='unknown', & - action='write') - - do j=1,nr - write(10)(tileid_index(i,j),i=1,nc) - end do - - close (10,status='keep') - - !----------------------------------------------------------- - - ! write catchment.def and *.til files - - print *,'Writing ... ', trim(gfile)//'.til and catchment.def' - - open (11,file='clsm/catchment.def', & - form='formatted',status='unknown') - write(11,*)l_index - - open (10, file ='til/'//trim(gfile)//'.til',form='formatted',status='unknown',action='write') - write (10,*) i_index,SRTM_maxcat, nc, nr - write (10,*) 1 - write (10,*) EASELabel - write (10,*) nc_ease - write (10,*) nr_ease - - ! write (10,*)'NO-OCEAN' - ! write (10,*) -9999 - ! write (10,*) -9999 - - dx_ease = 180./real(nc_ease) + clon = -180. + float(i-1)*dx + dx/2. ! center lon of raster grid cell (*,j) + + do j =nr ,1 ,-1 + + lats = -90._8 + (j - 0.5_8)*dy ! center lat of raster grid cell (*,j) -- lats declared REAL ?!?!?! same as clat ?!?!?! + clat = -90. + float(j-1)*dy + dy/2. ! center lat of raster grid cell (*,j) + + ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) + + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) + + ind_col = nint(r_ease) + 1 + ind_row = nint(s_ease) + 1 ! NOTE: can be zero or negative or greater than nr_ease (lat near N/S pole) + + if( (ind_row.ge.1) .and. (ind_row.le.nr_ease) ) then + + ! raster grid cell (i,j) is located within the lat band covered by the (cylindrical) EASE grid (approx 85S-85N) + + kkE = (ind_row-1)*nc_ease + ind_col ! 1-dim index for EASE grid cells + + kkR = (j -1)*nc + i ! 1-dim index for raster grid cells + + pix_area = ( sin(d2r*(lats+0.5*dy)) - sin(d2r*(lats-0.5*dy)) )*(dx*d2r) ! area of 30-arcsec raster grid cell + + if (veg(i,j).ne.OceanType) then + + ! set tile ID and compute tile area + + select case (veg(i,j)) + + case (LakeType) ! raster grid cell (i,j) is Lake + + if(water_id(kkE)==0) then + ! raster grid cell is first Lake type seen for this EASE grid cell, tile ID has not yet been set + ! recall: tile IDs for Lake: [(n_land+1):(n_land+n_lake)], and w_index was initalized to n_land above + w_index = w_index + 1 + water_id(kkE) = w_index ! needed so if condition will be false for next raster grid cell of same type + end if + + tileid_index(i,j) = water_id(kkE) + + case (IceType) ! raster grid cell (i,j) is Landice + + if(ice_id(kkE)==0) then + ! raster grid cell is first Landice type seen for this EASE grid cell, tile ID has not yet been set + ! recall: tile IDs for Landice: [(n_land+n_lake+1):n_landlakelandice], and i_index was initalized to n_land+n_lake above + i_index = i_index + 1 + ice_id(kkE) = i_index ! needed so if condition will be false for next raster grid cell of same type + end if + + tileid_index(i,j) = ice_id(kkE) + + case (LandType) ! raster grid cell (i,j) is Land + + if(land_id(kkE)==0) then + ! raster grid cell is first Land type seen for this EASE grid cell, tile ID has not yet been set + ! recall: tile IDs for Land: [1:n_land], and l_index was initalized to 0 above + l_index = l_index + 1 + land_id(kkE) = l_index ! needed so if condition will be false for next raster grid cell of same type + end if + + tileid_index(i,j) = land_id(kkE) + + ! sum up area and (area-weighted) elevation (only over raster grid cells of type land!) + + tile_elev( tileid_index(i,j)) = tile_elev( tileid_index(i,j)) + q0(i,j) * pix_area ! q0 = elevation + + ! tile_area_land should be obsolete because identical to tile_area(1:n_land) + !tile_area_land(tileid_index(i,j)) = tile_area_land(tileid_index(i,j)) + pix_area + + case default + + print *,'ERROR: unknown tile type value in veg(i,j): ', veg(i,j), ' STOPPING.' + stop + + end select + + ! sum up area of raster grid cells contributing to each tile (land or lake or landice) + + tile_area(tileid_index(i,j)) = tile_area(tileid_index(i,j)) + pix_area + + my_land(tileid_index(i,j)) = kkE ! for this tile, store 1-dim index for EASE grid cells + all_id( tileid_index(i,j)) = kkR ! for this tile, store 1-dim index for raster grid cells - last in prevails! + + ! BUG??? all_id stores only the 1-dim index of the 30-arcsec raster grid cells that is the + ! last one to contribute to the EASE tile specified by tileid_index(i,j) + ! This does not seem to be the desired dominant ID across the EASE grid cell + + endif ! raster grid cell has tile type other than ocean + + ! compute total area of raster grid cells of any tile type that contribute to EASE grid cell + + ease_grid_area(kkE) = ease_grid_area(kkE) + pix_area + + endif ! raster grid cell is located within lat band of EASE grid + + end do + end do + + ! verify l_index, w_index, i_index against tile counts from first loop + + if (l_index/= n_land ) then + print *, 'ERROR: l_index = ', l_index, ' must now match # of Land tiles = ', n_land + print *, 'STOPPING.' + stop + end if + + if (w_index/=(n_landlake) ) then + print *, 'ERROR: w_index = ', w_index, ' must now match # of Land+Lake tiles = ', n_landlake + print *, 'STOPPING.' + stop + + end if + + if (i_index/= n_landlakelandice ) then + print *, 'ERROR: i_index = ', i_index, ' must now match # of Land+Lake+Landice tiles = ', n_landlakelandice + print *, 'STOPPING.' + stop + + end if + + deallocate(land_id, q0) + deallocate(water_id) + deallocate(ice_id ) + + tile_elev = tile_elev/tile_area(1:n_land) ! finalize tile elevation + - do l=1,i_index + ! --------------------------------------------------------------------------------- + ! + ! adjust tile elevation to match Target_mean_land_elev (parameter defined above) + ! + ! mean land(?) elevation in GTOPO30 = 615.662 m + + total_land_area = sum(tile_area(1:n_land)) ! in radians squared - !ig = my_land(l)-NDND*(my_land(l)/NDND) - !jg = my_land(l)/NDND + mean_land_elev=0. - ! for EASE grid: get row and column indices from 1-dim index (that is, invert l=(ind_row-1)*NDND+ind_col) - - jg = (my_land(l)-1)/NDND + 1 ! = ind_row (1-based) [note integer division] - - ig = my_land(l) - NDND*(jg-1) ! = ind_col (1-based) - - ! extract original PfafID from catid_index - ! - ! BUG??? catid_index is in 30-arcsec raster space and contains the dominant PfafID - ! (or LakeID, or LandiceID) within the 30-arcsec raster grid cell, - ! where dominant is w.r.t. the 10-arcsec mask with PfafIDs - ! It seems that here we just pick the PfafID associated with a single - ! 30-arcsec raster grid cell within the EASE tile, which happens to be the - ! last of the 30-arcsec raster grid cells that contribute to the EASE tile - ! - - ! for raster grid: get row and column indices from 1-dim index (that is, invert all_id(l)=(j-1)*ND_raster + i + do j=1,n_land + mean_land_elev = mean_land_elev + tile_elev(j)*tile_area(j) + enddo + + mean_land_elev = mean_land_elev/total_land_area + + if(mean_land_elev .ne. Target_mean_land_elev ) then - j = (all_id(l) -1)/ND_raster + 1 - - i = all_id(l) - ND_raster*(j-1) - - !cindex= catid_index(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) - - cindex = catid_index(i,j) - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - pfaf = cindex - else - print *,'MaskFile = ', trim(MaskFile) - print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' - stop - ! pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) - endif - - if ((l>l_index) .and. (l<=w_index)) typ = 19 ! Lake tile - - if ( l>w_index ) typ = 20 ! Landice tile - - if (l<=l_index) then + print *, 'Global mean land elevation before adjustment [m]: ', mean_land_elev - typ = 100 ! Land tile + tile_elev = tile_elev*(Target_mean_land_elev/mean_land_elev) - ! get min/max lat/lon of EASE grid cell - ! BUG: This is *not* the desired min/max lat/lon of the land tile!!! - - call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) - - mnx = clon - dx_ease - mxx = clon + dx_ease - - jgv = real(jg-1) + 0.5 - - call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) + ! print adjusted elevation to log file + mean_land_elev=0. + do j=1,n_land + mean_land_elev = mean_land_elev + tile_elev(j)*tile_area(j) + enddo + print *, 'Global mean land elevation after scaling to SRTM [m]: ', mean_land_elev/total_land_area - mny = clat - - jgv = real(jg-1) - 0.5 - - call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) + endif - mxy = clat - - ! write tile properties into catchment.def file + print *, 'Total land area [km^2]: ', total_land_area*MAPL_RADIUS*MAPL_RADIUS/1000./1000. + + !------------------------------------------- + ! + ! write *.rst file + + print *, 'Writing ... ', trim(gfile)//'.rst' + + open( 10, file ='rst/'//trim(gfile)//'.rst', form='unformatted', status='unknown', action='write' ) + + do j=1,nr + write(10) (tileid_index(i,j),i=1,nc) + end do + + close (10,status='keep') + + !----------------------------------------------------------- + + ! write catchment.def and *.til files + + print *, 'Writing ... ', trim(gfile)//'.til and catchment.def' + + open( 11, file='clsm/catchment.def', form='formatted', status='unknown' ) - write (11,'(i10,i8,5(2x,f9.4), i4)') l, pfaf, mnx, mxx, mny, mxy, tile_ele(l) + write( 11,*) n_land + + open( 10, file ='til/'//trim(gfile)//'.til', form='formatted', status='unknown', action='write' ) - endif + write( 10,*) n_landlakelandice, SRTM_maxcat, nc, nr + write( 10,*) 1 + write( 10,*) EASELabel + write( 10,*) nc_ease + write( 10,*) nr_ease - ! get area fraction of tile within EASE grid cell - ! NOTE: the area of the EASE grid cell here has to be the sum of the areas of the - ! contributing raster grid cells, which is *not* the same for all EASE grid cells - - call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) - - fr_gcm = tile_area(l) / ease_grid_area((jg-1)*NDND+ig) - - ! write tile properties into *.til file - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - - ! Note: If-condition uses MaskFile as a proxy for format & specs of *.til file (columns written) - - ! pfaf = running *index* of Pfafstetter catchment (1:SRTM_maxcat) - ! SRTM_catid = 12-digit Pfafstetter code (encoding network/routing info) - - write(10,'(i10,i9,2f10.4,2i6,f19.12,i10,i15,e13.4)') & - typ, pfaf, clon, clat, ig-1, jg-1, fr_gcm, pfaf, SRTM_catid(cindex) - else + ! obsolete header lines: + ! + ! write (10,*)'NO-OCEAN' + ! write (10,*) -9999 + ! write (10,*) -9999 + + dx_ease = 180./real(nc_ease) + + do nn=1,n_landlakelandice + + ! get *EASE* grid row and column indices from 1-dim index (that is, invert kkE=(ind_row-1)*nc_ease+ind_col) - ! write statement below was used with MaskFile prior to availability of SRTM-based Pfafstetter catchments - ! obsolete with EASE grid - ! - reichle, 15 Jun 2023 + kkE = my_land(nn) + + jg = ((kkE-1)/nc_ease) + 1 ! = ind_row (1-based) [note integer division] + + ig = kkE - nc_ease*(jg-1) ! = ind_col (1-based) + + ! extract original PfafID from catid_index + ! + ! BUG??? catid_index is in 30-arcsec raster space and contains the dominant PfafID + ! (or LakeID, or LandiceID) within the 30-arcsec raster grid cell, + ! where dominant is w.r.t. the 10-arcsec mask with PfafIDs + ! It seems that here we just pick the PfafID associated with a single + ! 30-arcsec raster grid cell within the EASE tile, which happens to be the + ! last of the 30-arcsec raster grid cells that contribute to the EASE tile + ! + + ! get *raster* grid i and j indices from 1-dim index (that is, invert kkR=(j-1)*nc + i + + kkR = all_id(nn) - print *,'MaskFile = ', trim(MaskFile) - print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' - stop -! -! write(10,'(i10,i9,2f10.4,2i5,f19.12,i10,e13.4,i8)') & + j = ((kkR - 1)/nc) + 1 ! (1-based) [note integer division] + + i = kkR - nc*(j-1) ! (1-based) + + pfaf_index = catid_index(i,j) + + if ((nn>n_land) .and. (nn<=n_landlake)) typ = 19 ! Lake tile + + if ( nn>n_landlake ) typ = 20 ! Landice tile + + if ( nn<=n_land) then + + typ = 100 ! Land tile + + ! get min/max lat/lon of EASE grid cell + ! BUG: This is *not* the desired min/max lat/lon of the land tile!!! + + call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) + + mnx = clon - dx_ease + mxx = clon + dx_ease + + jgv = real(jg-1) + 0.5 + + call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) + + mny = clat + + jgv = real(jg-1) - 0.5 + + call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) + + mxy = clat + + ! write tile properties into catchment.def file + + write (11,'(i10,i8,5(2x,f9.4), i4)') nn, pfaf_index, mnx, mxx, mny, mxy, tile_elev(nn) + + endif + + ! get area fraction of tile within EASE grid cell + ! + ! NOTE: The area of the EASE grid cell here has to be the sum of the areas of the + ! contributing raster grid cells, which is *not* the same for all EASE grid cells; + ! that is, cannot use exact (globally constant) area of EASE grid cell. + + call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) + + fr_gcm = tile_area(nn) / ease_grid_area((jg-1)*nc_ease+ig) + + ! write tile properties into *.til file + + if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then + + ! Note: If-condition uses MaskFile as a proxy for format & specs of *.til file (columns written) + + ! pfaf_index = running *index* of Pfafstetter catchment (1:SRTM_maxcat) + ! SRTM_catid = 12-digit Pfafstetter code (encoding network/routing info) + + write(10,'(i10,i9,2f10.4,2i6,f19.12,i10,i15,e13.4)') & + typ, pfaf_index, clon, clat, ig-1, jg-1, fr_gcm, pfaf_index, SRTM_catid(pfaf_index) + else + + ! write statement below was used with MaskFile prior to availability of SRTM-based Pfafstetter catchments + ! obsolete with EASE grid + ! - reichle, 15 Jun 2023 + + print *,'MaskFile = ', trim(MaskFile) + print *,'ERROR: Selected mask file not supported for creating tiles on the EASE grid, stopping.' + stop + ! + ! write(10,'(i10,i9,2f10.4,2i5,f19.12,i10,e13.4,i8)') & ! typ,pfaf,clon,clat,ig-1,jg-1,fr_gcm ,cindex - endif - - end do - - close(10,status='keep') - close(11,status='keep') - - deallocate( tileid_index, catid_index,veg ) - deallocate( tile_area, ease_grid_area, tile_ele, my_land, all_id ) + endif + + end do + + close(10,status='keep') + close(11,status='keep') + + deallocate( tileid_index, catid_index,veg ) + deallocate( tile_area, ease_grid_area, tile_elev, my_land, all_id ) - ! Commented out "empty" if-block. -rreichle, 15 Jun 2023 + ! Commented out "empty" if-block. -rreichle, 15 Jun 2023 ! ! if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then ! @@ -991,24 +965,24 @@ PROGRAM mkEASETilesParam ! ! endif - ! create Grid2Catch transfer file - ! ------------------------------- - - ! CALL CREATE_ROUT_PARA_FILE (NC, NR, trim(gfile), MGRID=MGRID) - - ! now run mkCatchParam - ! -------------------- - - ! WY Note: now mkCatchParam is run in the make_bcs script, not here - ! and nthread will be reset to run mkCatchParam - - ! tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBCSV) - ! write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr - ! tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1) - ! print *,trim(tmpstring) - - ! call execute_command_line (tmpstring) - + ! create Grid2Catch transfer file + ! ------------------------------- + + ! CALL CREATE_ROUT_PARA_FILE (NC, NR, trim(gfile), MGRID=MGRID) + + ! now run mkCatchParam + ! -------------------- + + ! WY Note: now mkCatchParam is run in the make_bcs script, not here + ! and nthread will be reset to run mkCatchParam + + ! tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBCSV) + ! write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr + ! tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1) + ! print *,trim(tmpstring) + + ! call execute_command_line (tmpstring) + !!! commented out. It may be used in the future for irrigation tiles !!! contains @@ -1091,5 +1065,7 @@ PROGRAM mkEASETilesParam !!! close (11, status = 'keep') !!! !!! END SUBROUTINE write_tilfile - END PROGRAM + + +END PROGRAM mkEASETilesParam From 212bbbed76c7b9ef7f5a0e876a2b9d55f1ff70e6 Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Thu, 22 Jun 2023 17:57:01 -0400 Subject: [PATCH 14/16] fixed old error in comments (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 2fcff34f0..5893db488 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -206,8 +206,8 @@ PROGRAM mkEASETilesParam ! Coarsen 10-arcsec mask to 30-arcsec "raster" grid - nc = 43200 ! Number of rows of raster grid - nr = 21600 ! Number of columns of raster grid + nc = 43200 ! Number of columns of raster grid + nr = 21600 ! Number of rows of raster grid regrid = .true. From 7e065c144eedd5e849d99669fb00cd5199a5d17c Mon Sep 17 00:00:00 2001 From: Rolf Reichle <54944691+gmao-rreichle@users.noreply.github.com> Date: Tue, 27 Jun 2023 09:28:17 -0400 Subject: [PATCH 15/16] remove unnecessary dependency on rmTinyCatchParaMod (mkEASETilesParam.F90) --- .../Utils/Raster/makebcs/mkEASETilesParam.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 5893db488..d64bf0967 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -33,7 +33,6 @@ PROGRAM mkEASETilesParam use EASE_conv, only : EASE_extent, EASE_convert, EASE_inverse use rmTinyCatchParaMod, only : i_raster, j_raster, SRTM_maxcat use rmTinyCatchParaMod, only : RegridRasterReal - use rmTinyCatchParaMod, only : MAKE_BCS_INPUT_DIR use process_hres_data, only : histogram use MAPL_SortMod use MAPL_ConstantsMod @@ -109,7 +108,8 @@ PROGRAM mkEASETilesParam character(len=128) :: usage1, usage2 character(len=128) :: Iam = "mkEASETilesParam" character(len=512) :: fname_mask - + character(len=400) :: MAKE_BCS_INPUT_DIR + ! -------------------------------------------------------------------------------------- call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR ) From 3fffc8f9ea56b81197bb4a505c5503bb369a3fa7 Mon Sep 17 00:00:00 2001 From: Rolf Reichle Date: Tue, 27 Jun 2023 12:24:36 -0400 Subject: [PATCH 16/16] added documentation, cleanup, fixed indentation (mkLandRaster.F90) --- .../Utils/Raster/makebcs/mkLandRaster.F90 | 934 +++++++++--------- 1 file changed, 481 insertions(+), 453 deletions(-) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 index 5c2494f28..1ecb9a0c3 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkLandRaster.F90 @@ -1,427 +1,468 @@ #define I_AM_MAIN #include "MAPL_ErrLog.h" - Program MakeLandRaster - use MAPL_ExceptionHandling - use LogRectRasterizeMod - use MAPL_HashMod - use process_hres_data - use MAPL_SortMod - use rmTinyCatchParaMod, ONLY: SRTM_maxcat, MAKE_BCS_INPUT_DIR - use MAPL_Constants, only: PI=>MAPL_PI_R8 -! Program to create a surface raster file at 2.5' that has -! the ocean divided with a regular lat-lon DE grid. Its inputs -! are Sarith's formatted 2.5' raster of the Pfafstetter catchments with -! ocean points set to 0. Output is an unformatted integer*4 2.5" raster -! with the ocean points set to negative numbers. - - implicit none - - - integer, parameter :: InUNIT = 20 - integer, parameter :: nx0 = 8640 - integer, parameter :: ny0 = 4320 - - integer :: command_argument_count - integer :: i,j,k,n, status,ncid, ncid2 - integer :: ip, nxt - integer :: type, maxtiles, nx, ny - integer :: count0,count1,count_rate - - real(kind=8) :: dx, dy, d2r ! Grid spacing of raster grid - real(kind=8) :: xmin, ymin, xmax, ymax, xs, ys, da - - real(kind=8), allocatable :: cc(:), ss(:) - real(kind=8) , allocatable :: rTable(:,:) - - integer, pointer :: Raster(:,:) - integer, allocatable, target :: Raster0(:,:) - integer, allocatable :: iTable(:,:) - - integer :: Pix, IceID, LakeID - integer :: CreateHash, IncrementHash, Hash - - logical :: DoZip - logical :: Verb - logical :: regrid, reynolds_sst = .false. - - real(kind=8) :: VV(4) - - ! ESA/SRTM ocean/land/ice/lake mask parameters - ! -------------------------------------------- - - integer, parameter :: nc_esa = 129600, nr_esa = 64800 - integer, allocatable, target, dimension (:,:) :: geos_msk, geos_msk2 - REAL, allocatable, DIMENSION (:) :: loc_val - INTEGER, ALLOCATABLE, DIMENSION (:) :: density, loc_int - REAL, allocatable, DIMENSION (:) :: loc_val2 - INTEGER, ALLOCATABLE, DIMENSION (:) :: density2, loc_int2 - logical, dimension (:), allocatable :: unq_mask, unq_mask2 - integer, pointer , dimension (:,:) :: subset, subset2 - integer :: dx_esa, dy_esa, NBINS, NPLUS,NPLUS2, NBINS2, NonReynold - - integer*8, allocatable, dimension (:) :: SRTM_catid - - character*1 :: Opt - character*128 :: arg - character*128 :: TilFile, RstFile - character*128 :: Tildir, Rstdir - character*128 :: GridName - character*128 :: InputFile - character*128 :: MaskFile - character*128 :: & - Usage = "mkLandRaster -x nx -y ny -v -h -z -t maxtiles -l LandFile -g GridName" - character*128 :: Iam = "MakeLandRaster" - include 'netcdf.inc' - call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) - -! Process Arguments -!------------------ - - nx = nx0 - ny = ny0 - GridName = "Pfafstetter" - DoZip = .false. ! No zipping - Verb = .false. ! Run quiet - tildir = 'til/' ! Write in current dir - rstdir = 'rst/' ! Write in current dir - maxtiles = 50000 - InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/global.cat_id.catch.DL' +Program MakeLandRaster + + use MAPL_ExceptionHandling + use LogRectRasterizeMod + use MAPL_HashMod + use process_hres_data + use MAPL_SortMod + use rmTinyCatchParaMod, ONLY: SRTM_maxcat, RegridRaster + use MAPL_Constants, only: PI=>MAPL_PI_R8 + + ! Program to create a surface raster file that has + ! the ocean divided with a regular lat-lon DE grid. Its inputs + ! are Sarith's formatted raster of the Pfafstetter catchments with + ! ocean points set to 0. Output is an unformatted integer*4 raster + ! with the ocean points set to negative numbers. + ! + ! Usage = "mkLandRaster -x nx -y ny -v -h -z -t maxtiles -g GridName" + ! + ! -x nx : # raster grid cells in longitude direction + ! -y ny : # raster grid cells in latitude direction + ! -v : verbose + ! -t maxtiles : max # tiles expected + ! + ! [The following command line arguments are not used by makebcs as of Jun 2023.] + ! + ! -h : "here" [??] (write *.til and *.rst files to current directory) + ! -z : zip + ! -g GridName : "Pfafstetter" [??] + ! + ! Since Icarus-NLv3 (ca. 2020), typically used with: + ! - GEOS5_10arcsec_mask* file + ! - 30-arcsec raster (nx=43200, ny=21600) + ! + ! Requires the following environment variables: + ! - MAKE_BCS_INPUT_DIR + ! - MASKFILE + ! + ! ---------------------------------------------------------------------------- + + implicit none + + integer, parameter :: InUNIT = 20 + integer, parameter :: nx0 = 8640 ! default long resolution (typically overwritten by command line) + integer, parameter :: ny0 = 4320 ! default lat resolution (typically overwritten by command line) + + integer :: command_argument_count + integer :: i,j,k,n, status,ncid, ncid2 + integer :: ip, nxt + integer :: type, maxtiles, nx, ny + integer :: count0,count1,count_rate + + real(kind=8) :: dx, dy, d2r ! Grid spacing of raster grid + real(kind=8) :: xmin, ymin, xmax, ymax, xs, ys, da + + real(kind=8), allocatable :: cc(:), ss(:) + + real(kind=8) , allocatable :: rTable( :,:) + + integer, pointer :: Raster( :,:) + integer, allocatable, target :: Raster0(:,:) + integer, allocatable :: iTable( :,:) + + integer :: Pix, IceID, LakeID + integer :: CreateHash, IncrementHash, Hash + + logical :: DoZip + logical :: Verb + logical :: regrid=.false., reynolds_sst=.false. + + real(kind=8) :: VV(4) + + ! ESA/SRTM ocean/land/ice/lake mask parameters + ! -------------------------------------------- + + integer, parameter :: nc_esa = 129600, nr_esa = 64800 ! 10-arcsec resolution of GEOS5_10arcsec_mask* + + integer, allocatable, dimension(:,:), target :: geos_msk, geos_msk2 + REAL, allocatable, DIMENSION(:) :: loc_val + INTEGER, ALLOCATABLE, DIMENSION(:) :: density, loc_int + REAL, allocatable, DIMENSION(:) :: loc_val2 + INTEGER, ALLOCATABLE, DIMENSION(:) :: density2, loc_int2 + logical, allocatable, dimension(:) :: unq_mask, unq_mask2 + integer, dimension(:,:), pointer :: subset, subset2 + + integer :: dx_esa, dy_esa, NBINS, NPLUS,NPLUS2, NBINS2, NonReynold + + integer*8, allocatable, dimension(:) :: SRTM_catid + + character*1 :: Opt + character*128 :: arg + character*128 :: TilFile, RstFile + character*128 :: Tildir, Rstdir + character*128 :: GridName + character*512 :: InputFile + character*128 :: MaskFile + character*400 :: MAKE_BCS_INPUT_DIR + + character*128 :: & + Usage = "mkLandRaster -x nx -y ny -v -h -z -t maxtiles -g GridName" + character*128 :: Iam = "MakeLandRaster" + + ! ------------------------------------------------------------------------ + + include 'netcdf.inc' + + call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR ) + + ! Process command line arguments + + ! set defaults (typically overwritten by command line) + + nx = nx0 + ny = ny0 + GridName = "Pfafstetter" + DoZip = .false. ! No zipping + Verb = .false. ! Run quiet + tildir = 'til/' ! Write in current dir + rstdir = 'rst/' ! Write in current dir + maxtiles = 50000 + + ! read command line arguments + + I = command_argument_count() + + if(I > 13) then + print *, "Wrong Number of arguments: ", i + print *, trim(Usage) + call exit(1) + end if + + nxt = 1 + + do while(nxt<=I) + call get_command_argument(nxt,arg) + if(arg(1:1)/='-') then + print *, trim(Usage) + call exit(1) + end if + opt=arg(2:2) + if(len(trim(arg))==2) then + if(scan(opt,'zvh')==0) then + nxt = nxt + 1 + call get_command_argument(nxt,arg) + end if + else + arg = arg(3:) + end if + select case (opt) + case ('x') + read(arg,'(i6)') nx + case ('y') + read(arg,'(i6)') ny + case ('g') + Gridname = arg + case ('v') + Verb = .true. + case ('z') + DoZip = .true. + case ('h') + tildir = '' + rstdir = '' + case ('t') + read(arg,*) maxtiles + case default + print *, trim(Usage) + call exit(1) + end select - I = command_argument_count() - - if(I > 13) then - print *, "Wrong Number of arguments: ", i - print *, trim(Usage) - call exit(1) - end if - - nxt = 1 - - do while(nxt<=I) - call get_command_argument(nxt,arg) - if(arg(1:1)/='-') then - print *, trim(Usage) - call exit(1) - end if - opt=arg(2:2) - if(len(trim(arg))==2) then - if(scan(opt,'zvh')==0) then - nxt = nxt + 1 - call get_command_argument(nxt,arg) - end if - else - arg = arg(3:) - end if - select case (opt) - case ('x') - read(arg,'(i6)') nx - case ('y') - read(arg,'(i6)') ny - case ('l') - InputFile = arg - case ('g') - Gridname = arg - case ('v') - Verb = .true. - case ('z') - DoZip = .true. - case ('h') - tildir = '' - rstdir = '' - case ('t') - read(arg,*) maxtiles - case default - print *, trim(Usage) - call exit(1) - end select - - nxt = nxt + 1 - end do - -! Check for the 10 arc-sec MaskFile (SM) - - call get_environment_variable ("MASKFILE" ,MaskFile ) - - print *,'Using Mask file : ',trim(MaskFile) - - call system_clock(count0,count_rate) - - if(DoZip) then - TilFile = trim(tildir)//trim(Gridname)//'.til.gz' - RstFile = trim(rstdir)//trim(Gridname)//'.rst.gz' - else - TilFile = trim(tildir)//trim(Gridname)//'.til' - RstFile = trim(rstdir)//trim(Gridname)//'.rst' - end if - - allocate(cc(nx), stat=STATUS); VERIFY_(STATUS) - allocate(ss(nx), stat=STATUS); VERIFY_(STATUS) - - allocate(iTable(0:3,maxtiles),stat=status) - VERIFY_(STATUS) - allocate(rTable(1:4,maxtiles),stat=status) - VERIFY_(STATUS) - -! Compute ys and xs at centers of raster cell - - xmin = -180.0_8 - xmax = 180.0_8 - ymin = -90.0_8 - ymax = 90.0_8 - - dx = (xmax-xmin)/nx - dy = (ymax-ymin)/ny - d2r = (2._8*PI)/360.0_8 - - do i = 1,nx - xs = (xmin + dx*(i-0.5_8))*d2r - cc(i) = cos(xs) - ss(i) = sin(xs) - enddo - - InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/'//trim(MaskFile) - - if (index(trim(MaskFile),'GEOS5_10arcsec_mask')/=0) then - ! 10 arcsec new mask - LakeID = 190000000 - IceID = 200000000 - dx_esa = nc_esa / nx - dy_esa = nr_esa / ny - regrid = .false. - regrid = nx/=nc_esa .or. ny/=nr_esa - allocate(SRTM_catid (1:SRTM_maxcat)) - allocate(geos_msk (1:nc_esa,1:dy_esa)) - allocate (raster (1:nx, 1:ny)) - - InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/'//trim(MaskFile) - - status = NF_OPEN (InputFile, NF_NOWRITE, ncid) - - if(status /=0) then - PRINT *, NF_STRERROR(STATUS) - print *, 'Problem with NF_OPEN',InputFile - endif - - status = NF_GET_VARA_INT64 (ncid,3,(/1/),(/SRTM_maxcat/),SRTM_catid ) - - if(index(trim(MaskFile),'GEOS5_10arcsec_mask_freshwater-lakes')/=0) then - - ! Special case for Reynolds SST - i.e. Great lakes and the Caspian - ! sea are treated as freshwater bodies. We make sure both land-tiles don't - ! when the 30 arc-sec mask (raster file) is created from the 10 arc-sec - ! GEOS5_10arcsec_mask_freshwater-lakes.nc - - print *, 'Using Reynolds SSTs MASKFILE',trim(MaskFile) - reynolds_sst = .true. - - InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/GEOS5_10arcsec_mask.nc' - status = NF_OPEN (InputFile, NF_NOWRITE, ncid2) - allocate(geos_msk2 (1:nc_esa,1:dy_esa)) - endif - - if(Verb) then - call system_clock(count1) - print *, 'Opened land file. Time = ', (count1-count0)/float(count_rate) - count0 = count1 - end if - - do j=1,ny - status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' - - if(reynolds_sst) & - status = NF_GET_VARA_INT (ncid2,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk2) ! Read 10-arcsec MERRA2/OSTIA mask file - - if (regrid) then - if(status /=0) then - PRINT *, NF_STRERROR(STATUS) - print *, 'Problem with NF_GET_VARA_INT',InputFile,status - endif - - do i = 1,nx - if (associated (subset)) NULLIFY (subset) - subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains 10-arcsec pixels that lie within the raster grid cell at i,j - if(maxval (subset) > SRTM_maxcat) then - where (subset == LakeID)subset = SRTM_maxcat + 1 - where (subset == IceID) subset = SRTM_maxcat + 2 - endif - - if (maxval (subset) > 0) then ! check whether there are Non-ocean 10-arcsec pixels - - NPLUS = count(subset >= 1 .and. subset <= SRTM_maxcat + 2) ! Count non-ocean 10-arcsec pixels within - allocate (loc_int (1:NPLUS)) - allocate (unq_mask(1:NPLUS)) - loc_int = pack(subset,mask = (subset >= 1 .and. subset <= SRTM_maxcat + 2)) ! loc_int contains catch_indices of non-ocean 10-arcsec pixels - call MAPL_Sort (loc_int) - unq_mask = .true. - do n = 2,NPLUS - unq_mask(n) = .not.(loc_int(n) == loc_int(n-1)) ! count number of unique numbers in loc_int for binning - end do - NBINS = count(unq_mask) - - if(reynolds_sst) then - if (associated (subset2)) NULLIFY (subset2) - subset2 => geos_msk2 ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) - if(maxval (subset2) > SRTM_maxcat) then - where (subset2 == LakeID)subset2 = SRTM_maxcat + 1 - where (subset2 == IceID) subset2 = SRTM_maxcat + 2 - endif - NPLUS2 = count(subset2 >= 1 .and. subset2 <= SRTM_maxcat + 2) ! Count non-ocean 10-arcsec pixels within - allocate (loc_int2 (1:NPLUS2)) - allocate (unq_mask2(1:NPLUS2)) - loc_int2 = pack(subset2,mask = (subset2 >= 1 .and. subset2 <= SRTM_maxcat + 2)) ! loc_int contains catch_indices of non-ocean 10-arcsec pixels - call MAPL_Sort (loc_int2) - unq_mask2 = .true. - do n = 2,NPLUS2 - unq_mask2(n) = .not.(loc_int2(n) == loc_int2(n-1)) ! count number of unique numbers in loc_int for binning - end do - NBINS2 = count(unq_mask2) - endif - - if (NBINS > 1) then - allocate(loc_val (1:NBINS)) - allocate(density (1:NBINS)) - loc_val = 1.*pack(loc_int,mask =unq_mask) ! loc_val contains available non-ocean catch_indices within the i,j grid cell, - ! Those numbers will be used as bin values - call histogram (dx_esa*dy_esa, NBINS, density, loc_val, real(subset)) ! density is the pixel count for each bin value - raster (i,j) = loc_val (maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j - - if(reynolds_sst) then - allocate(loc_val2 (1:NBINS2)) - allocate(density2 (1:NBINS2)) - loc_val2 = 1.*pack(loc_int2,mask =unq_mask2) - call histogram (dx_esa*dy_esa, NBINS2, density2, loc_val2, real(subset2)) - NonReynold = loc_val2 (maxloc(density2,1)) - if(NonReynold /= raster (i,j)) then - if(NonReynold >= 1 .and. NonReynold <= SRTM_maxcat) raster (i,j) = NonReynold ! reset to NonReynold land - endif - endif - - deallocate (loc_val, density) - if(reynolds_sst) deallocate (loc_val2, density2) - else - raster (i,j) = loc_int (1) - endif - deallocate (loc_int, unq_mask) - if(reynolds_sst) deallocate (loc_int2, unq_mask2) - - if(raster (i,j) == SRTM_maxcat + 1) raster (i,j) = LakeID - if(raster (i,j) == SRTM_maxcat + 2) raster (i,j) = IceID - endif - end do - else - raster (:,j) = geos_msk (:,1) - endif - end do - - status = NF_CLOSE (ncid) - if(reynolds_sst) status = NF_CLOSE (ncid2) - else - LakeID = 6190000 - IceID = 6200000 - ! Input file: (formatted for now) - - allocate(raster0(nx0,ny0),stat=STATUS); VERIFY_(STATUS) - - regrid = nx/=nx0 .or. ny/=ny0 + nxt = nxt + 1 + end do + + ! --------------------------------------------------------------------------------- + ! + ! remap mask file onto raster grid + + call get_environment_variable( "MASKFILE", MaskFile ) + + print *, 'Using Mask file : ', trim(MaskFile) + + call system_clock(count0,count_rate) + + InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/'//trim(MaskFile) + + if (index(trim(MaskFile),'GEOS5_10arcsec_mask')/=0) then + + ! 10 arcsec new mask (nc_esa-by-nr_esa) + + LakeID = 190000000 + IceID = 200000000 - if(regrid) then - allocate(raster(nx,ny),stat=STATUS); VERIFY_(STATUS) - else - raster => raster0 - end if + dx_esa = nc_esa / nx ! # of ESA columns within a raster grid cell in x-dim (=3 for 30-arcsec raster; nx=nc=43200) + dy_esa = nr_esa / ny ! # of ESA rows within a raster grid cell in y-dim (=3 for 30-arcsec raster; ny=nr=21600) - open (InUnit,file=InputFile, form='formatted',status='old') + regrid = nx/=nc_esa .or. ny/=nr_esa - if(Verb) then - call system_clock(count1) - print *, 'Opened land file. Time = ', (count1-count0)/float(count_rate) - count0 = count1 - end if + allocate(SRTM_catid(1:SRTM_maxcat )) + allocate(geos_msk (1:nc_esa, 1:dy_esa)) + allocate(raster (1:nx, 1:ny )) - do j=1,ny0 - read (InUnit,*) raster0(:,j) - enddo + InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/'//trim(MaskFile) -! call ReadRaster(InputFile,raster) ! reead unformatted + status = NF_OPEN (InputFile, NF_NOWRITE, ncid) - if(regrid) then - call RegridRaster(Raster0,raster) - endif - - endif + if(status /=0) then + PRINT *, NF_STRERROR(STATUS) + print *, 'Problem with NF_OPEN',InputFile + endif + + status = NF_GET_VARA_INT64 (ncid,3,(/1/),(/SRTM_maxcat/),SRTM_catid ) + + if(index(trim(MaskFile),'GEOS5_10arcsec_mask_freshwater-lakes')/=0) then - if(Verb) then - call system_clock(count1) - print *, 'Read land file. Time = ', (count1-count0)/float(count_rate) - count0 = count1 - end if + ! Special case for Reynolds SST - i.e. Great lakes and the Caspian + ! sea are treated as freshwater bodies. We make sure both land-tiles don't + ! when the 30 arc-sec mask (raster file) is created from the 10 arc-sec + ! GEOS5_10arcsec_mask_freshwater-lakes.nc + + print *, 'Using Reynolds SSTs MASKFILE',trim(MaskFile) + reynolds_sst = .true. + + InputFile = trim(MAKE_BCS_INPUT_DIR)//'/shared/mask/GEOS5_10arcsec_mask.nc' + status = NF_OPEN (InputFile, NF_NOWRITE, ncid2) + allocate(geos_msk2 (1:nc_esa,1:dy_esa)) + endif + + if(Verb) then + call system_clock(count1) + print *, 'Opened mask file. Time = ', (count1-count0)/float(count_rate) + count0 = count1 + end if + + do j=1,ny + status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' + + if(reynolds_sst) & + status = NF_GET_VARA_INT (ncid2,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk2) ! Read 10-arcsec MERRA2/OSTIA mask file + + if (regrid) then + if(status /=0) then + PRINT *, NF_STRERROR(STATUS) + print *, 'Problem with NF_GET_VARA_INT',InputFile,status + endif + + do i = 1,nx + if (associated (subset)) NULLIFY (subset) + subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains 10-arcsec pixels that lie within the raster grid cell at i,j + if(maxval (subset) > SRTM_maxcat) then + where (subset == LakeID)subset = SRTM_maxcat + 1 + where (subset == IceID) subset = SRTM_maxcat + 2 + endif + + if (maxval (subset) > 0) then ! check whether there are Non-ocean 10-arcsec pixels + + NPLUS = count(subset >= 1 .and. subset <= SRTM_maxcat + 2) ! Count non-ocean 10-arcsec pixels within + allocate (loc_int (1:NPLUS)) + allocate (unq_mask(1:NPLUS)) + loc_int = pack(subset,mask = (subset >= 1 .and. subset <= SRTM_maxcat + 2)) ! loc_int contains catch_indices of non-ocean 10-arcsec pixels + call MAPL_Sort (loc_int) + unq_mask = .true. + do n = 2,NPLUS + unq_mask(n) = .not.(loc_int(n) == loc_int(n-1)) ! count number of unique numbers in loc_int for binning + end do + NBINS = count(unq_mask) + + if(reynolds_sst) then + if (associated (subset2)) NULLIFY (subset2) + subset2 => geos_msk2 ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) + if(maxval (subset2) > SRTM_maxcat) then + where (subset2 == LakeID)subset2 = SRTM_maxcat + 1 + where (subset2 == IceID) subset2 = SRTM_maxcat + 2 + endif + NPLUS2 = count(subset2 >= 1 .and. subset2 <= SRTM_maxcat + 2) ! Count non-ocean 10-arcsec pixels within + allocate (loc_int2 (1:NPLUS2)) + allocate (unq_mask2(1:NPLUS2)) + loc_int2 = pack(subset2,mask = (subset2 >= 1 .and. subset2 <= SRTM_maxcat + 2)) ! loc_int contains catch_indices of non-ocean 10-arcsec pixels + call MAPL_Sort (loc_int2) + unq_mask2 = .true. + do n = 2,NPLUS2 + unq_mask2(n) = .not.(loc_int2(n) == loc_int2(n-1)) ! count number of unique numbers in loc_int for binning + end do + NBINS2 = count(unq_mask2) + endif + + if (NBINS > 1) then + allocate(loc_val (1:NBINS)) + allocate(density (1:NBINS)) + loc_val = 1.*pack(loc_int,mask =unq_mask) ! loc_val contains available non-ocean catch_indices within the i,j grid cell, + ! Those numbers will be used as bin values + call histogram (dx_esa*dy_esa, NBINS, density, loc_val, real(subset)) ! density is the pixel count for each bin value + raster (i,j) = loc_val (maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j + + if(reynolds_sst) then + allocate(loc_val2 (1:NBINS2)) + allocate(density2 (1:NBINS2)) + loc_val2 = 1.*pack(loc_int2,mask =unq_mask2) + call histogram (dx_esa*dy_esa, NBINS2, density2, loc_val2, real(subset2)) + NonReynold = loc_val2 (maxloc(density2,1)) + if(NonReynold /= raster (i,j)) then + if(NonReynold >= 1 .and. NonReynold <= SRTM_maxcat) raster (i,j) = NonReynold ! reset to NonReynold land + endif + endif + + deallocate (loc_val, density) + if(reynolds_sst) deallocate (loc_val2, density2) + else + raster (i,j) = loc_int (1) + endif + deallocate (loc_int, unq_mask) + if(reynolds_sst) deallocate (loc_int2, unq_mask2) + + if(raster (i,j) == SRTM_maxcat + 1) raster (i,j) = LakeID + if(raster (i,j) == SRTM_maxcat + 2) raster (i,j) = IceID + endif + end do + else + raster (:,j) = geos_msk (:,1) + endif + end do + + status = NF_CLOSE (ncid) + if(reynolds_sst) status = NF_CLOSE (ncid2) + + else + + ! old mask (before Icarus-NLv3): 'global.cat_id.catch.DL' + + LakeID = 6190000 + IceID = 6200000 + + ! Input file: (formatted for now) + + allocate(raster0(nx0,ny0),stat=STATUS); VERIFY_(STATUS) + + regrid = nx/=nx0 .or. ny/=ny0 + + if(regrid) then + allocate(raster(nx,ny),stat=STATUS); VERIFY_(STATUS) + else + raster => raster0 + end if -! Loop over ocean latitudes + open (InUnit,file=InputFile, form='formatted',status='old') if(Verb) then - write (6, '(A)', advance='NO') 'Processing catchment data' + call system_clock(count1) + print *, 'Opened land file. Time = ', (count1-count0)/float(count_rate) + count0 = count1 end if - ip = 0 - Hash = CreateHash(4*1024) - - LATITUDES: do j=1,ny - - where(raster(:,j)>=LakeID) - raster(:,j) = (raster(:,j)/10000)*10000 - end where - - ys = ymin + dy*(j-0.5_8) - da = (sin(d2r*(ys+0.5*dy)) - & - sin(d2r*(ys-0.5*dy)) )*(dx*d2r) - - vv(3) = da*ys - vv(4) = da - - LONGITUDES: do i=1,nx - - vv(1) = ss(i)*da - vv(2) = cc(i)*da - - Pix = raster(i,j) - - k = MAPL_HashIncrement(Hash,Pix); - - raster(i,j) = k - - if(k<=ip) then ! Bump the counter and the lon, lat, and area sums - iTable(1,k) = iTable(1,k) + 1 - rTable(1,k) = rTable(1,k) + vv(1) - rTable(2,k) = rTable(2,k) + vv(2) - rTable(3,k) = rTable(3,k) + vv(3) - rTable(4,k) = rTable(4,k) + vv(4) - else ! We have a new tile in the exchange grid - ip = ip + 1 - if(ip>maxtiles) then - print *, "Exceeded maxtiles = ", maxtiles, & - ". Use -t option to increase." - call exit(1) - end if - - select case (Pix) - case (0) - type = 0 - case(190000000,6190000) - type = 19 - case(200000000,6200000) - type = 20 - case default - type = 100 - end select - - iTable(0 ,ip) = type - iTable(1 ,ip) = 1 - rTable(:4,ip) = vv - - iTable( 2,ip) = Pix - iTable( 3,ip) = 1 - end if + do j=1,ny0 + read (InUnit,*) raster0(:,j) + enddo + + ! call ReadRaster(InputFile,raster) ! reead unformatted + + if(regrid) then + call RegridRaster(Raster0,raster) + endif + + endif ! MaskFile==GEOS5_10arcsec_mask + + if(Verb) then + call system_clock(count1) + print *, 'Done reading land file. Time = ', (count1-count0)/float(count_rate) + count0 = count1 + end if + + ! ------------------------------------------------------------------------------- + ! + ! Loop through raster grid and assign tile IDs [??] + + if(Verb) then + write (6, '(A)', advance='NO') 'Processing catchment data' + end if + + ip = 0 + Hash = CreateHash(4*1024) + + allocate(cc(nx), stat=STATUS); VERIFY_(STATUS) + allocate(ss(nx), stat=STATUS); VERIFY_(STATUS) + + allocate(iTable(0:3,maxtiles),stat=status) + VERIFY_(STATUS) + allocate(rTable(1:4,maxtiles),stat=status) + VERIFY_(STATUS) + + ! lat/lon coords of raster grid cells + + xmin = -180.0_8 ! western edge of grid + xmax = 180.0_8 ! eastern edge of grid + ymin = -90.0_8 ! southern edge of grid + ymax = 90.0_8 ! northern edge of grid + + dx = (xmax-xmin)/nx ! raster grid spacing (longitude) + dy = (ymax-ymin)/ny ! raster grid spacing (latitude) + + d2r = (2._8*PI)/360.0_8 ! degree-to-radians conversion + + ! Compute sin and cos at longitude centers of raster cells + + do i = 1,nx + xs = (xmin + dx*(i-0.5_8))*d2r ! radians + cc(i) = cos(xs) + ss(i) = sin(xs) + enddo + + LATITUDES: do j=1,ny + + where(raster(:,j)>=LakeID) + raster(:,j) = (raster(:,j)/10000)*10000 + end where + + ys = ymin + dy*(j-0.5_8) + da = (sin(d2r*(ys+0.5*dy)) - & + sin(d2r*(ys-0.5*dy)) )*(dx*d2r) + + vv(3) = da*ys + vv(4) = da + + LONGITUDES: do i=1,nx + + vv(1) = ss(i)*da + vv(2) = cc(i)*da + + Pix = raster(i,j) + + k = MAPL_HashIncrement(Hash,Pix); + + raster(i,j) = k + + if(k<=ip) then ! Bump the counter and the lon, lat, and area sums + iTable(1,k) = iTable(1,k) + 1 + rTable(1,k) = rTable(1,k) + vv(1) + rTable(2,k) = rTable(2,k) + vv(2) + rTable(3,k) = rTable(3,k) + vv(3) + rTable(4,k) = rTable(4,k) + vv(4) + else ! We have a new tile in the exchange grid + ip = ip + 1 + if(ip>maxtiles) then + print *, "Exceeded maxtiles = ", maxtiles, & + ". Use -t option to increase." + call exit(1) + end if + + select case (Pix) + case (0) + type = 0 ! ocean + case(190000000,6190000) + type = 19 ! lake + case(200000000,6200000) + type = 20 ! landice + case default + type = 100 ! land + end select + + iTable(0 ,ip) = type + iTable(1 ,ip) = 1 + rTable(:4,ip) = vv + + iTable( 2,ip) = Pix + iTable( 3,ip) = 1 + end if end do LONGITUDES @@ -430,9 +471,9 @@ Program MakeLandRaster write (6, '(A)', advance='NO') '.' end if end if - + end do LATITUDES - + close(InUnit) if(Verb) then @@ -444,12 +485,12 @@ Program MakeLandRaster type = sum(iTable(1,:ip)) print *, type, nx*ny, nx*ny-type end if - + call DestroyHash(Hash) - -! Compute proper longitude and latitude in degrees and compress -! the real table for WriteTiling. - + + ! Compute proper longitude and latitude in degrees and compress + ! the real table for WriteTiling. + if(Verb) print *, "Computing weighted lons and lats..." do k=1,ip @@ -457,63 +498,50 @@ Program MakeLandRaster rTable(2,k) = rTable(3,k)/rTable(4,k) rTable(3,k) = rTable(4,k) end do - -! Sort table by the first grid type in descending order. + + ! Sort table by the first grid type in descending order. if(Verb) print *, "Sorting..." call SortTiling(Raster,rTable(:,:ip),iTable(:,:ip)) - + if(Verb) then call system_clock(count1) - print *, "Done Sorting. Time = ", (count1-count0)/float(count_rate) + print *, "Done sorting. Time = ", (count1-count0)/float(count_rate) count0 = count1 end if - + do k=1,ip iTable( 3,k) = 1 end do + + ! Write *.rst and *.til files + + if(DoZip) then + TilFile = trim(tildir)//trim(Gridname)//'.til.gz' + RstFile = trim(rstdir)//trim(Gridname)//'.rst.gz' + else + TilFile = trim(tildir)//trim(Gridname)//'.til' + RstFile = trim(rstdir)//trim(Gridname)//'.rst' + end if if(Verb) print *, 'Writing til file...' call WriteTiling(TilFile, (/Gridname/), (/ip/), (/1/), (/ip/), & - nx, ny, iTable(:,:ip), rTable(:4,:ip), Dozip, Verb ) - + nx, ny, iTable(:,:ip), rTable(:4,:ip), Dozip, Verb ) + if(Verb) print *, 'Writing raster file...' call WriteRaster( RstFile, Raster, DoZip) if(Verb) then call system_clock(count1) - print *, "Done Writing. Time = ", (count1-count0)/float(count_rate) + print *, "Done writing. Time = ", (count1-count0)/float(count_rate) count0 = count1 end if - -! All done - - if(Verb) print * , 'Terminated Normally' + + ! All done + + if(Verb) print * , 'Terminated normally' call exit(0) - -contains - -subroutine RegridRaster(Rin,Rout) - - integer, intent(IN) :: Rin(:,:) - integer, intent(OUT) :: Rout(:,:) - - real(kind=8) :: xx, yy - integer :: i,j,ii,jj - - xx = size(Rin ,1)/float(size(Rout,1)) - yy = size(Rin ,2)/float(size(Rout,2)) - - do j=1,size(Rout,2) - jj = (j-1)*yy + 1 - do i=1,size(Rout,1) - ii = (i-1)*xx + 1 - Rout(i,j) = Rin(ii,jj) - end do - end do - -end subroutine RegridRaster - - + end Program MakeLandRaster - + +! =========================== EOF ========================================================