diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 index 82879861c..9619c906d 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 @@ -3,25 +3,26 @@ PROGRAM mkCatchParam -! !INTERFACE: -! -! !ARGUMENTS: -! -! Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV " -! -x: Size of longitude dimension of input raster. DEFAULT: 8640 -! -y: Size of latitude dimension of input raster. DEFAULT: 4320 -! -b: Position of dateline w.r.t. first grid cell boundaries. DEFAULT: DC (dateline-on-center) -! -g: Gridname (name of the .til or .rst file without file extension) -! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09, v11, v12 and v13 ) -! -! -! This program is good to generate -! model, vegetation, soil, and MODIS albedo parameter files for the -! catchment model implementation -! -! Sarith Mahanama - March 23, 2012 -! Email: sarith.p.mahanama@nasa.gov - use MAPL, only: MAPL_ease_extent, MAPL_ReadTilingNC4 + ! !INTERFACE: + ! + ! !ARGUMENTS: + ! + ! Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV " + ! -x: Size of longitude dimension of input raster. DEFAULT: 8640 + ! -y: Size of latitude dimension of input raster. DEFAULT: 4320 + ! -b: Position of dateline w.r.t. first grid cell boundaries. DEFAULT: DC (dateline-on-center) + ! -g: Gridname (name of the .til or .rst file without file extension) + ! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09, v11, ...) + ! + ! + ! This program is good to generate + ! model, vegetation, soil, and MODIS albedo parameter files for the + ! catchment model implementation + ! + ! Sarith Mahanama - March 23, 2012 + ! Email: sarith.p.mahanama@nasa.gov + + use MAPL, only: MAPL_ease_extent, MAPL_ReadTilingNC4 use rmTinyCatchParaMod use process_hres_data use MAPL_ExceptionHandling @@ -29,9 +30,9 @@ PROGRAM mkCatchParam ! use module_irrig_params, ONLY : create_irrig_params implicit none - + include 'netcdf.inc' - + ! The default is NC=i_raster=8640, NR=j_raster=4320 via "use rmTinyCatchParaMod", but ! NC and NR are typically overwritten through command-line arguments "-x nx -y ny". ! @@ -75,114 +76,113 @@ PROGRAM mkCatchParam real :: minlon, minlat, maxlon, maxlat, elev integer :: tindex1, pfaf1, n, status -! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------ -! -! NOTE: "!$" is for conditional compilation -! -logical :: running_omp = .false. -! -!$ integer :: omp_get_thread_num, omp_get_num_threads -! -integer :: n_threads=1 - -! ----------- OpenMP PARALLEL ENVIRONMENT ---------------------------- -! -! FIND OUT WHETHER -omp FLAG HAS BEEN SET DURING COMPILATION -! -!$ running_omp = .true. ! conditional compilation -! -! ECHO BASIC OMP VARIABLES -! -!$OMP PARALLEL DEFAULT(NONE) SHARED(running_omp,n_threads) -! -!$OMP SINGLE -! -!$ n_threads = omp_get_num_threads() -! -!$ write (*,*) 'running_omp = ', running_omp -!$ write (*,*) -!$ write (*,*) 'parallel OpenMP with ', n_threads, 'threads' -!$ write (*,*) -!$OMP ENDSINGLE -! -!$OMP CRITICAL -!$ write (*,*) 'thread ', omp_get_thread_num(), ' alive' -!$OMP ENDCRITICAL -! -!$OMP BARRIER -! -!$OMP ENDPARALLEL - - - USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV " - USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " - USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " - USAGE(4) =" -g: Gridname (name of the .til or .rst file *without* file extension) " - USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " - USAGE(6) =" -v: Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08 v09 ) " - USAGE(7) =" -p: if no, it creates catchment_def and nc4 tile files then exits " - -! Process Arguments -!------------------ - - CALL system_clock(count_rate=clock_rate) - CALL get_command (cmd) - inquire(file='clsm/mkCatchParam.log', exist=file_exists) - if(file_exists) then - open (log_file, file ='clsm/mkCatchParam.log', status='old', position='append', form='formatted',action='write') - else - open (log_file, file ='clsm/mkCatchParam.log', status='unknown', form='formatted',action='write') - write (log_file,'(a)')trim(cmd) - write (log_file,'(a)')' ' - endif - withbcs = 'yes' - I = command_argument_count() - if(I < 1 .or. I > 10) then - write (log_file,'(a)') "Wrong Number of arguments: ", i - do j = 1,size(usage) - print "(sp,a100)", Usage(j) - end do - call exit(1) - end if - - nxt = 1 - call get_command_argument(nxt,arg) - do while(arg(1:1)=='-') - opt=arg(2:2) - if(len(trim(arg))==2) then - if(scan(opt,'zh')==0) then - nxt = nxt + 1 - call get_command_argument(nxt,arg) - endif - else - arg = arg(3:) - end if - select case (opt) - case ('x') - read(arg,'(i6)') nc - case ('y') - read(arg,'(i6)') nr - case ('g') - Gridname = trim(arg) - case ('v') - LBCSV = trim(arg) - if (trim(arg).eq."F25") F25Tag = .true. - call init_bcs_config (trim(LBCSV)) ! get bcs details from version string - case ('b') - DL = trim(arg) - case ('p') - withbcs = trim(arg) - case default - do j = 1,size(usage) - print "(sp,a100)", Usage(j) - end do - call exit(1) - end select - nxt = nxt + 1 - call get_command_argument(nxt,arg) - end do + ! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------ + ! + ! NOTE: "!$" is for conditional compilation + ! + logical :: running_omp = .false. + ! + !$ integer :: omp_get_thread_num, omp_get_num_threads + ! + integer :: n_threads=1 - call get_environment_variable ("MASKFILE" ,MaskFile ) + ! ----------- OpenMP PARALLEL ENVIRONMENT ---------------------------- + ! + ! FIND OUT WHETHER -omp FLAG HAS BEEN SET DURING COMPILATION + ! + !$ running_omp = .true. ! conditional compilation + ! + ! ECHO BASIC OMP VARIABLES + ! + !$OMP PARALLEL DEFAULT(NONE) SHARED(running_omp,n_threads) + ! + !$OMP SINGLE + ! + !$ n_threads = omp_get_num_threads() + ! + !$ write (*,*) 'running_omp = ', running_omp + !$ write (*,*) + !$ write (*,*) 'parallel OpenMP with ', n_threads, 'threads' + !$ write (*,*) + !$OMP ENDSINGLE + ! + !$OMP CRITICAL + !$ write (*,*) 'thread ', omp_get_thread_num(), ' alive' + !$OMP ENDCRITICAL + ! + !$OMP BARRIER + ! + !$OMP ENDPARALLEL + + USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV " + USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " + USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " + USAGE(4) =" -g: Gridname (name of the .til or .rst file *without* file extension) " + USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " + USAGE(6) =" -v: Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08 v09 ) " + USAGE(7) =" -p: if no, it creates catchment_def and nc4 tile files then exits " + + ! Process Arguments + !------------------ + + CALL system_clock(count_rate=clock_rate) + CALL get_command (cmd) + inquire(file='clsm/mkCatchParam.log', exist=file_exists) + if(file_exists) then + open (log_file, file ='clsm/mkCatchParam.log', status='old', position='append', form='formatted',action='write') + else + open (log_file, file ='clsm/mkCatchParam.log', status='unknown', form='formatted',action='write') + write (log_file,'(a)')trim(cmd) + write (log_file,'(a)')' ' + endif + withbcs = 'yes' + I = command_argument_count() + if(I < 1 .or. I > 10) then + write (log_file,'(a)') "Wrong Number of arguments: ", i + do j = 1,size(usage) + print "(sp,a100)", Usage(j) + end do + call exit(1) + end if + + nxt = 1 + call get_command_argument(nxt,arg) + do while(arg(1:1)=='-') + opt=arg(2:2) + if(len(trim(arg))==2) then + if(scan(opt,'zh')==0) then + nxt = nxt + 1 + call get_command_argument(nxt,arg) + endif + else + arg = arg(3:) + end if + select case (opt) + case ('x') + read(arg,'(i6)') nc + case ('y') + read(arg,'(i6)') nr + case ('g') + Gridname = trim(arg) + case ('v') + LBCSV = trim(arg) + if (trim(arg).eq."F25") F25Tag = .true. + call init_bcs_config (trim(LBCSV)) ! get bcs details from version string + case ('b') + DL = trim(arg) + case ('p') + withbcs = trim(arg) + case default + do j = 1,size(usage) + print "(sp,a100)", Usage(j) + end do + call exit(1) + end select + nxt = nxt + 1 + call get_command_argument(nxt,arg) + end do + + call get_environment_variable ("MASKFILE" ,MaskFile ) if(trim(Gridname) == '') then write (log_file,'(a)')'Unable to create parameters without til/rst files.... !' @@ -298,540 +298,568 @@ PROGRAM mkCatchParam call exit(0) endif - call MAPL_ReadTilingNC4( trim(fnameTil)//".nc4", iTable = iTable) - - N_land = count(iTable(:,0) == 100) ! n_land = number of land tiles - allocate(tile_j_dum, source = iTable(1:n_land,7)) ! possible used in cti_stats.dat - deallocate (iTable) + tmpstring = 'Step 01: Supplemental tile attributes and nc4-formatted tile file' + fname_tmp = 'clsm/catchment.def' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + if(.not.ease_grid) then + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating catchment def and nc4 tile file...' + call system_clock(clock1) + call supplemental_tile_attributes(nc,nr,regrid,dl,fnameTil, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + else + write (log_file,'(a)')'Skipping step for EASE grid. ' + write (log_file,'(a)')'catchment.def file and tile file should already have been created by mkEASETilesParam.x ' + endif + write (log_file,'(a)')' ' - ! reading from catchment to preserve zero-diff - open (newunit=unit,file='clsm/catchment.def',status='old',action='read',form='formatted', IOSTAT=status) - if (status /=0) then - write (log_file,'(a)')' clsm/cathment.def cannot be opened, exit ' - call exit(1) - endif - read(unit,*) N - if (n /= n_land) then - write (log_file,'(a)')'n_land not consistent between tile file and catchment.def file, exit ' - write (log_file,*) n_land, n - call exit(1) - endif + if (trim(withbcs) == 'no') then + write (log_file,'(a)')'Skipping MOM BCs. BCs will be extracted from the corresponding BCs. ' + close (log_file,status='keep') + call exit(0) + endif - allocate(min_lon(n_land), max_lon(n_land), min_lat(n_land), max_lat(n_land)) - allocate(tile_lat(n_land), tile_lon(n_land)) - allocate(tile_pfs(n_land)) - - do n = 1, N_land - read (unit,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat, elev - min_lon(n) = minlon - max_lon(n) = maxlon - min_lat(n) = minlat - max_lat(n) = maxlat - tile_lon(n)= (minlon + maxlon)/2.0 - tile_lat(n)= (minlat + maxlat)/2.0 - tile_pfs(n)= pfaf1 - end do - close (unit,status='keep') - - inquire(file='clsm/catch_params.nc4', exist=file_exists) - if (.not.file_exists) CALL open_landparam_nc4_files(N_land,process_snow_albedo) - - ! Creating cti_stats.dat - ! ---------------------- + call MAPL_ReadTilingNC4( trim(fnameTil)//".nc4", iTable = iTable) - tmpstring = 'Step 02: Compound Topographic Index (CTI) stats' - fname_tmp = 'clsm/cti_stats.dat' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call cti_stat_file (MaskFile, n_land, tile_pfs, tile_j_dum) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' - - ! Creating vegetation classification files - !----------------------------------------- - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then + N_land = count(iTable(:,0) == 100) ! n_land = number of land tiles + allocate(tile_j_dum, source = iTable(1:n_land,7)) ! possible used in cti_stats.dat + deallocate (iTable) - tmpstring = 'Step 03: Vegetation types using ESA land cover (MOSAIC/Catch)' - fname_tmp = 'clsm/mosaic_veg_typs_fracs' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call ESA2MOSAIC (nc,nr, n_land, tile_pfs, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' + ! reading from catchment to preserve zero-diff + open (newunit=unit,file='clsm/catchment.def',status='old',action='read',form='formatted', IOSTAT=status) + if (status /=0) then + write (log_file,'(a)')' clsm/cathment.def cannot be opened, exit ' + call exit(1) + endif + read(unit,*) N + if (n /= n_land) then + write (log_file,'(a)')'n_land not consistent between tile file and catchment.def file, exit ' + write (log_file,*) n_land, n + call exit(1) + endif - tmpstring = 'Step 04: Vegetation types using ESA land cover (CatchCNCLM40)' - fname_tmp = 'clsm/CLM_veg_typs_fracs' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call ESA2CLM (nc,nr, n_land, tile_lat, tile_pfs, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' + allocate(min_lon(n_land), max_lon(n_land), min_lat(n_land), max_lat(n_land)) + allocate(tile_lat(n_land), tile_lon(n_land)) + allocate(tile_pfs(n_land)) + + do n = 1, N_land + read (unit,*) tindex1,pfaf1,minlon,maxlon,minlat,maxlat, elev + min_lon(n) = minlon + max_lon(n) = maxlon + min_lat(n) = minlat + max_lat(n) = maxlat + tile_lon(n)= (minlon + maxlon)/2.0 + tile_lat(n)= (minlat + maxlat)/2.0 + tile_pfs(n)= pfaf1 + end do + close (unit,status='keep') + + inquire(file='clsm/catch_params.nc4', exist=file_exists) + if (.not.file_exists) CALL open_landparam_nc4_files(N_land,process_snow_albedo) + + ! Creating cti_stats.dat + ! ---------------------- + + tmpstring = 'Step 02: Compound Topographic Index (CTI) stats' + fname_tmp = 'clsm/cti_stats.dat' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call cti_stat_file (MaskFile, n_land, tile_pfs, tile_j_dum) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + ! Creating vegetation classification files + !----------------------------------------- + + if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then + + tmpstring = 'Step 03: Vegetation types using ESA land cover (MOSAIC/Catch)' + fname_tmp = 'clsm/mosaic_veg_typs_fracs' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call ESA2MOSAIC (nc,nr, n_land, tile_pfs, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + tmpstring = 'Step 04: Vegetation types using ESA land cover (CatchCNCLM40)' + fname_tmp = 'clsm/CLM_veg_typs_fracs' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call ESA2CLM (nc,nr, n_land, tile_lat, tile_pfs, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' - else + else - tmpstring = 'Step 03: Vegetation types using IGBP SiB2 land cover (MOSAIC/Catch)' - fname_tmp = 'clsm/mosaic_veg_typs_fracs' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call compute_mosaic_veg_types (nc, nr, regrid, n_land, tile_pfs, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' - - ! Per make_bcs, it looks like there are four possible mask files: - ! - ! GEOS5_10arcsec_mask.nc - ! global.cat_id.catch.DL - ! global.cat_id.catch.GreatLakesCaspian_Updated.DL - ! GEOS5_10arcsec_mask_freshwater-lakes.nc - ! - ! If we are in this else block, we must be using one of the latter three masks. - ! It looks like these latter masks only work for Catchment and not CatchCNCLM[xx] - ! - ! - reichle, 11 Jan 2022 - - write (log_file,'(a)')'NOTE: The selected mask works only for the Catchment model.' - write (log_file,'(a)')' Vegetation types *not* created for CatchCNCLM[xx].' - write (log_file,'(a)')' SKIPPING Step 04 and Step 05 !!!' - write (log_file,'(a)')' ' - - endif - - ! Processing Vegetation Climatology - ! --------------------------------- - - ! creating mapping arrays if necessary - - tmpstring = 'Step 05: Vegetation climatologies' - write (log_file,'(a,a,a)') trim(tmpstring),' ', trim(LAIBCS) - - if((trim(LAIBCS) == 'MODGEO').or.(trim(LAIBCS) == 'GEOLAND2')) then - fname_tmp = 'clsm/lai.GEOLAND2_10-DayClim' - write (log_file,'(a,a)')' --> ', trim(fname_tmp) - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - !allocate (mapgeoland2 (1:40320,1:20160)) - call system_clock(clock1) - call create_mapping (nc,nr,40320,20160,mapgeoland2, n_land, tile_id ) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done create mapping mapgeoland2. Spent ', seconds, " seconds" - lai_name = 'GEOLAND2_10-DayClim/geoland2_' - - write (log_file,'(a)')' Creating '//lai_name - call system_clock(clock1) - if(trim(LAIBCS) == 'GEOLAND2') then - call hres_lai_no_gswp (40320,20160,mapgeoland2, lai_name, n_land, tile_lon, tile_lat) - else - call hres_lai_no_gswp (40320,20160,mapgeoland2, lai_name, n_land, tile_lon, tile_lat, merge=1) - endif - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - ! if(allocated(mapgeoland2)) deallocate (mapgeoland2) - deallocate (mapgeoland2%map) - deallocate (mapgeoland2%ij_index) - write (log_file,'(a)')' Done.' - else - write (log_file,'(a)')' Using existing file.' - endif - endif - - if ((LAIBCS == 'MODGEO').or.(LAIBCS == 'MODIS').or.(MODALB == 'MODIS2')) then - ! allocate (maparc30 (1:43200,1:21600)) - call system_clock(clock1) - call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done create mapping maparc30. Spent ', seconds, " seconds" - endif - - fname_tmp = 'clsm/green.dat' - write (log_file,'(a,a)')' --> ', trim(fname_tmp) - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' - - call system_clock(clock1) - if (trim(LAIBCS) == 'GSWP2') then - call process_gswp2_veg (nc,nr,regrid,'grnFrac',n_land, tile_id) - else - if (size(maparc30%ij_index,1) /= 43200) then - ! allocate (maparc30 (1:43200,1:21600)) - call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) - endif - call hres_gswp2 (43200,21600, maparc30, 'green', n_land, tile_lon, tile_lat) - endif - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif + tmpstring = 'Step 03: Vegetation types using IGBP SiB2 land cover (MOSAIC/Catch)' + fname_tmp = 'clsm/mosaic_veg_typs_fracs' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call compute_mosaic_veg_types (nc, nr, regrid, n_land, tile_pfs, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + ! Per make_bcs, it looks like there are four possible mask files: + ! + ! GEOS5_10arcsec_mask.nc + ! global.cat_id.catch.DL + ! global.cat_id.catch.GreatLakesCaspian_Updated.DL + ! GEOS5_10arcsec_mask_freshwater-lakes.nc + ! + ! If we are in this else block, we must be using one of the latter three masks. + ! It looks like these latter masks only work for Catchment and not CatchCNCLM[xx] + ! + ! - reichle, 11 Jan 2022 + + write (log_file,'(a)')'NOTE: The selected mask works only for the Catchment model.' + write (log_file,'(a)')' Vegetation types *not* created for CatchCNCLM[xx].' + write (log_file,'(a)')' SKIPPING Step 04 and Step 05 !!!' + write (log_file,'(a)')' ' - fname_tmp = 'clsm/lai.dat' - write (log_file,'(a,a)')' --> ', trim(fname_tmp) - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' - redo_modis = .true. - - call system_clock(clock1) - if (trim(LAIBCS) == 'GSWP2') call process_gswp2_veg (nc,nr,regrid,'LAI', n_land, tile_id) - if (trim(LAIBCS) == 'GSWPH') then - if (size(maparc30%ij_index,1) /= 43200) then - ! allocate (maparc30 (1:43200,1:21600)) - call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) - endif - inquire(file='clsm/lai.MODIS_8-DayClim', exist=file_exists) - if (.not.file_exists) call hres_gswp2 (43200,21600, maparc30, 'lai', n_land, tile_lon, tile_lat) - endif - - if (trim(LAIBCS) == 'MODIS') then - lai_name = 'MODIS_8-DayClim/MODIS_' - call hres_lai_no_gswp (43200,21600,maparc30,lai_name, n_land, tile_lon, tile_lon) - endif - - if (trim(LAIBCS) == 'MODGEO') then - lai_name = 'MODIS_8-DayClim/MODIS_' - inquire(file='clsm/lai.MODIS_8-DayClim', exist=file_exists) - if (.not.file_exists)call hres_lai_no_gswp (43200,21600,maparc30,lai_name, n_land, tile_lon, tile_lat, merge=1) - call merge_lai_data (MaskFile, n_land, tile_pfs) - endif - - if (trim(LAIBCS) == 'MODISV6') then - lai_name = 'MCD15A2H.006/MODIS_' - call grid2tile_modis6 (86400,43200,nc,nr,n_land, tile_lon, tile_lat, tile_id, lai_name) - endif + endif + + ! Processing Vegetation Climatology + ! --------------------------------- + + ! creating mapping arrays if necessary + + tmpstring = 'Step 05: Vegetation climatologies' + write (log_file,'(a,a,a)') trim(tmpstring),' ', trim(LAIBCS) + + if((trim(LAIBCS) == 'MODGEO').or.(trim(LAIBCS) == 'GEOLAND2')) then + fname_tmp = 'clsm/lai.GEOLAND2_10-DayClim' + write (log_file,'(a,a)')' --> ', trim(fname_tmp) + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + !allocate (mapgeoland2 (1:40320,1:20160)) + call system_clock(clock1) + call create_mapping (nc,nr,40320,20160,mapgeoland2, n_land, tile_id ) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done create mapping mapgeoland2. Spent ', seconds, " seconds" + lai_name = 'GEOLAND2_10-DayClim/geoland2_' + + write (log_file,'(a)')' Creating '//lai_name + call system_clock(clock1) + if(trim(LAIBCS) == 'GEOLAND2') then + call hres_lai_no_gswp (40320,20160,mapgeoland2, lai_name, n_land, tile_lon, tile_lat) + else + call hres_lai_no_gswp (40320,20160,mapgeoland2, lai_name, n_land, tile_lon, tile_lat, merge=1) + endif + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + ! if(allocated(mapgeoland2)) deallocate (mapgeoland2) + deallocate (mapgeoland2%map) + deallocate (mapgeoland2%ij_index) + write (log_file,'(a)')' Done.' + else + write (log_file,'(a)')' Using existing file.' + endif + endif - if (trim(LAIBCS) == 'GLASSA') then - lai_name = 'GLASS-LAI/AVHRR.v4/GLASS01B02.V04.AYYYY' - call grid2tile_glass (nc,nr, tile_id,lai_name, n_land, tile_lon, tile_lat) - endif + if ((LAIBCS == 'MODGEO').or.(LAIBCS == 'MODIS').or.(MODALB == 'MODIS2')) then + ! allocate (maparc30 (1:43200,1:21600)) + call system_clock(clock1) + call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done create mapping maparc30. Spent ', seconds, " seconds" + endif - if (trim(LAIBCS) == 'GLASSM') then - lai_name = 'GLASS-LAI/MODIS.v4/GLASS01B01.V04.AYYYY' - call grid2tile_glass (nc,nr,tile_id,lai_name, n_land, tile_lon, tile_lat) - endif - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif + fname_tmp = 'clsm/green.dat' + write (log_file,'(a,a)')' --> ', trim(fname_tmp) + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' + + call system_clock(clock1) + if (trim(LAIBCS) == 'GSWP2') then + call process_gswp2_veg (nc,nr,regrid,'grnFrac',n_land, tile_id) + else + if (size(maparc30%ij_index,1) /= 43200) then + ! allocate (maparc30 (1:43200,1:21600)) + call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) + endif + call hres_gswp2 (43200,21600, maparc30, 'green', n_land, tile_lon, tile_lat) + endif + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif - fname_tmp = 'clsm/ndvi.dat' - write (log_file,'(a,a)')' --> ', trim(fname_tmp) - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' - call system_clock(clock1) - call gimms_clim_ndvi (nc,nr, n_land, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif + fname_tmp = 'clsm/lai.dat' + write (log_file,'(a,a)')' --> ', trim(fname_tmp) + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' + redo_modis = .true. + + call system_clock(clock1) + if (trim(LAIBCS) == 'GSWP2') call process_gswp2_veg (nc,nr,regrid,'LAI', n_land, tile_id) + if (trim(LAIBCS) == 'GSWPH') then + if (size(maparc30%ij_index,1) /= 43200) then + ! allocate (maparc30 (1:43200,1:21600)) + call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) + endif + inquire(file='clsm/lai.MODIS_8-DayClim', exist=file_exists) + if (.not.file_exists) call hres_gswp2 (43200,21600, maparc30, 'lai', n_land, tile_lon, tile_lat) + endif + + if (trim(LAIBCS) == 'MODIS') then + lai_name = 'MODIS_8-DayClim/MODIS_' + call hres_lai_no_gswp (43200,21600,maparc30,lai_name, n_land, tile_lon, tile_lon) + endif + + if (trim(LAIBCS) == 'MODGEO') then + lai_name = 'MODIS_8-DayClim/MODIS_' + inquire(file='clsm/lai.MODIS_8-DayClim', exist=file_exists) + if (.not.file_exists)call hres_lai_no_gswp (43200,21600,maparc30,lai_name, n_land, tile_lon, tile_lat, merge=1) + call merge_lai_data (MaskFile, n_land, tile_pfs) + endif + + if (trim(LAIBCS) == 'MODISV6') then + lai_name = 'MCD15A2H.006/MODIS_' + call grid2tile_modis6 (86400,43200,nc,nr,n_land, tile_lon, tile_lat, tile_id, lai_name) + endif + + if (trim(LAIBCS) == 'GLASSA') then + lai_name = 'GLASS-LAI/AVHRR.v4/GLASS01B02.V04.AYYYY' + call grid2tile_glass (nc,nr, tile_id,lai_name, n_land, tile_lon, tile_lat) + endif + + if (trim(LAIBCS) == 'GLASSM') then + lai_name = 'GLASS-LAI/MODIS.v4/GLASS01B01.V04.AYYYY' + call grid2tile_glass (nc,nr,tile_id,lai_name, n_land, tile_lon, tile_lat) + endif + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif - write (log_file,'(a)')' ' - - ! ------------------------------------------------- - - ! call modis_alb_on_tiles (nc,nr,ease_grid,regrid,fnameTil,fnameRst) - ! call modis_scale_para (ease_grid,fnameTil) - ! NOTE: modis_alb_on_tiles uses monthly climatological raster data on 8640x4320 to produce - ! MODIS albedo on tile space. The subroutine was replaced with "modis_alb_on_tiles_high" that process - ! MODIS1 data on native grid and produces 8/16-day MODIS Albedo climatology + fname_tmp = 'clsm/ndvi.dat' + write (log_file,'(a,a)')' --> ', trim(fname_tmp) + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' + call system_clock(clock1) + call gimms_clim_ndvi (nc,nr, n_land, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + ! ------------------------------------------------- + + ! call modis_alb_on_tiles (nc,nr,ease_grid,regrid,fnameTil,fnameRst) + ! call modis_scale_para (ease_grid,fnameTil) + ! NOTE: modis_alb_on_tiles uses monthly climatological raster data on 8640x4320 to produce + ! MODIS albedo on tile space. The subroutine was replaced with "modis_alb_on_tiles_high" that process + ! MODIS1 data on native grid and produces 8/16-day MODIS Albedo climatology + + tmpstring = 'Step 06: Albedo climatologies' + write (log_file,'(a,a,a)') trim(tmpstring),' ', trim(MODALB) + + if(MODALB == 'MODIS1') then + fname_tmp = 'clsm/AlbMap.WS.16-day.tile.0.7_5.0.dat' + write (log_file,'(a,a)')' --> ', trim(fname_tmp) + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + if(F25Tag) then + call create_mapping (nc,nr,21600,10800,maparc60, n_land, tile_id) + call modis_alb_on_tiles_high (21600,10800,maparc60,MODALB, n_land) + deallocate (maparc60%map) + deallocate (maparc60%ij_index) + else + ! This option is for legacy sets like Fortuna 2.1 + call modis_alb_on_tiles (nc,nr,regrid, n_land, tile_id) + endif + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + endif - tmpstring = 'Step 06: Albedo climatologies' - write (log_file,'(a,a,a)') trim(tmpstring),' ', trim(MODALB) - - if(MODALB == 'MODIS1') then - fname_tmp = 'clsm/AlbMap.WS.16-day.tile.0.7_5.0.dat' - write (log_file,'(a,a)')' --> ', trim(fname_tmp) - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - if(F25Tag) then - call create_mapping (nc,nr,21600,10800,maparc60, n_land, tile_id) - call modis_alb_on_tiles_high (21600,10800,maparc60,MODALB, n_land) - deallocate (maparc60%map) - deallocate (maparc60%ij_index) - else - ! This option is for legacy sets like Fortuna 2.1 - call modis_alb_on_tiles (nc,nr,regrid, n_land, tile_id) - endif - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - endif - - if(MODALB == 'MODIS2') then - fname_tmp = 'clsm/AlbMap.WS.8-day.tile.0.3_0.7.dat' - fname_tmp2 = 'clsm/AlbMap.WS.8-day.tile.0.7_5.0.dat' - write (log_file,'(a,a,a,a)')' --> ', trim(fname_tmp), ', ', trim(fname_tmp2) - inquire(file=trim(fname_tmp ), exist=file_exists ) - inquire(file=trim(fname_tmp2), exist=file_exists2) - if ((.not.file_exists).or.(.not.file_exists2)) then - call system_clock(clock1) - write (log_file,'(a)')' Creating files...' - call modis_alb_on_tiles_high (43200,21600,maparc30,MODALB, n_land) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - endif - write (log_file,'(a)')' ' - - ! --------------------------------------------- + if(MODALB == 'MODIS2') then + fname_tmp = 'clsm/AlbMap.WS.8-day.tile.0.3_0.7.dat' + fname_tmp2 = 'clsm/AlbMap.WS.8-day.tile.0.7_5.0.dat' + write (log_file,'(a,a,a,a)')' --> ', trim(fname_tmp), ', ', trim(fname_tmp2) + inquire(file=trim(fname_tmp ), exist=file_exists ) + inquire(file=trim(fname_tmp2), exist=file_exists2) + if ((.not.file_exists).or.(.not.file_exists2)) then + call system_clock(clock1) + write (log_file,'(a)')' Creating files...' + call modis_alb_on_tiles_high (43200,21600,maparc30,MODALB, n_land) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + endif + write (log_file,'(a)')' ' + + ! --------------------------------------------- + + tmpstring = 'Step 07: Albedo scale factors' + write (log_file,'(a,a,a)') trim(tmpstring),' ', trim(MODALB) + + ! NOTE: There are two files with albedo scale factors: "visdf.dat" and "nirdf.dat". + ! Added check for "nirdf.dat", which was missing before. - reichle, 13 Jan 2022 + + fname_tmp = 'clsm/visdf.dat' + fname_tmp2 = 'clsm/nirdf.dat' + write (log_file,'(a,a,a,a)')' --> ', trim(fname_tmp), ', ', trim(fname_tmp2) + inquire(file=trim(fname_tmp ), exist=file_exists ) + inquire(file=trim(fname_tmp2), exist=file_exists2) + if ((redo_modis).or.(.not.file_exists).or.(.not.file_exists2)) then + ! if(.not.F25Tag) then + write (log_file,'(a)')' Creating files... (resolution will be added to file name later)' + call system_clock(clock1) + call modis_scale_para_high (MODALB, n_land) + ! else + ! This option is for legacy sets like Fortuna 2.1 + ! inquire(file='clsm/modis_scale_factor.albvf.clim', exist=file_exists) + ! if ((redo_modis).or.(.not.file_exists)) then + ! call modis_scale_para (ease_grid,fnameTil) + ! call REFORMAT_VEGFILES + ! endif + ! endif + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing files.' + endif + write (log_file,'(a)')' ' - tmpstring = 'Step 07: Albedo scale factors' - write (log_file,'(a,a,a)') trim(tmpstring),' ', trim(MODALB) - - ! NOTE: There are two files with albedo scale factors: "visdf.dat" and "nirdf.dat". - ! Added check for "nirdf.dat", which was missing before. - reichle, 13 Jan 2022 - - fname_tmp = 'clsm/visdf.dat' - fname_tmp2 = 'clsm/nirdf.dat' - write (log_file,'(a,a,a,a)')' --> ', trim(fname_tmp), ', ', trim(fname_tmp2) - inquire(file=trim(fname_tmp ), exist=file_exists ) - inquire(file=trim(fname_tmp2), exist=file_exists2) - if ((redo_modis).or.(.not.file_exists).or.(.not.file_exists2)) then - ! if(.not.F25Tag) then - write (log_file,'(a)')' Creating files... (resolution will be added to file name later)' - call system_clock(clock1) - call modis_scale_para_high (MODALB, n_land) - ! else - ! This option is for legacy sets like Fortuna 2.1 - ! inquire(file='clsm/modis_scale_factor.albvf.clim', exist=file_exists) - ! if ((redo_modis).or.(.not.file_exists)) then - ! call modis_scale_para (ease_grid,fnameTil) - ! call REFORMAT_VEGFILES - ! endif - ! endif - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing files.' - endif - write (log_file,'(a)')' ' - - ! tmpstring1 = '-e EASE -g '//trim(gfile) - ! write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr - ! tmpstring = 'bin/mkCatchParam_openmp '//trim(tmpstring2)//' '//trim(tmpstring1) - -! else - - ! this block is for n_threads>1 - !============================== - - if(trim(SOILBCS)=='NGDC') then - write (log_file,'(a)')'Creating (intermediate) NGDC soil types file...' - call system_clock(clock1) - call create_soil_types_files (nc,nr, n_land, tile_pfs, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - write (log_file,'(a)')' ' - endif - - ! Creating soil_param.first and tau_param.dat files that has 2 options: - ! 1) NGDC soil properties, 2) HWSD-STATSGO2 Soil Properties - ! --------------------------------------------------------------------- - - tmpstring = 'Step 08: Soil parameters ' // trim(SOILBCS) - fname_tmp = 'clsm/soil_param.first' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - if(trim(SOILBCS)=='NGDC') then - if( F25Tag) call soil_para_high (nc,nr,regrid, n_land, tile_id,F25Tag=F25Tag) - if(.not.F25Tag) call soil_para_high (nc,nr,regrid, n_land, tile_id) - endif - if(SOILBCS(1:4)=='HWSD') call soil_para_hwsd (nc,nr, n_land, tile_pfs, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a,a)')' Using existing file.' - endif - write (log_file,'(a)')' ' - - tmpstring = 'Step 09: CLSM model parameters ' // trim(SOILBCS) - fname_tmp = 'clsm/ar.new' - fname_tmp2 = 'clsm/bf.dat' - fname_tmp3 = 'clsm/ts.dat' - fname_tmp4 = 'clsm/soil_param.dat' - tmpstring1 = trim(fname_tmp) // ', ' // trim(fname_tmp2) // ', ' // trim(fname_tmp3) // ', ' // trim(fname_tmp4) - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(tmpstring1), ')' - inquire(file=trim(fname_tmp ), exist=file_exists ) - inquire(file=trim(fname_tmp2), exist=file_exists2) - inquire(file=trim(fname_tmp3), exist=file_exists3) - inquire(file=trim(fname_tmp4), exist=file_exists4) - if ((.not.file_exists).or.(.not.file_exists2).or.(.not.file_exists3).or.(.not.file_exists4)) then - write (log_file,'(a)')' Creating files...' - call system_clock(clock1) - if(trim(SOILBCS)=='NGDC') call create_model_para( MaskFile, n_land, tile_lon, tile_lat, tile_pfs) - if(SOILBCS(1:4) =='HWSD') call create_model_para_woesten(MaskFile, n_land, tile_lon, tile_lat, tile_pfs) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) ' Done. Spent ', seconds, " seconds" - else - write (log_file,'(a,a)')' Using existing files.' - endif - write (log_file,'(a)')' ' - - ! Commented out this call because 7.5-minute raster file is only used - ! for plotting purposes - ! call make_75 (nc,nr,regrid,c_data,fnameRst) - ! write (log_file,'(a)')'Done creating 7.5 minute raster file ......................' - write (log_file,'(a)')'NOTE: 7.5 minute raster file not created (only needed for diagnostic plotting).' - write (log_file,'(a)')' Uncomment associated lines in source to generate 7.5 minute raster file.' - write (log_file,'(a)')' ' - - tmpstring = 'Step 10: CatchCNCLM40 NDep T2m SoilAlb parameters' - fname_tmp = 'clsm/CLM_NDep_SoilAlb_T2m' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - ! create this file only if matching veg types file already exists - inquire(file='clsm/CLM_veg_typs_fracs', exist=file_exists) - if (file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call grid2tile_ndep_t2m_alb (nc,nr, n_land,tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) 'Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')'Skipping step for lack of matching veg types file.' - endif - write (log_file,'(a)')' ' - - tmpstring = 'Step 11: CatchCNCLM45 abm peatf gdp hdm fc parameters' - fname_tmp = 'clsm/CLM4.5_abm_peatf_gdp_hdm_fc' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call CLM45_fixed_parameters (nc,nr, n_land, tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) 'Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' - - tmpstring = 'Step 12: CatchCNCLM45 lightning frequency' - fname_tmp = 'clsm/lnfm.dat' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' - call system_clock(clock1) - call CLM45_clim_parameters (nc,nr,n_land,tile_id) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) 'Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' + ! tmpstring1 = '-e EASE -g '//trim(gfile) + ! write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr + ! tmpstring = 'bin/mkCatchParam_openmp '//trim(tmpstring2)//' '//trim(tmpstring1) - tmpstring = 'Step 13: Country and state codes' - fname_tmp = 'clsm/country_and_state_code.data' - write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' - inquire(file=trim(fname_tmp), exist=file_exists) - if (.not.file_exists) then - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - call map_country_codes (nc,nr,n_land, tile_lon, tile_lat) - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) 'Done. Spent ', seconds, " seconds" - else - write (log_file,'(a)')' Using existing file.' - endif - write (log_file,'(a)')' ' + ! else + + ! this block is for n_threads>1 + !============================== + + if(trim(SOILBCS)=='NGDC') then + write (log_file,'(a)')'Creating (intermediate) NGDC soil types file...' + call system_clock(clock1) + call create_soil_types_files (nc,nr, n_land, tile_pfs, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + write (log_file,'(a)')' ' + endif - if(process_snow_albedo)then - tmpstring = 'Step 14: Static snow albedo from MODIS' - write (log_file,'(a)') trim(tmpstring) - write (log_file,'(a)')' Creating file...' - call system_clock(clock1) - if (trim(SNOWALB)=='MODC061') then - call MODIS_snow_alb (n_land, min_lon, max_lon, min_lat, max_lat) - elseif (trim(SNOWALB)=='MODC061v2') then - if (size(maparc30%ij_index,1) /= 43200) then - call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) - end if - call MODIS_snow_alb_v2(43200,21600,maparc30, n_land) - else - write (log_file,'(a)')'Unknown SNOWALB... stopping!' - stop - endif - call system_clock(clock2) - seconds = (clock2-clock1)/real(clock_rate) - write (log_file, *) 'Done. Spent ', seconds, " seconds" - write (log_file,'(a)')' ' - endif + ! Creating soil_param.first and tau_param.dat files that has 2 options: + ! 1) NGDC soil properties, 2) HWSD-STATSGO2 Soil Properties + ! --------------------------------------------------------------------- + + tmpstring = 'Step 08: Soil parameters ' // trim(SOILBCS) + fname_tmp = 'clsm/soil_param.first' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + if(trim(SOILBCS)=='NGDC') then + if( F25Tag) call soil_para_high (nc,nr,regrid, n_land, tile_id,F25Tag=F25Tag) + if(.not.F25Tag) call soil_para_high (nc,nr,regrid, n_land, tile_id) + endif + if(SOILBCS(1:4)=='HWSD') call soil_para_hwsd (nc,nr, n_land, tile_pfs, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a,a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + tmpstring = 'Step 09: CLSM model parameters ' // trim(SOILBCS) + fname_tmp = 'clsm/ar.new' + fname_tmp2 = 'clsm/bf.dat' + fname_tmp3 = 'clsm/ts.dat' + fname_tmp4 = 'clsm/soil_param.dat' + tmpstring1 = trim(fname_tmp) // ', ' // trim(fname_tmp2) // ', ' // trim(fname_tmp3) // ', ' // trim(fname_tmp4) + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(tmpstring1), ')' + inquire(file=trim(fname_tmp ), exist=file_exists ) + inquire(file=trim(fname_tmp2), exist=file_exists2) + inquire(file=trim(fname_tmp3), exist=file_exists3) + inquire(file=trim(fname_tmp4), exist=file_exists4) + if ((.not.file_exists).or.(.not.file_exists2).or.(.not.file_exists3).or.(.not.file_exists4)) then + write (log_file,'(a)')' Creating files...' + call system_clock(clock1) + if(trim(SOILBCS)=='NGDC') call create_model_para( MaskFile, n_land, tile_lon, tile_lat, tile_pfs) + if(SOILBCS(1:4) =='HWSD') call create_model_para_woesten(MaskFile, n_land, tile_lon, tile_lat, tile_pfs) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) ' Done. Spent ', seconds, " seconds" + else + write (log_file,'(a,a)')' Using existing files.' + endif + write (log_file,'(a)')' ' + + ! Commented out this call because 7.5-minute raster file is only used + ! for plotting purposes + ! call make_75 (nc,nr,regrid,c_data,fnameRst) + ! write (log_file,'(a)')'Done creating 7.5 minute raster file ......................' + write (log_file,'(a)')'NOTE: 7.5 minute raster file not created (only needed for diagnostic plotting).' + write (log_file,'(a)')' Uncomment associated lines in source to generate 7.5 minute raster file.' + write (log_file,'(a)')' ' + + tmpstring = 'Step 10: CatchCNCLM40 NDep T2m SoilAlb parameters' + fname_tmp = 'clsm/CLM_NDep_SoilAlb_T2m' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + ! create this file only if matching veg types file already exists + inquire(file='clsm/CLM_veg_typs_fracs', exist=file_exists) + if (file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call grid2tile_ndep_t2m_alb (nc,nr, n_land,tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) 'Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')'Skipping step for lack of matching veg types file.' + endif + write (log_file,'(a)')' ' + + tmpstring = 'Step 11: CatchCNCLM45 abm peatf gdp hdm fc parameters' + fname_tmp = 'clsm/CLM4.5_abm_peatf_gdp_hdm_fc' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call CLM45_fixed_parameters (nc,nr, n_land, tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) 'Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + tmpstring = 'Step 12: CatchCNCLM45 lightning frequency' + fname_tmp = 'clsm/lnfm.dat' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file... (resolution will be added to file name later)' + call system_clock(clock1) + call CLM45_clim_parameters (nc,nr,n_land,tile_id) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) 'Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + tmpstring = 'Step 13: Country and state codes' + fname_tmp = 'clsm/country_and_state_code.data' + write (log_file,'(a,a,a,a)') trim(tmpstring), ' (', trim(fname_tmp), ')' + inquire(file=trim(fname_tmp), exist=file_exists) + if (.not.file_exists) then + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + call map_country_codes (nc,nr,n_land, tile_lon, tile_lat) + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) 'Done. Spent ', seconds, " seconds" + else + write (log_file,'(a)')' Using existing file.' + endif + write (log_file,'(a)')' ' + + if(process_snow_albedo)then + tmpstring = 'Step 14: Static snow albedo from MODIS' + write (log_file,'(a)') trim(tmpstring) + write (log_file,'(a)')' Creating file...' + call system_clock(clock1) + if (trim(SNOWALB)=='MODC061') then + call MODIS_snow_alb (n_land, min_lon, max_lon, min_lat, max_lat) + elseif (trim(SNOWALB)=='MODC061v2') then + if (size(maparc30%ij_index,1) /= 43200) then + call create_mapping (nc,nr,43200,21600,maparc30, n_land, tile_id) + end if + call MODIS_snow_alb_v2(43200,21600,maparc30, n_land) + else + write (log_file,'(a)')'Unknown SNOWALB... stopping!' + stop + endif + call system_clock(clock2) + seconds = (clock2-clock1)/real(clock_rate) + write (log_file, *) 'Done. Spent ', seconds, " seconds" + write (log_file,'(a)')' ' + endif - ! inquire(file='clsm/irrig.dat', exist=file_exists) - ! if (.not.file_exists) call create_irrig_params (nc,nr,fnameRst) - ! write (log_file,'(a)')'Done computing irrigation model parameters ...............13' - - write (log_file,'(a)')'============================================================' - write (log_file,'(a)')'DONE creating CLSM data files...............................' - write (log_file,'(a)')'============================================================' - write (log_file,'(a)')' ' - - ! call execute_command_line ('chmod 755 bin/create_README.csh ; bin/create_README.csh') -! endif + ! inquire(file='clsm/irrig.dat', exist=file_exists) + ! if (.not.file_exists) call create_irrig_params (nc,nr,fnameRst) + ! write (log_file,'(a)')'Done computing irrigation model parameters ...............13' - close (log_file,status='keep') + write (log_file,'(a)')'============================================================' + write (log_file,'(a)')'DONE creating CLSM data files...............................' + write (log_file,'(a)')'============================================================' + write (log_file,'(a)')' ' + + ! call execute_command_line ('chmod 755 bin/create_README.csh ; bin/create_README.csh') + ! endif + + close (log_file,status='keep') END PROGRAM mkCatchParam + +! ==================== EOF ================================================================ 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 cd716f111..d6c0b3021 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 @@ -29,10 +29,11 @@ PROGRAM mkEASETilesParam ! - removed repetition of identical operations ! - added comments ! - white-space changes for improved readability - use, intrinsic :: iso_fortran_env, only: REAL64 + use, intrinsic :: iso_fortran_env, only: REAL64 use rmTinyCatchParaMod, only : i_raster, j_raster use rmTinyCatchParaMod, only : RegridRasterReal use rmTinyCatchParaMod, only : Target_mean_land_elev + use rmTinyCatchParaMod, only : LakeTopoCat_on_tiles_from_raster, AppendLakeTypeToTileNC4 use process_hres_data, only : histogram use LogRectRasterizeMod, only : SRTM_maxcat, MAPL_UNDEF_R8 ! rasterize.F90 use MAPL_SortMod @@ -112,7 +113,10 @@ PROGRAM mkEASETilesParam character(len=512) :: fname_mask character(len=400) :: MAKE_BCS_INPUT_DIR - ! -------------------------------------------------------------------------------------- + ! LakeTopoCat / ReachTopoCat LakeType + integer, allocatable :: tile_lake_type(:) + integer :: rc_lake + ! ---------------------------------------------- call get_environment_variable( "MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR ) @@ -976,7 +980,37 @@ PROGRAM mkEASETilesParam call MAPL_WriteTilingNC4('til/'//trim(gfile)//'.nc4', [EASELabel],[nc_ease],[nr_ease], & nc, nr, iTable, rTable) + + ! GEOS lake tiles include coastal ocean surfaces such as fjords and estuaries + ! + ! To distinguish between these surfaces and "real" lakes (inland water), add + ! a tile-space flag (field) to the nc4 tile file with information based on LakeTopoCat v1.1 data. + ! + ! LakeTopoCat / ReachTopoCat: + ! compute encoded LakeType in tile space and append to nc4 tile file. + ! + ! tile_lake_type: + ! -9999 = UNDEF / excluded, e.g. typ==100 + ! 0 = no LakeTopoCat lake or ReachTopoCat reach touch + ! 1 = lake touch only + ! 2 = reach touch only + ! 3 = lake + reach touch + ! + ! This runs only with nc=43200, nr=21600 for GEOS5_10arcsec_mask workflows; + ! returns rc_lake>0 otherwise. + + allocate(tile_lake_type(n_landlakelandice)) + call LakeTopoCat_on_tiles_from_raster(n_landlakelandice, nc, nr, tileid_index, & + iTable(1:n_landlakelandice,0), tile_lake_type, rc_lake) + + if (rc_lake == 0) call AppendLakeTypeToTileNC4( & + 'til/'//trim(gfile)//'.nc4', n_landlakelandice, tile_lake_type) + + deallocate(tile_lake_type) + + ! -------------------------------------------- + deallocate( tileid_index, catid_index,veg ) deallocate( tile_area, ease_grid_area, tile_elev, my_land, all_id ) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 index b3df41f23..c31b2ead0 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/rmTinyCatchParaMod.F90 @@ -2,8 +2,6 @@ #define ASSERT_(A) if(.not.A)then;print *,'Error:',__FILE__,__LINE__;stop;endif ! ! A Collection subroutines needed by mkCatchParam.F90 -! Contact: Sarith Mahanama sarith.p.mahanama@nasa.gov -! Email : sarith.p.mahanama@nasa.gov module rmTinyCatchParaMod @@ -12,8 +10,10 @@ module rmTinyCatchParaMod use MAPL_Base, ONLY: MAPL_UNDEF use MAPL, only: MAPL_WriteTilingNC4, MAPL_ease_extent use lsm_routines, ONLY: sibalb - use LogRectRasterizeMod, ONLY: SRTM_maxcat, MAPL_UNDEF_R8 + use LogRectRasterizeMod, ONLY: SRTM_maxcat, MAPL_UNDEF_R8 + use, intrinsic :: iso_fortran_env, only: REAL64 + implicit none logical, parameter :: error_file=.true. @@ -31,7 +31,7 @@ module rmTinyCatchParaMod logical, parameter :: bug =.false. include 'netcdf.inc' - + logical :: preserve_soiltype = .false. ! Bugfix for Target_mean_land_elev: @@ -55,6 +55,7 @@ module rmTinyCatchParaMod public REFORMAT_VEGFILES public Get_MidTime, Time_Interp_Fac public ascat_r0, jpl_canoph, NC_VarID, init_bcs_config + public LakeTopoCat_on_tiles_from_raster, AppendLakeTypeToTileNC4 ! The following variables define the details of the BCs versions (data sources). ! Initialize to dummy values here and set to desired values in init_bcs_config(). @@ -924,6 +925,279 @@ END SUBROUTINE modis_alb_on_tiles ! subroutine pick_cat(sam,clr) ! !---------------------------------------------------------------------- + !---------------------------------------------------------------------- + !---------------------------------------------------------------------- + SUBROUTINE LakeTopoCat_on_tiles_from_raster( n_tile, nc_rst, nr_rst, tile_id, & + tile_typ, tile_lake_type, rc) + + ! GEOS lake tiles include coastal ocean surfaces such as fjords and estuaries + ! + ! To distinguish between these surfaces and "real" lakes (inland water), create + ! a tile-space flag (field) to the nc4 tile file with information based on LakeTopoCat v1.1 data. + + ! Map preprocessed LakeTopoCat / ReachTopoCat any-touch support + ! (30 arcsec resolution, 10 deg x 10 deg chunks) + ! to tile space using the 30 arcsec raster tile-id map. + ! + ! Why NOT use rmap/create_mapping here? + ! ------------------------------------ + ! 1) In mkCatchParam, the "tile creation" step (Step 01) happens BEFORE any + ! precomputed rmap exists. Those mappings are built later for other products. + ! 2) create_mapping() is designed for land tiles only (1..n_land), + ! which is appropriate for land-only products such as MODIS snow albedo, + ! but this LakeTopoCat / ReachTopoCat diagnostic must be computed for + ! all candidate tile types in the tile file, including lake-like and + ! ocean tiles where present. + ! 3) Lake/reach support inputs are already on the same global 30-arcsec + ! grid as the raster tile-id map (43200x21600), so each 30" pixel can + ! be mapped directly to a tile via tile_id(iG,jG). For this diagnostic + ! we only need any-touch support, not area-weighted remapping. + ! + ! Output tile_lake_type: + ! -9999 = UNDEF / excluded, e.g. land typ==100 or landice typ==20 + ! 0 = no LakeTopoCat lake or ReachTopoCat reach touch + ! 1 = lake touch only + ! 2 = reach touch only + ! 3 = lake + reach touch + ! + + implicit none + + integer, intent(in) :: n_tile + integer, intent(in) :: nc_rst, nr_rst + integer, intent(in) :: tile_id(1:nc_rst, 1:nr_rst) + integer, intent(in) :: tile_typ(1:n_tile) + + integer, intent(out) :: tile_lake_type(1:n_tile) + integer, intent(out), optional :: rc + + integer, parameter :: LAKETYPE_UNDEF = -9999 + + integer(kind=4), parameter :: nc_10 = 1200 + integer(kind=4), parameter :: nr_10 = 1200 + + integer :: status, ncid, varid + integer :: ix, jx, ii, jj + integer :: iLL, jLL, iG, jG, tid + character*512 :: lake_fname, reach_fname + character*2 :: HH, VV + character*512 :: MAKE_BCS_INPUT_DIR + + integer, allocatable :: lake_any(:,:) ! 10 deg x 10 deg chunk, 30 arcsec + integer, allocatable :: reach_any(:,:) ! 10 deg x 10 deg chunk, 30 arcsec + + call get_environment_variable("MAKE_BCS_INPUT_DIR", MAKE_BCS_INPUT_DIR) + + ! ------------------------------------------------------------------ + ! IMPORTANT: + ! LakeTopoCat / ReachTopoCat inputs are preprocessed onto a + ! 30-arcsec raster grid (43200 x 21600). This routine assumes that + ! the tile_id raster (nc_rst x nr_rst) is also on that same 30 arcsec grid. + ! + ! For workflows using GEOS5_10arcsec_mask, nx=43200 and ny=21600, + ! lake/reach presence can be mapped directly via tile_id(iG,jG). + ! + ! For coarser or alternative masks (e.g., 8640x4320), there is no + ! consistent 30" raster-to-tile alignment here. In those cases we + ! skip LakeTopoCat / ReachTopoCat generation rather than attempting + ! an implicit remap that could introduce silent errors. + ! + ! tile_lake_type is therefore only populated when the raster resolution + ! matches the 30 arcsec grid. If skipped, tile_lake_type remains UNDEF. + ! ------------------------------------------------------------------ + + ! Start as UNDEF. If we skip generation, do NOT report candidate tiles as + ! LakeType=0, because that would incorrectly mean "checked and no touch". + + tile_lake_type(:) = LAKETYPE_UNDEF + + if (nc_rst /= 43200 .or. nr_rst /= 21600) then + print *, 'NOTE: Skipping LakeTopoCat / ReachTopoCat LakeType generation.' + print *, ' Requires 43200x21600 raster. Got ', nc_rst, nr_rst + if (present(rc)) rc = 1 + return + endif + + ! Define LakeType only for candidate MAPL tile types (water surfaces). Set default to 0 (no touch). + do tid = 1, n_tile + if (tile_typ(tid) == 0 .or. tile_typ(tid) == 19 ) then + tile_lake_type(tid) = 0 + endif + enddo + + allocate(lake_any (1:nc_10, 1:nr_10)) + allocate(reach_any(1:nc_10, 1:nr_10)) + + ! Loop through 36x18 10-degree chunks. + do jx = 1, 18 + do ix = 1, 36 + + write(HH,'(i2.2)') ix + write(VV,'(i2.2)') jx + + lake_fname = trim(MAKE_BCS_INPUT_DIR)// & + '/lake/lake_mask/v1/LakeTopoCat_30arcsec_H'//HH//'V'//VV//'.nc4' + + reach_fname = trim(MAKE_BCS_INPUT_DIR)// & + '/lake/reach_mask/v1/ReachTopoCat_30arcsec_H'//HH//'V'//VV//'.nc4' + + ! Read lake_presence_any. + status = NF_OPEN(trim(lake_fname), NF_NOWRITE, ncid) + if (status /= NF_NOERR) then + print *, 'ERROR: Lake mask tile file not found: ', trim(lake_fname) + print *, 'STOPPING.' + stop + endif + + status = NF_INQ_VARID(ncid, 'lake_presence_any', varid) ; VERIFY_(status) + status = NF_GET_VARA_INT(ncid, varid, (/1,1/), (/nc_10,nr_10/), lake_any) ; VERIFY_(status) + status = NF_CLOSE(ncid) ; VERIFY_(status) + + ! Read reach_presence_any. + status = NF_OPEN(trim(reach_fname), NF_NOWRITE, ncid) + if (status /= NF_NOERR) then + print *, 'ERROR: Reach mask tile file not found: ', trim(reach_fname) + print *, 'STOPPING.' + stop + endif + + status = NF_INQ_VARID(ncid, 'reach_presence_any', varid) ; VERIFY_(status) + status = NF_GET_VARA_INT(ncid, varid, (/1,1/), (/nc_10,nr_10/), reach_any) ; VERIFY_(status) + status = NF_CLOSE(ncid) ; VERIFY_(status) + + ! Lower-left global indices of this 10 deg x 10 deg chunk, 1-based. + iLL = (ix-1)*nc_10 + 1 + jLL = (jx-1)*nr_10 + 1 + + do jj = 1, nr_10 + do ii = 1, nc_10 + + iG = ii + iLL - 1 + jG = jj + jLL - 1 + + tid = tile_id(iG, jG) + + if (tid < 1 .or. tid > n_tile) cycle + + ! Only candidate tile types get LakeType. + if (tile_lake_type(tid) == LAKETYPE_UNDEF) cycle + + ! Bit 1: lake touch. + if (lake_any(ii,jj) > 0) then + tile_lake_type(tid) = IOR(tile_lake_type(tid), 1) + endif + + ! Bit 2: reach touch. + if (reach_any(ii,jj) > 0) then + tile_lake_type(tid) = IOR(tile_lake_type(tid), 2) + endif + + enddo + enddo + + enddo + enddo + + if (present(rc)) rc = 0 + + if (allocated(lake_any)) deallocate(lake_any) + if (allocated(reach_any)) deallocate(reach_any) + + END SUBROUTINE LakeTopoCat_on_tiles_from_raster + + !---------------------------------------------------------------------- + + SUBROUTINE AppendLakeTypeToTileNC4( tilefile, n_tile, lake_type ) + + implicit none + + character(*), intent(in) :: tilefile + integer, intent(in) :: n_tile + integer, intent(in) :: lake_type(1:n_tile) + + integer, parameter :: LAKETYPE_UNDEF = -9999 + + integer :: status, ncid, ndims, nvars, ngatts, unlimdimid + integer :: dimid_tile, dimlen, dd + integer :: varid_lake_type + integer :: fill_value(1) + integer :: flag_values(5) + + character(len=NF_MAX_NAME) :: dname + character(len=256) :: mylongname + character(len=256) :: flagmeanings + character(len=512) :: comment + + fill_value(1) = LAKETYPE_UNDEF + flag_values = (/ LAKETYPE_UNDEF, 0, 1, 2, 3 /) + + ! append data into nc4 tile file + + status = NF_OPEN(trim(tilefile), NF_WRITE, ncid) ; VERIFY_(status) + + ! find tile dimension by matching length + + status = NF_INQ(ncid, ndims, nvars, ngatts, unlimdimid) ; VERIFY_(status) + + dimid_tile = -1 + do dd = 1, ndims + status = NF_INQ_DIM(ncid, dd, dname, dimlen) ; VERIFY_(status) + if (dimlen == n_tile) then + dimid_tile = dd + exit + endif + enddo + + if (dimid_tile < 0) then + print *, 'ERROR: could not find tile dimension of length ', n_tile, ' in ', trim(tilefile) + stop + endif + + ! Define lake_type in the freshly written nc4 tile file. + ! In the normal make_bcs workflow this variable should not already exist; + ! if it does, NF_DEF_VAR will fail rather than silently replacing data. + + status = NF_REDEF(ncid) ; VERIFY_(status) + + status = NF_DEF_VAR(ncid, 'lake_type', NF_INT, 1, (/dimid_tile/), & + varid_lake_type) ; VERIFY_(status) + + ! _FillValue should be defined before leaving define mode. + + status = NF_PUT_ATT_INT(ncid, varid_lake_type, '_FillValue', NF_INT, 1, & + fill_value) ; VERIFY_(status) + + status = NF_PUT_ATT_INT(ncid, varid_lake_type, 'missing_value', NF_INT, 1, & + fill_value) ; VERIFY_(status) + + mylongname = 'flag identifying intersection of lake tile with a lake polygon or reach line from LakeTopoCat v1.1 data' + status = NF_PUT_ATT_TEXT(ncid, varid_lake_type, 'long_name', & + len_trim(mylongname), mylongname) ; VERIFY_(status) + + status = NF_PUT_ATT_TEXT(ncid, varid_lake_type, 'units', 1, '1') ; VERIFY_(status) + + status = NF_PUT_ATT_INT(ncid, varid_lake_type, 'flag_values', NF_INT, 5, & + flag_values) ; VERIFY_(status) + + flagmeanings = 'undefined, no_lake_or_reach_intersection, lake_intersection_only, reach_intersection_only, lake_and_reach_intersection' + status = NF_PUT_ATT_TEXT(ncid, varid_lake_type, 'flag_meanings', & + len_trim(flagmeanings), flagmeanings) ; VERIFY_(status) + + comment = 'Defined for typ==0 (ocean) and typ==19 (lake); set to -9999 otherwise (land or landice).' + status = NF_PUT_ATT_TEXT(ncid, varid_lake_type, 'comment', & + len_trim(comment), comment) ; VERIFY_(status) + + status = NF_ENDDEF(ncid) ; VERIFY_(status) + + ! Write lake_type data and close file. + + status = NF_PUT_VARA_INT(ncid, varid_lake_type, (/1/), (/n_tile/), & + lake_type) ; VERIFY_(status) + + status = NF_CLOSE(ncid) ; VERIFY_(status) + + END SUBROUTINE AppendLakeTypeToTileNC4 + !---------------------------------------------------------------------- SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id) @@ -980,6 +1254,10 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id) integer :: ip_keep, k, nc_tmp, nr_tmp real :: EASE_LAT_MAX ! + ! LakeTopoCat / ReachTopoCat LakeType + integer, allocatable :: lake_type(:) + integer :: rc_lake + ! ----------------------------------------------------- ! ! get elevation (q0) from "gtopo30" raster file ("srtm30_withKMS_2.5x2.5min.data") @@ -1254,6 +1532,7 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id) endwhere fname=trim(fnameTil)//'.nc4' + if (two_EASE) then ! remove tiles outside EASE grid domain call MAPL_ease_extent(gName(1), nc_tmp, nr_tmp, ur_lat = EASE_LAT_MAX) @@ -1274,19 +1553,44 @@ SUBROUTINE supplemental_tile_attributes(nx,ny,regrid,dateline,fnameTil, Rst_id) endif enddo - call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), & - nx, ny, iTable_keep, rTable_keep, N_PfafCat=SRTM_maxcat, rc=status) + call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), & + nx, ny, iTable_keep, rTable_keep, N_PfafCat=SRTM_maxcat, rc=status) deallocate(iTable_keep, rTable_keep, keep_tile) else - call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), & - nx, ny, iTable, rTable, N_PfafCat=SRTM_maxcat, rc=status) + call MAPL_WriteTilingNC4(fname, gName(1:n_grid), im(1:n_grid), jm(1:n_grid), & + nx, ny, iTable, rTable, N_PfafCat=SRTM_maxcat, rc=status) endif + ! LakeTopoCat / ReachTopoCat: + ! compute encoded LakeType in tile space and append to nc4 tile file. + ! + ! lake_type: + ! -9999 = UNDEF / excluded, e.g. typ==100 + ! 0 = no LakeTopoCat lake or ReachTopoCat reach touch + ! 1 = lake touch only + ! 2 = reach touch only + ! 3 = lake + reach touch + ! + ! This runs only with nx=43200, ny=21600 for GEOS5_10arcsec_mask workflows; + ! returns rc_lake>0 otherwise. + + allocate(lake_type(ip)) ! ip = n_tile + + call LakeTopoCat_on_tiles_from_raster(ip, nx, ny, Rst_id, & + iTable(1:ip,0), lake_type, rc_lake) + + if (rc_lake == 0) call AppendLakeTypeToTileNC4(fname, ip, lake_type) + + deallocate(lake_type) + + ! -------------------------------------------- + deallocate (rTable, iTable) deallocate (limits) deallocate (catid) deallocate (q0) + if(regrid) then deallocate(raster) endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/README.md b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/README.md new file mode 100755 index 000000000..a07450b6f --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/README.md @@ -0,0 +1,224 @@ +## LakeTopoCat Lake and Reach Mask Preprocessing for GEOS BCS + +**February 2026** + +## Overview + +This repository contains preprocessing scripts that convert **HydroLAKES-TopoCat v1.1 (2023)** global **lake polygons** and **reach line features** into **30-arcsecond raster products** compatible with `make_bcs`. + +The workflow produces: + +- A native global **30″ raster** (43200 × 21600) +- 36 × 18 GEOS **10°×10° tiles** (1200 × 1200 each) +- Separate tiled products for **lakes** and **reaches** +- Tile-level lake and reach information that can be read later during GEOS tile aggregation + +The output format is consistent with other GEOS static datasets (soil properties, snow albedo, etc.). + +## Workflow Summary + +## Step 1 — Rasterization & Aggregation + +### Lakes + +**Script:** `preproc_lake_to_30arcsec.py` + +#### Method + +- Rasterize HydroLAKES-TopoCat lake polygons at **10 arc-second resolution** +- Aggregate 3×3 blocks → **30 arc-second grid** +- Process in **5° latitude bands** +- Combine bands into global file: + `LakeTopoCat_Global_30arcsec.nc4` + +#### Why 10″ → 30″? + +Oversampling at 10″ improves shoreline representation and allows estimation of **fractional lake coverage** within each 30″ pixel. + +Direct 30″ rasterization would produce only a binary mask. + +## Global Lake Output File + +**Dimensions** + +- 43200 (lon) × 21600 (lat) +- Resolution: 30 arc-seconds + +### Variables + +| Variable | Type | Description | +|----------|------|-------------| +| `lake_area_frac` | float32 | Fractional lake coverage in each 30″ cell [0–1] | +| `lake_presence_any` | uint8 | Binary lake presence (0/1) | + +### Data Completeness + +- Data are complete everywhere (no NaNs). +- `_FillValue` attributes exist but should not appear in valid data. + +## Step 1b — Reach Rasterization & Aggregation + +### Reaches + +**Script:** `preproc_reach_to_30arcsec.py` + +#### Method + +- Rasterize HydroLAKES-TopoCat reach line features at **10 arc-second resolution** +- Use touch-based rasterization for line geometry +- Aggregate 3×3 blocks → **30 arc-second grid** +- Process in **5° latitude bands** +- Combine bands into global file: + `ReachTopoCat_Global_30arcsec.nc4` + +#### Why a separate reach product? + +Reaches are **linear features**, not polygons, so they do not represent true area coverage in the same way as lakes. + +For this reason, the reach product is stored separately from the lake product and is intended for later experimentation during GEOS tile aggregation. + +## Global Reach Output File + +**Dimensions** + +- 43200 (lon) × 21600 (lat) +- Resolution: 30 arc-seconds + +### Variables + +| Variable | Type | Description | +|----------|------|-------------| +| `reach_occupancy_frac` | float32 | Fraction of 10″ subpixels within each 30″ cell touched by a reach [0–1] | +| `reach_presence_any` | uint8 | Binary reach presence (0/1) | + +## Step 2 — Split into GEOS 10°×10° Tiles + +**Scripts:** `split_lake_30arcsec_to_HV_generic.py` + +The configuration block is switched between lake and reach inputs. + +Produces: + +- `LakeTopoCat_30arcsec_H##V##.nc4` +- `ReachTopoCat_30arcsec_H##V##.nc4` + +### Per-Tile Properties + +For both lake and reach tile products: + +- 1200 × 1200 pixels +- 30 arc-second resolution +- 1D coordinates: + - `latitude(N_lat)` + - `longitude(N_lon)` + +### Per-Tile Variables + +**Lake tiles** (`LakeTopoCat_30arcsec_H##V##.nc4`) + +| Variable | Type | Description | +|----------|------|-------------| +| `lake_area_frac` | float32 | Fractional lake coverage in each 30″ cell [0–1] | +| `lake_presence_any` | uint8 | Binary lake presence (0/1) | + +**Reach tiles** (`ReachTopoCat_30arcsec_H##V##.nc4`) + +| Variable | Type | Description | +|----------|------|-------------| +| `reach_occupancy_frac` | float32 | Fraction of 10″ subpixels within each 30″ cell touched by a reach [0–1] | +| `reach_presence_any` | uint8 | Binary reach presence (0/1) | + +### Global Attributes Included + +- `N_lon_global = 43200` +- `N_lat_global = 21600` +- `i_ind_offset_LL` +- `j_ind_offset_LL` +- `CellSize_arc_Secs = 30` + +## Source Datasets + +**HydroLAKES-TopoCat v1.1 (2023)** + +- Source product is vector data, not a regular raster grid. +- Constructed from the HydroLAKES v1.0 lake mask and MERIT Hydro v1.0.1 hydrography. +- MERIT Hydro v1.0.1 is a 3 arc-second hydrography dataset, approximately 90 m at the equator. +- CRS: EPSG:4326 + +### Lake input + +- Input: `Lakes_pfaf_*.shp` +- Geometry used: polygon geometry only +- HydroLAKES includes global lake/reservoir shoreline polygons for lakes with surface area of at least 10 ha. +- No filtering by lake size, permanence, or type is applied in this preprocessing workflow. + +### Reach input + +- Input: `Reaches_pfaf_*.shp` +- Geometry used: line geometry only +- Reaches are rasterized as touch-based linear features, not polygon area + +### GEOS rasterization used here + +The GEOS preprocessing does not use the source vector data directly in `make_bcs`. +Instead, the LakeTopoCat / ReachTopoCat vector features are rasterized first at 10 arc-second resolution and then aggregated to the 30 arc-second products used by the Fortran tile aggregation. + +## GEOS Usage (mkCatchParam Step 01) + +### Current implementation + +During `mkCatchParam` Step 01, the LakeTopoCat and ReachTopoCat products are mapped from the 30 arcsecond H/V raster files to GEOS tile space using the raster `tile_id` grid. + +The Fortran aggregation reads only the binary any-touch fields: + +- `lake_presence_any` +- `reach_presence_any` + +No interpolation or separate spatial remapping is used. Each 30″ raster cell is assigned directly to a GEOS tile using `tile_id(iG,jG)`. + +This processing is enabled only when the raster tile-id grid is `43200 × 21600`. +For coarser or alternative masks, LakeTopoCat / ReachTopoCat tile aggregation is skipped to avoid an implicit remap. + + +### Tile-Level Variable Created + +| Variable | Type | Description | +|----------|------|-------------| +| `lake_type` | int | Encoded LakeTopoCat / ReachTopoCat intersection type for candidate water tiles | + +### `lake_type` Coding + +| Value | Meaning | +|-------|---------| +| `-9999` | undefined / excluded tile | +| `0` | candidate water tile with no LakeTopoCat lake intersection and no ReachTopoCat reach intersection | +| `1` | LakeTopoCat lake intersection only | +| `2` | ReachTopoCat reach intersection only | +| `3` | both LakeTopoCat lake and ReachTopoCat reach intersection | + +### Candidate Tile Types + +`lake_type` is defined only for candidate water tile types: + +- `typ == 0` — ocean +- `typ == 19` — lake + +All other tile types, including land and land ice, are assigned `-9999`. + +For the current standard EASE file, `typ == 0` is absent, so the defined candidate space is effectively `typ == 19`. For CF and future tile files, `typ == 0` may be present and will use the same coding. + +### Notes + +The intermediate preprocessing products still include fractional/occupancy diagnostics: + +- `lake_area_frac` +- `reach_occupancy_frac` + +These are retained for QA and diagnostics, but the Fortran `lake_type` product is based on the binary any-touch fields: + +- `lake_presence_any` +- `reach_presence_any` + +The variable is appended during supplemental tile attribute generation when the 30 arcsecond raster tile-id grid is available. +In the normal `make_bcs` workflow, `lake_type` is created once in the freshly written nc4 tile file; rerunning on an already modified tile file is not expected. + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/preproc_lake_to_30arcsec.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/preproc_lake_to_30arcsec.py new file mode 100755 index 000000000..c64d99020 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/preproc_lake_to_30arcsec.py @@ -0,0 +1,227 @@ +#!/usr/bin/env python3 +""" +Create HydroLAKES-based lake mask for GEOS/BCS tiling from Lake-TopoCat (HydroLAKES-TopoCat v1.1 2023): + +Output (the only file makebcs Fortran needs): + LakeTopoCat_Global_30arcsec.nc4 + +Variables: + - lake_area_frac(lat, lon): float32 in [0, 1] + Fractional lake coverage of each 30 arc-sec cell, estimated by oversampling + at 10 arc-sec and averaging 3x3 subpixels. + - lake_presence_any(lat, lon): uint8 0/1 (optional but cheap and handy) + 1 if any 10" subpixel in the 30" cell is lake. + +Method: + 1) Build 10 arc-sec binary raster in memory for a latitude band (5 degrees). + 2) Immediately aggregate 10" -> 30" in that band (3x3 mean/max). + 3) Write 30" band NetCDF. + 4) Combine all 30" bands into one global NetCDF. + +Why oversample at 10"? + Rasterization at 30" directly gives only a binary pixel mask. Oversampling at 10" + lets us estimate fractional coverage in each 30" pixel at low cost. +""" + +import os +import glob +import numpy as np +import geopandas as gpd +import rasterio.features +import xarray as xr +from affine import Affine + +# ------------------------------------------------------------------ +# Input lake polygon dataset +# +# Source: HydroLAKES-TopoCat v1.1 Y2023 processing) +# Files: Lakes_pfaf_*.shp +# +# Description: +# Global lake polygon dataset (EPSG:4326) partitioned by Pfafstetter +# basin code. Each shapefile contains HydroLAKES polygons with attributes +# such as lake area, lake type, and permanence. +# +# Usage in this script: +# Only polygon geometry is used. +# All polygons are rasterized as lake presence (value = 1). +# No filtering by lake type, size, or permanence is applied. +# ------------------------------------------------------------ +# Paths (edit as needed) +# ------------------------------------------------------------------ +shapefile_dir = "/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/lake/Lake_TopoCat/v1/Lakes/" +out_dir = "/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/lake/Lake_TopoCat/v1/lake_mask_build/" +out_global = "/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/lake/lake_mask/v1/before_splitting_nc4/LakeTopoCat_Global_30arcsec.nc4" + + +os.makedirs(out_dir, exist_ok=True) + +# ------------------------------------------------------------ +# Grid definitions +# ------------------------------------------------------------ +# Oversampling grid: 10 arc-sec +res10 = 10 / 3600.0 +lon10 = np.arange(-180 + res10/2, 180, res10) # 129600 +lat10 = np.arange(-90 + res10/2, 90, res10) # 64800 + +# Target grid: 30 arc-sec = 3 x 10 arc-sec +lon30 = lon10[1::3] # 43200 +# lat30 varies by band; global lat30 would be lat10[1::3] = 21600 + +# Process in latitude bands to limit memory +band_deg = 5.0 + +# ------------------------------------------------------------ +# Shapefile list +# ------------------------------------------------------------ +shp_list = sorted(glob.glob(os.path.join(shapefile_dir, "Lakes_pfaf_*.shp"))) +if not shp_list: + raise RuntimeError(f"No shapefiles found: {shapefile_dir}/Lakes_pfaf_*.shp") +print(f"Found {len(shp_list)} shapefiles") + +# ------------------------------------------------------------ +# Helper: 10" -> 30" aggregation (3x3 blocks) +# ------------------------------------------------------------ +def agg10_to_30(binary10: np.ndarray): + """ + binary10: uint8 array (nlat10, nlon10) with values 0/1 + Returns: + frac30: float32 (nlat30, nlon30) = mean of 3x3 block + any30 : uint8 (nlat30, nlon30) = max of 3x3 block + """ + nlat, nlon = binary10.shape + nlat3 = (nlat // 3) * 3 + nlon3 = (nlon // 3) * 3 + binary10 = binary10[:nlat3, :nlon3] + + # reshape into 3x3 blocks + blk = binary10.reshape(nlat3//3, 3, nlon3//3, 3) + frac30 = blk.mean(axis=(1, 3)).astype(np.float32) + any30 = blk.max(axis=(1, 3)).astype(np.uint8) + return frac30, any30 + +# ------------------------------------------------------------ +# Build 30" band files (each band is small enough to write fast) +# ------------------------------------------------------------ +band_files = [] + +for lat_start in np.arange(-90, 90, band_deg): + lat_end = lat_start + band_deg + print(f"\n=== Latitude band {lat_start} to {lat_end} ===") + + # 10" lat centers in this band + lat10_band = lat10[(lat10 >= lat_start) & (lat10 < lat_end)] + if lat10_band.size == 0: + continue + + out_shape10 = (lat10_band.size, lon10.size) + + # Affine transform for this band (pixel -> lon/lat) + transform10 = ( + Affine.translation(lon10[0] - res10/2, lat10_band[0] - res10/2) + * Affine.scale(res10, res10) + ) + + # In-memory 10" band raster (binary 0/1) + lake10 = np.zeros(out_shape10, dtype=np.uint8) + + # Rasterize all shapefiles, but only features intersecting this lat band + for shp in shp_list: + gdf = gpd.read_file(shp).to_crs("EPSG:4326") + gdf_band = gdf.cx[:, lat_start:lat_end] + if gdf_band.empty: + continue + + shapes = ((geom, 1) for geom in gdf_band.geometry) + + tmp = rasterio.features.rasterize( + shapes, + out_shape=out_shape10, + transform=transform10, + fill=0, + dtype=np.uint8, + ) + + # Union across shapefiles + lake10 = np.maximum(lake10, tmp) + + # Aggregate to 30" for this band + frac30, any30 = agg10_to_30(lake10) + + # 30" lat centers for this band (aligned with 10" centers) + lat30_band = lat10_band[1::3] + + # Sanity check dimensions + if frac30.shape != (lat30_band.size, lon30.size): + raise RuntimeError(f"Shape mismatch: frac30={frac30.shape} vs coords={(lat30_band.size, lon30.size)}") + + # Write band netcdf + ds = xr.Dataset( + { + "lake_area_frac": (("lat", "lon"), frac30), + "lake_presence_any": (("lat", "lon"), any30), # optional but useful + }, + coords={"lat": lat30_band, "lon": lon30}, + ) + + band_out = os.path.join(out_dir, f"LakeTopoCat_30arcsec_{lat_start}_{lat_end}.nc4") + encoding = { + "lake_area_frac": {"zlib": True, "complevel": 4, "dtype": "float32", "_FillValue": -9999.0}, + "lake_presence_any": {"zlib": True, "complevel": 4, "dtype": "uint8", "_FillValue": 255}, + } + + ds.to_netcdf(band_out, format="NETCDF4", engine="netcdf4", encoding=encoding) + print("Wrote band:", band_out) + band_files.append(band_out) + +# ------------------------------------------------------------ +# Combine band files into the single global 30" product +# ------------------------------------------------------------ +print(f"\nCombining {len(band_files)} 30\" band files into:\n {out_global}") + +ds = xr.open_mfdataset( + sorted(band_files), + combine="by_coords", + parallel=False +).sortby("lat") + +# Force arrays into memory to avoid dask/lazy-write issues +frac = ds["lake_area_frac"].load().astype("float32") +anyv = ds["lake_presence_any"].load().astype("uint8") + +ds_out = xr.Dataset( + { + "lake_area_frac": frac, + "lake_presence_any": anyv, + }, + coords={ + "lat": ds["lat"].values, + "lon": ds["lon"].values, + }, +) + +encoding = { + "lake_area_frac": { + "zlib": True, + "complevel": 4, + "dtype": "float32", + "_FillValue": -9999.0, + }, + "lake_presence_any": { + "zlib": True, + "complevel": 4, + "dtype": "uint8", + "_FillValue": 255, + }, +} + +ds_out.to_netcdf( + out_global, + format="NETCDF4", + engine="netcdf4", + encoding=encoding, +) + +print("Done. Global file written:") +print(" ", out_global) + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/preproc_reach_to_30arcsec.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/preproc_reach_to_30arcsec.py new file mode 100755 index 000000000..97807fada --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/preproc_reach_to_30arcsec.py @@ -0,0 +1,265 @@ +""" +Create HydroLAKES-TopoCat reach mask for GEOS/BCS tiling from Lake-TopoCat Reaches. + +Output: + ReachTopoCat_Global_30arcsec.nc4 + +Variables: + - reach_occupancy_frac(lat, lon): float32 in [0, 1] + Fraction of 10 arc-sec subpixels touched by any reach geometry within + the 30 arc-sec cell. This is a line-occupancy metric, NOT true channel area. + - reach_presence_any(lat, lon): uint8 0/1 + 1 if any 10 arc-sec subpixel in the 30 arc-sec cell is touched by a reach. + +Method: + 1) Build 10 arc-sec binary raster in memory for a latitude band (5 degrees). + 2) Rasterize reach lines with all_touched=True to preserve thin features. + 3) Aggregate 10" -> 30" in that band (3x3 mean/max). + 4) Write 30" band NetCDF. + 5) Combine all 30" bands into one global NetCDF. + +Notes: + This script is intentionally separate from the lake preprocessor because + reaches are line geometries, not polygons. + reach_occupancy_frac is useful as a density-style metric, but it should not + be interpreted as water area fraction. +""" + +import os +import glob +import numpy as np +import geopandas as gpd +import rasterio.features +import xarray as xr +from affine import Affine + +# ------------------------------------------------------------------ +# Input reach line dataset +# +# Source: HydroLAKES-TopoCat v1.1 (2023) +# Files: Reaches_pfaf_*.shp +# +# Usage in this script: +# Only line geometry is used. +# All reaches are rasterized as binary reach presence (value = 1). +# ------------------------------------------------------------------ +# Paths (edit as needed) +# ------------------------------------------------------------------ +shapefile_dir = "/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/lake/Lake_TopoCat/v1/Reaches/" +out_dir = "/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/lake/Lake_TopoCat/v1/reach_mask_build/" +out_global = "/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/lake/reach_mask/v1/before_splitting_nc4/ReachTopoCat_Global_30arcsec.nc4" + +os.makedirs(out_dir, exist_ok=True) +os.makedirs(os.path.dirname(out_global), exist_ok=True) + +# ------------------------------------------------------------------ +# Grid definitions +# ------------------------------------------------------------------ +# Oversampling grid: 10 arc-sec +res10 = 10.0 / 3600.0 +lon10 = np.arange(-180.0 + res10 / 2.0, 180.0, res10) # 129600 +lat10 = np.arange(-90.0 + res10 / 2.0, 90.0, res10) # 64800 + +# Target grid: 30 arc-sec = 3 x 10 arc-sec +lon30 = lon10[1::3] # 43200 +# lat30 varies by band; global lat30 would be lat10[1::3] = 21600 + +# Process in latitude bands to limit memory +band_deg = 5.0 + +# Preserve thin line features when rasterizing +ALL_TOUCHED = True + +# ------------------------------------------------------------------ +# Shapefile list +# ------------------------------------------------------------------ +shp_list = sorted(glob.glob(os.path.join(shapefile_dir, "Reaches_pfaf_*.shp"))) +if not shp_list: + raise RuntimeError(f"No shapefiles found: {shapefile_dir}/Reaches_pfaf_*.shp") +print(f"Found {len(shp_list)} reach shapefiles") + + +# ------------------------------------------------------------------ +# Helper: 10" -> 30" aggregation (3x3 blocks) +# ------------------------------------------------------------------ +def agg10_to_30(binary10: np.ndarray): + """ + binary10: uint8 array (nlat10, nlon10) with values 0/1 + Returns: + frac30: float32 (nlat30, nlon30) = mean of 3x3 block + interpreted here as subpixel occupancy by reach lines + any30 : uint8 (nlat30, nlon30) = max of 3x3 block + """ + nlat, nlon = binary10.shape + nlat3 = (nlat // 3) * 3 + nlon3 = (nlon // 3) * 3 + binary10 = binary10[:nlat3, :nlon3] + + blk = binary10.reshape(nlat3 // 3, 3, nlon3 // 3, 3) + frac30 = blk.mean(axis=(1, 3)).astype(np.float32) + any30 = blk.max(axis=(1, 3)).astype(np.uint8) + return frac30, any30 + + +# ------------------------------------------------------------------ +# Build 30" band files +# ------------------------------------------------------------------ +band_files = [] + +for lat_start in np.arange(-90.0, 90.0, band_deg): + lat_end = lat_start + band_deg + print(f"\n=== Latitude band {lat_start} to {lat_end} ===") + + # 10" lat centers in this band + lat10_band = lat10[(lat10 >= lat_start) & (lat10 < lat_end)] + if lat10_band.size == 0: + continue + + out_shape10 = (lat10_band.size, lon10.size) + + # Affine transform for this band (pixel -> lon/lat) + # Keep consistent with the existing lake preprocessor. + transform10 = ( + Affine.translation(lon10[0] - res10 / 2.0, lat10_band[0] - res10 / 2.0) + * Affine.scale(res10, res10) + ) + + # In-memory 10" band raster (binary 0/1) + reach10 = np.zeros(out_shape10, dtype=np.uint8) + + # Rasterize all shapefiles, but only features intersecting this latitude band + for shp in shp_list: + gdf = gpd.read_file(shp).to_crs("EPSG:4326") + gdf_band = gdf.cx[:, lat_start:lat_end] + if gdf_band.empty: + continue + + # Drop null / empty geometries + gdf_band = gdf_band[gdf_band.geometry.notnull()] + gdf_band = gdf_band[~gdf_band.geometry.is_empty] + if gdf_band.empty: + continue + + shapes = ((geom, 1) for geom in gdf_band.geometry) + + tmp = rasterio.features.rasterize( + shapes, + out_shape=out_shape10, + transform=transform10, + fill=0, + dtype=np.uint8, + all_touched=ALL_TOUCHED, + ) + + # Union across shapefiles + reach10 = np.maximum(reach10, tmp) + + # Aggregate to 30" + frac30, any30 = agg10_to_30(reach10) + + # 30" lat centers for this band (aligned with 10" centers) + lat30_band = lat10_band[1::3] + + if frac30.shape != (lat30_band.size, lon30.size): + raise RuntimeError( + f"Shape mismatch: frac30={frac30.shape} vs coords={(lat30_band.size, lon30.size)}" + ) + + ds = xr.Dataset( + { + "reach_occupancy_frac": (("lat", "lon"), frac30), + "reach_presence_any": (("lat", "lon"), any30), + }, + coords={"lat": lat30_band, "lon": lon30}, + ) + + ds.attrs["Source"] = "HydroLAKES-TopoCat v1.1 (2023) Reaches" + ds.attrs["GeometryType"] = "line" + ds.attrs["Rasterization"] = "10 arc-sec oversampling, aggregated to 30 arc-sec" + ds.attrs["AllTouched"] = int(ALL_TOUCHED) + ds.attrs["reach_occupancy_frac_note"] = ( + "Fraction of 10 arc-sec subpixels touched by reach lines; not true river area fraction." + ) + + band_out = os.path.join(out_dir, f"ReachTopoCat_30arcsec_{lat_start}_{lat_end}.nc4") + encoding = { + "reach_occupancy_frac": { + "zlib": True, + "complevel": 4, + "dtype": "float32", + "_FillValue": -9999.0, + }, + "reach_presence_any": { + "zlib": True, + "complevel": 4, + "dtype": "uint8", + "_FillValue": 255, + }, + } + + ds.to_netcdf(band_out, format="NETCDF4", engine="netcdf4", encoding=encoding) + print("Wrote band:", band_out) + band_files.append(band_out) + +# ------------------------------------------------------------------ +# Combine band files into the single global 30" product +# ------------------------------------------------------------------ +print(f"\nCombining {len(band_files)} 30\" band files into:\n {out_global}") + +if not band_files: + raise RuntimeError("No band files were created. Check input shapefiles and paths.") + +ds = xr.open_mfdataset( + sorted(band_files), + combine="by_coords", + parallel=False, +).sortby("lat") + +# Force arrays into memory to avoid dask/lazy-write issues +frac = ds["reach_occupancy_frac"].load().astype("float32") +anyv = ds["reach_presence_any"].load().astype("uint8") + +ds_out = xr.Dataset( + { + "reach_occupancy_frac": frac, + "reach_presence_any": anyv, + }, + coords={ + "lat": ds["lat"].values, + "lon": ds["lon"].values, + }, +) + +ds_out.attrs["Title"] = "HydroLAKES-TopoCat reaches rasterized to 30 arc-sec" +ds_out.attrs["Source"] = "HydroLAKES-TopoCat v1.1 (2023) Reaches" +ds_out.attrs["GeometryType"] = "line" +ds_out.attrs["AllTouched"] = int(ALL_TOUCHED) +ds_out.attrs["reach_occupancy_frac_note"] = ( + "Fraction of 10 arc-sec subpixels touched by reach lines; not true river area fraction." +) + +encoding = { + "reach_occupancy_frac": { + "zlib": True, + "complevel": 4, + "dtype": "float32", + "_FillValue": -9999.0, + }, + "reach_presence_any": { + "zlib": True, + "complevel": 4, + "dtype": "uint8", + "_FillValue": 255, + }, +} + +ds_out.to_netcdf( + out_global, + format="NETCDF4", + engine="netcdf4", + encoding=encoding, +) + +print("Done. Global file written:") +print(" ", out_global) + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/split_lake_30arcsec_to_HV_generic.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/split_lake_30arcsec_to_HV_generic.py new file mode 100755 index 000000000..e8735b466 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/lake/v1/split_lake_30arcsec_to_HV_generic.py @@ -0,0 +1,199 @@ +#!/usr/bin/env python3 +import os +import numpy as np +import xarray as xr +from datetime import datetime + +""" +Split a global 30 arc-sec mask into GEOS 10°x10° H/V tiles. + +This is a lightly generalized version of the HydroLAKES-TopoCat lake splitter so +that the same script can be used for lakes or reaches by changing only the +configuration block below. + +Default configuration is set for lakes. + +Expected global input: + - Dimensions: lat, lon + - Size: 21600 x 43200 (30 arc-sec global) + - Variables: + * one float field in [0,1] (fraction-like field) + * one uint8 field in {0,1} (binary presence field) + +Per-tile output: + - Dimensions: N_lat, N_lon (1200 x 1200) + - Coordinates: latitude(N_lat), longitude(N_lon) + - Variables: copied using names given in CONFIG + - Global attrs: + N_lon_global, N_lat_global + i_ind_offset_LL, j_ind_offset_LL + CellSize_arc_Secs + +Notes: + - No missing values are written anywhere; NaNs are converted to zero. + - Ocean / non-feature cells remain zero. + - For reaches, the float field is typically an occupancy-like fraction, not + a true area fraction. +""" + +# ----------------------------------------------------------------------------- +# CONFIGURATION +# Change only this block to switch between lakes and reaches. +# ----------------------------------------------------------------------------- +### LAKES #### +# ----------------------------------------------------------------------------- +CONFIG = { + # Input/output locations + "in_global": "/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/lake/lake_mask/v1/before_splitting_nc4/LakeTopoCat_Global_30arcsec.nc4", + "out_dir": "/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/lake/lake_mask/v1/", + + # Output filename prefix + "file_prefix": "LakeTopoCat_30arcsec", + + # Variable names in the global input file + "frac_var_in": "lake_area_frac", + "any_var_in": "lake_presence_any", + + # Variable names to write in each tile file + "frac_var_out": "lake_area_frac", + "any_var_out": "lake_presence_any", + + # Metadata strings + "region_label": "LakeTopoCat", + "source_attr": 'HydroLAKES-TopoCat v1.1 (2023); rasterized at 10" and aggregated to 30".', + "created_by": "borescan", +} + +# ----------------------------------------------------------------------------- +### REACHES #### +# ----------------------------------------------------------------------------- +#CONFIG = { + # Input/output locations +# "in_global": "/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/lake/reach_mask/v1/before_splitting_nc4/ReachTopoCat_Global_30arcsec.nc4", +# "out_dir": "/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/lake/reach_mask/v1/", + + # Output filename prefix +# "file_prefix": "ReachTopoCat_30arcsec", + + # Variable names in the global input file +# "frac_var_in": "reach_occupancy_frac", +# "any_var_in": "reach_presence_any", + + # Variable names to write in each tile file +# "frac_var_out": "reach_occupancy_frac", +# "any_var_out": "reach_presence_any", + + # Metadata strings +# "region_label": "ReachTopoCat", +# "source_attr": 'HydroLAKES-TopoCat v1.1 (2023) reaches; rasterized at 10" and aggregated to 30".', +# "created_by": "borescan", +# } + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + +# 30 arc-sec global grid +N_LON_GLOBAL = 43200 +N_LAT_GLOBAL = 21600 + +# 10° x 10° tile size at 30 arc-sec +NC_10 = 1200 +NR_10 = 1200 + + +def main(): + in_global = CONFIG["in_global"] + out_dir = CONFIG["out_dir"] + file_prefix = CONFIG["file_prefix"] + frac_var_in = CONFIG["frac_var_in"] + any_var_in = CONFIG["any_var_in"] + frac_var_out = CONFIG["frac_var_out"] + any_var_out = CONFIG["any_var_out"] + region_label = CONFIG["region_label"] + source_attr = CONFIG["source_attr"] + created_by = CONFIG["created_by"] + + os.makedirs(out_dir, exist_ok=True) + + ds = xr.open_dataset(in_global) + + # Get global 1-D coords + lat_g = ds["lat"].astype("float32") + lon_g = ds["lon"].astype("float32") + + # Data vars + frac_g = ds[frac_var_in].astype("float32") + any_g = ds[any_var_in].astype("uint8") + + # Basic sanity checks + assert frac_g.sizes["lon"] == N_LON_GLOBAL and frac_g.sizes["lat"] == N_LAT_GLOBAL, \ + f"Unexpected grid for {frac_var_in}: {dict(frac_g.sizes)}" + assert any_g.sizes["lon"] == N_LON_GLOBAL and any_g.sizes["lat"] == N_LAT_GLOBAL, \ + f"Unexpected grid for {any_var_in}: {dict(any_g.sizes)}" + + created = datetime.now().strftime("%Y-%m-%d %H:%M:%S") + + for jx in range(1, 18 + 1): # V 01..18 + for ix in range(1, 36 + 1): # H 01..36 + i0 = (ix - 1) * NC_10 # 0-based lon start + j0 = (jx - 1) * NR_10 # 0-based lat start + + hh = f"{ix:02d}" + vv = f"{jx:02d}" + out = os.path.join(out_dir, f"{file_prefix}_H{hh}V{vv}.nc4") + + # Subset coords (1-D) + lat_sub = lat_g.isel(lat=slice(j0, j0 + NR_10)).values + lon_sub = lon_g.isel(lon=slice(i0, i0 + NC_10)).values + + # Subset data (2-D), enforce no-missing + valid ranges + frac_sub = frac_g.isel(lat=slice(j0, j0 + NR_10), lon=slice(i0, i0 + NC_10)).values + any_sub = any_g.isel(lat=slice(j0, j0 + NR_10), lon=slice(i0, i0 + NC_10)).values + + frac_sub = np.nan_to_num(frac_sub, nan=0.0).astype("float32") + frac_sub = np.clip(frac_sub, 0.0, 1.0) + + # ensure binary presence is exactly 0/1 uint8 + any_sub = np.nan_to_num(any_sub, nan=0).astype("uint8") + any_sub = np.clip(any_sub, 0, 1).astype("uint8") + + # Build soil-style dataset + sub = xr.Dataset( + data_vars={ + frac_var_out: (("N_lat", "N_lon"), frac_sub), + any_var_out: (("N_lat", "N_lon"), any_sub), + }, + coords={ + "latitude": ("N_lat", lat_sub), + "longitude": ("N_lon", lon_sub), + }, + ) + + # Global attributes (soil-style) + sub.attrs["N_lon_global"] = N_LON_GLOBAL + sub.attrs["N_lat_global"] = N_LAT_GLOBAL + sub.attrs["i_ind_offset_LL"] = i0 + 1 # 1-based like other products + sub.attrs["j_ind_offset_LL"] = j0 + 1 + sub.attrs["CellSize_arc_Secs"] = 30.0 + sub.attrs["Region"] = f"{region_label} tile H{hh}V{vv}" + sub.attrs["CreatedBy"] = created_by + sub.attrs["Date"] = created + sub.attrs["Source"] = source_attr + + # No fill values (data are complete everywhere after nan_to_num) + encoding = { + frac_var_out: {"zlib": True, "complevel": 4, "dtype": "float32", "_FillValue": None}, + any_var_out: {"zlib": True, "complevel": 4, "dtype": "uint8", "_FillValue": None}, + "latitude": {"zlib": True, "complevel": 1, "dtype": "float32", "_FillValue": None}, + "longitude": {"zlib": True, "complevel": 1, "dtype": "float32", "_FillValue": None}, + } + + sub.to_netcdf(out, format="NETCDF4", engine="netcdf4", encoding=encoding) + print("Wrote", out) + + ds.close() + + +if __name__ == "__main__": + main() +