diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_cube.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_cube.py index 78506693c..50287a2a3 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_cube.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_cube.py @@ -53,7 +53,7 @@ /bin/mv rst/Pfafstetter-M.rst rst/Pfafstetter.rst bin/CombineRasters.x -f 0 -t {NT} {OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter >/dev/null bin/CombineRasters.x -t {NT} CF{NC}x6C{SGNAME} {OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter - bin/mk_runofftbl.x CF{NC}x6C{SGNAME}_{OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + bin/mk_runofftbl.x -g CF{NC}x6C{SGNAME}_{OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} setenv OMP_NUM_THREADS 1 if ({SKIPLAND} != True) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C{SGNAME}_{OCEAN_VERSION}{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_latlon.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_latlon.py index e6e7844e6..d552f7ae4 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_latlon.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_latlon.py @@ -41,7 +41,7 @@ /bin/mv rst/Pfafstetter-M.rst rst/Pfafstetter.rst bin/CombineRasters.x -f 0 -t {NT} {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter >/dev/null bin/CombineRasters.x -t {NT} DC{IM}xPC{JM} {DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter - bin/mk_runofftbl.x DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + bin/mk_runofftbl.x -g DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} setenv OMP_NUM_THREADS 1 if ( {SKIPLAND} != True ) bin/mkCatchParam.x -x {NX} -y {NY} -g DE{IMO}xPE{JMO}_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} setenv OMP_NUM_THREADS {NCPUS} diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 index 298895f5c..88ed9210c 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 @@ -1,235 +1,237 @@ -program Runoff +program mk_runofftbl +! !INTERFACE: +! +! !ARGUMENTS: +! +! Usage = "mk_runofftbl.x -g Gridname -v LBCSV" +! +! -g: Gridname: a string that describes the grids associated with the atmosphere and ocean model configuration +! eg, CF0180x6C_M6TP0072x0036-Pfafstetter +! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09, v10, v11, v12, ...) +! +! This program generates the runoff table *.trn and *.TRN files that are used in the Catchment model for +! directing runoff to its ocean sink. The inputs are (i) bcs geometry files associated with the Gridname +! and (ii) a binary file ("Outlet_latlon.43200x21600") that provides the land raster grid cells where the +! outlets are located. The latter file is either created by [..]/Raster/preproc/routing/run_routing_raster.py +! or from Randy's (Randal.d.koster@nasa.gov) old file under {MAKE_BCS_INPUT_DIR}/land/route/v1. The version +! of the outlet file used by mk_runofftbl.x is determined via the LBCSV argument. +! The program first moves the outlet locations from the land raster grid cells to the nearest ocean pixels +! by calling outlets_to_ocean() and then generates the runoff table files. +! The subroutine outlets_to_ocean() currently works only for the MOM5 and MOM6 tripolar ocean grids. +! +!========================================================= +! +! Yujin Zeng - June 17, 2024 +! Email: yujin.zeng@nasa.gov + use mapl_hashmod use mapl_sortmod + use rmTinyCatchParaMod, only : init_bcs_config, OUTLETV use netcdf - + implicit none + include 'netcdf.inc' integer :: nx, ny, pf integer, allocatable :: lats(:,:), lons(:,:) integer, pointer :: rst(:,:), SortArr(:,:), key(:) integer, pointer :: srctile(:), srcweight(:), dstweight(:), dsttile(:) - real, allocatable :: SrcFraction(:), area(:), in(:), out(:) + real, allocatable :: SrcFraction(:), area(:), in(:), out(:) integer :: i,j,k,l, Hash, HashC, ii,jj,kk integer :: type, np,lnd, is,ie,ww integer :: numtrans, numclosed integer :: status - character*100 :: file, fileT, fileR, fileO, fileB, fileBB + character*100 :: Gridname, fileT, fileR, fileO, fileB character*400 :: fileLL character*400 :: MAKE_BCS_INPUT_DIR - character*5 :: C_NX, C_NY - logical :: adjust_oceanLandSea_mask = .false. ! default is .false. integer :: nxt, command_argument_count - character*(128) :: arg, & - Usage = "mk_runofftbl.x CF0012x6C_TM0072xTM0036-Pfafstetter", & - mapl_tp_file - - + character*(128) :: arg + character*(128) :: Usage = "mk_runofftbl.x -g CF0180x6C_M6TP0072x0036-Pfafstetter -v v12" + character*5 :: LBCSV = 'UNDEF' + character*1 :: opt + + ! ------------------------------------------------------------------ + call get_environment_variable ("MAKE_BCS_INPUT_DIR",MAKE_BCS_INPUT_DIR) - fileLL=trim(MAKE_BCS_INPUT_DIR)//'/land/route/Outlet_latlon.' - + + ! Read inputs ----------------------------------------------------- -! Read inputs ----------------------------------------------------- I = command_argument_count() - if (I < 1 .or. I > 3) then - print *, " " - print *, "Wrong number of input arguments, got: ", I - print *, "Example usage with defaults: " - print *, " " - print *, trim(Usage) - print *, " " - call exit(1) + if (I /= 4) then + print *, " " + print *, "Wrong number of input arguments, got: ", I + print *, "Example usage with defaults: " + print *, " " + print *, trim(Usage) + print *, " " + call exit(1) end if nxt = 1 - call get_command_argument(nxt, file) - print *, " " - print*, "Working with input BCs string: ", file - print *, " " - if (I > 1) then - nxt = nxt + 1 - call get_command_argument(nxt, arg) - if ( trim(arg) .ne. 'yes') then - print *, "Incorrect optional second argument, should be: yes" - call exit(2) - else - adjust_oceanLandSea_mask = .true. - nxt = nxt + 1 - call get_command_argument(nxt, mapl_tp_file) - endif - endif -! ------------------------------------------------------------------ - - fileT = "til/"//trim(file)//".til" ! input - fileR = "rst/"//trim(file)//".rst" ! input - fileO = "til/"//trim(file)//".trn" ! output - fileB = "til/"//trim(file)//".TRN" ! output + call get_command_argument(nxt,arg) + do while(arg(1:1)=='-') + opt=arg(2:2) + if(len(trim(arg))==2) then + nxt = nxt + 1 + call get_command_argument(nxt,arg) + else + arg = arg(3:) + end if + select case (opt) + case ('g') + Gridname = trim(arg) + case ('v') + LBCSV = trim(arg) + call init_bcs_config (trim(LBCSV)) ! get bcs details from version string + case default + print *, "Wrong flag -", opt + print *, "Example usage with defaults: " + print *, trim(Usage) + call exit(1) + end select + nxt = nxt + 1 + call get_command_argument(nxt,arg) + end do -! Read I and J indeces of river outlets. -! These should all be ocean pixels -!--------------------------------------- + if (trim(OUTLETV)=="v1" .or. trim(OUTLETV)=="v2") then + fileLL=trim(MAKE_BCS_INPUT_DIR)//'/land/route/'//trim(OUTLETV)//'/Outlet_latlon.' + else + print *, "Routing files will not be produced with the selected land BCs version (too old)" + stop + endif -! print *, "Getting raster size from "//trim(fileT) + ! ------------------------------------------------------------------ + + fileT = "til/"//trim(Gridname)//".til" ! input + fileR = "rst/"//trim(Gridname)//".rst" ! input + fileO = "til/"//trim(Gridname)//".trn" ! output + fileB = "til/"//trim(Gridname)//".TRN" ! output + + ! Read I and J indices of river outlets. + ! These should all be ocean pixels + ! --------------------------------------- + + ! print *, "Getting raster size from "//trim(fileT) open(10,file=fileT, form="formatted", status="old") - + read(10,*) np, pf, nx, ny close(10) -! print *, nx, ny - + write (C_NX, '(i5.5)') NX write (C_NY, '(i5.5)') NY - + print *, "Reading outlets..." - + allocate(lats(nx,ny), lons(nx,ny),stat=status) if(status/=0) then print *, "Out of Memory" stop __LINE__ endif - + open (30,file=trim(fileLL)//C_NX//'x'//C_NY,form="unformatted",status="old") do j = 1, ny read (30) lons(:,j) read (30) lats(:,j) end do close(30) + + print *, " " + print *, "Determining river outlets to ocean:" + print *, "- Output file: ", fileB + print *, " " -! do j=1,ny -!! if (mod(j,100) == 0) print *,'J=',j -! do i=1,nx -! ii = Lons(i,j) -! jj = lats(i,j) -! -! if(ii==-999 .or. jj==-999) then -! ! ii = i -! ! jj = j -! cycle -! endif -! -! if(ii==i .and. jj==j) then -! print *, '>>> Inland Ocean Point ', ii, jj, lons(i,j), lats(i,j) -! stop -! end if -! -! end do -! end do -! stop "DONE" -! Count the number of Ocean and land tiles in the tile file -! All land tiles preceed the ocean tiles. -!---------------------------------------------------------- - -! If asked for, adjust tiles to be -! comptabile with ocean model land-sea mask and write ANOTHER output file -!------------------------------------------------------------------------- - - if (adjust_oceanLandSea_mask) then - fileBB = "til/"//trim(file)//"_oceanMask_adj.TRN" ! output - - print *, " " - print *, "Accounting for any mismatch between land-sea masks:" - print *, "- Of GEOS land and external ocean model." - print *, "- Output file: ", fileB - print *, " " - call read_oceanModel_mask( mapl_tp_file) -! ... some adjustment of following variable: `type` -! ... using ocean model land-sea mask should be done here - endif - -! print *, "Reading til file "//trim(fileT) - + call outlets_to_ocean(Gridname,lons,lats,nx,ny) + open(10,file=fileT, form="formatted", status="old") - + read(10,*) np allocate(area(np), in(np), out(np)) - + out = 0.0 in = 0.0 - + do i=1,7 read(10,*) enddo - + lnd=-1 do l=1,np read(10,*) type, area(l) if(type==0 .and. lnd==-1) lnd = l-1 if(lnd<0) in(l) = 1. end do - + print *, "area of sphere = ", sum(real(area,kind=8)) print *, "area of land = ", sum(real(area*in,kind=8)) - + close(10) print *, "Number of Tiles ", np print *, "Ocean Tiles ", np-lnd print *, "Land tiles ", lnd - -! Read the raster file -!--------------------- - - print *, "Reading rst file "//trim(fileR) - + + ! Read the raster file + ! --------------------- + + print *, "Reading raster (rst) file "//trim(fileR) + open(20,file=fileR,form="unformatted",status="old") - + allocate(rst(nx,ny),stat=status) if(status/=0) then print *, "Out of Memory" stop __LINE__ endif - + do j=1,ny read(20) rst(:,j) enddo - + close(20) - + allocate(SortArr(2*lnd,3)) DstTile => SortArr(:,1) SrcTile => SortArr(:,2) SrcWeight => SortArr(:,3) - -! Hash the raster -!---------------- - + + ! Hash the raster + ! ---------------- + print *, "Hashing raster... " - + Hash = MAPL_HashCreate(1024) HashC = MAPL_HashCreate(1024) - + NumTrans=0 - + do j=1,ny -! if (mod(j,100) == 0) print *,'J=',j do i=1,nx if(rst(i,j)<=lnd) then ii = Lons(i,j) jj = lats(i,j) - + if(ii==-999 .or. jj==-999) then -! ii = i -! jj = j cycle endif - + if(ii==i .and. jj==j) then print *, '>>> Inland Ocean Point ', ii, jj, rst(i,j) stop end if - + k = MAPL_HASHIncrement(HashC,rst(i,j)) k = MAPL_HASHIncrement(Hash,rst(ii,jj),rst(i,j)) - + SrcWeight(k) = SrcWeight(k) + 1 - + if(k>NumTrans) then if(k/=NumTrans+1) then print *, NumTrans, k @@ -239,52 +241,52 @@ program Runoff SrcTile(NumTrans) = rst(i ,j ) DstTile(NumTrans) = rst(ii,jj) endif - + end if end do end do - + DstTile => SortArr(:NumTrans,1) SrcTile => SortArr(:NumTrans,2) SrcWeight => SortArr(:NumTrans,3) - + print *, "Total Transactions ", NumTrans print *, MAPL_HashSize(Hash),MAPL_HashSize(HashC) call MAPL_HashDestroy(Hash) call MAPL_HashDestroy(HashC) - -! Allocate space for transanction lists -!-------------------------------------- - + + ! Allocate space for transanction lists + ! -------------------------------------- + allocate(key(numTrans),stat=status) - + if(status/=0) then print *, "Out of Memory" stop endif - -! Sort transactions by source tile number to compute source -! fractions going into each transaction. -!---------------------------------------------------------- - + + ! Sort transactions by source tile number to compute source + ! fractions going into each transaction. + ! ---------------------------------------------------------- + print *, "Sorting transactions by source..." - + Key = SrcTile call MAPL_Sort(Key,SortArr(:NumTrans,:),DIM=1) - + print *, "Computing weights..." - + deallocate(key) allocate(SrcFraction(numTrans)) - -! Compute fractions -!------------------ - + + ! Compute fractions + ! ------------------ + is = 1 ie = 1 - + do j=2,NumTrans if(SrcTile(j)/=SrcTile(is)) then SrcFraction(is:ie) = SrcWeight (is:ie) @@ -295,22 +297,24 @@ program Runoff ie = ie + 1 end if end do - + SrcFraction(is:ie) = SrcWeight (is:ie) SrcFraction(is:ie) = SrcFraction(is:ie) / float(sum(SrcWeight(is:ie))) - + print *,"SrcWeight", sum(SrcFraction), lnd - + print *, '<<<', sum(SrcFraction*area(SrcTile)) - + SrcFraction = SrcFraction * (area(SrcTile)/area(DstTile)) - + print *, '>>>', sum(SrcFraction*area(DstTile)) - -! Write output files -!------------------- - - print *, "Writing output file..." + + ! Write output files + ! ------------------- + + print *, "Writing output files..." + + ! ASCII file (*.trn) open(10,file=fileO, form="formatted", status="unknown") write(10,*) NumTrans @@ -319,81 +323,705 @@ program Runoff end do close(10) + ! binary file (*.TRN) + call write_route_file( fileB, NumTrans, SrcTile, DstTile, SrcFraction) - if (adjust_oceanLandSea_mask) & - call write_route_file( fileBB, NumTrans, SrcTile, DstTile, SrcFraction) - + do j=1,NumTrans Out(DstTile(j)) = Out(DstTile(j)) + In(SrcTile(J))*SrcFraction(J) enddo print *, "area of land = ", sum(real(area*out,kind=8)) - - + print *, "Completed successfully" - + deallocate( SrcFraction) deallocate( SortArr) deallocate( rst) deallocate( area) deallocate( lats, lons) - + call exit(0) + + ! ----------------------------------------------------------------- -! ----------------------------------------------------------------- contains - - subroutine read_oceanModel_mask( mask_file) - implicit none - character*128, intent(in) :: mask_file - integer :: nx, ny - real, allocatable :: wetMask(:,:) - - integer :: ncid, varid - - print *, "Reading ocean model mask from : ", mask_file - call check( nf90_open(mask_file, nf90_nowrite, ncid)) ! open nc file - - call check( nf90_inq_dimid(ncid, "n_center_x", varid)) ! read dimenstion (x) - call check( nf90_inquire_dimension(ncid, varid, len=nx)) - - call check( nf90_inq_dimid(ncid, "n_center_y", varid)) ! read dimenstion (y) - call check( nf90_inquire_dimension(ncid, varid, len=ny)) - - allocate( wetMask(nx, ny)) - call check( nf90_inq_varid(ncid, "mask", varid)) ! read mask - call check( nf90_get_var(ncid, varid, wetMask)) - - call check( nf90_close(ncid)) ! close nc file - - deallocate( wetMask) - end subroutine read_oceanModel_mask -! ---------------------- - + subroutine write_route_file( fileB, NumTrans, SrcTile, DstTile, SrcFraction) - implicit none - character*100, intent(in) :: fileB - integer, intent(in) :: NumTrans - integer, pointer, intent(in) :: srctile(:), dsttile(:) - real, intent(in) :: SrcFraction(:) - - open(10,file=fileB, form="unformatted", status="unknown") - write(10) NumTrans - write(10) SrcTile - write(10) DstTile - write(10) SrcFraction - close(10) + + character*100, intent(in) :: fileB + integer, intent(in) :: NumTrans + integer, pointer, intent(in) :: srctile(:), dsttile(:) + real, intent(in) :: SrcFraction(:) + + open(10,file=fileB, form="unformatted", status="unknown") + write(10) NumTrans + write(10) SrcTile + write(10) DstTile + write(10) SrcFraction + close(10) + end subroutine write_route_file -! ---------------------- - - subroutine check(status) - implicit none - integer, intent (in) :: status - if (status /= nf90_noerr) then - print *, trim(nf90_strerror(status)) - print *, "Error in reading ocean mask file." - stop - end if - end subroutine check -! ----------------------------------------------------------------- - -end program Runoff + + ! ------------------------------------------------------------------------ + ! The subroutine moves outlets of each land & Greenland gridcell (with endpoint to ocean only) + ! to the nearest ocean gridcell defiend by "Gridname" and ocean mask if the original outlets are not in the ocean. + + subroutine outlets_to_ocean(Gridname,lons,lats,nx,ny) + + integer, intent(in) :: nx,ny !number of lon and lat, eg: 43200 and 21600 in 30s resolution + character(len=*), intent(in) :: Gridname !name of the domain, eg: CF0180x6C_M6TP0072x0036-Pfafstetter + integer, intent(inout) :: lons(nx,ny),lats(nx,ny) + !lons(nx,ny): lon idx of outlets of each land&Greenland gridcell (with endpoint to ocean only) + !lats(nx,ny): lat idx of outlets of each land&Greenland gridcell (with endpoint to ocean only) + !lons(nx,ny) and lats(nx,ny) are got from input binary file, eg: Outlet_latlon.43200x21600 + ! ----------------------------------------------------------- + + integer, allocatable, dimension(:) :: lati_lnd,loni_lnd + integer, allocatable, dimension(:) :: msk1d + integer, allocatable, dimension(:,:) :: msk2d + integer, allocatable, dimension(:,:) :: mask + integer, allocatable, dimension(:,:) :: boundary + real*8, allocatable, dimension(:) :: lonsh,latsh + real*8, allocatable, dimension(:) :: lons_adj,lats_adj + integer, allocatable, dimension(:) :: lati_ocn,loni_ocn + character*100 :: Gridname_ocn + character*100 :: Gridname_ocn_lnd + character*100 :: res_MAPL + integer, allocatable, dimension(:,:) :: rst_ocn,rst_ocn_lnd + integer :: nt_ocn_lnd,nl_ocn_lnd,nt_ocn,nx_MAPL,ny_MAPL,nsh + integer, pointer, dimension(:) :: t2lati,t2loni + real*8, allocatable, dimension(:) :: lon30s,lat30s + integer :: ns + integer, allocatable, dimension(:,:) :: ns_map + real*8, allocatable, dimension(:) :: lat_lnd,lon_lnd + integer :: i,j,l,k,status,type,np,flag,flag2 + + call get_domain_name(Gridname,Gridname_ocn,Gridname_ocn_lnd,res_MAPL,nx_MAPL,ny_MAPL) + allocate(rst_ocn(nx,ny),rst_ocn_lnd(nx,ny)) + allocate(lon30s(nx),lat30s(ny)) + call read_rst_til_files(nx,ny,Gridname_ocn,Gridname_ocn_lnd,rst_ocn,rst_ocn_lnd,t2loni,t2lati,nt_ocn_lnd,nl_ocn_lnd,nt_ocn,lon30s,lat30s) + !print *,"running outlets_num() ..." + call outlets_num(rst_ocn_lnd,nl_ocn_lnd,nt_ocn_lnd,lons,lats,nx,ny,ns) + !print *,"outlets num is ",ns + allocate(loni_lnd(ns),lati_lnd(ns)) + allocate(lons_adj(ns),lats_adj(ns)) + allocate(loni_ocn(ns),lati_ocn(ns)) + allocate(ns_map(nx,ny)) + allocate(lon_lnd(ns),lat_lnd(ns)) + !print *,"running retrieve_outlets() ..." + call retrieve_outlets(lons,lats,lon30s,lat30s,loni_lnd,lati_lnd,lon_lnd,lat_lnd,ns_map,nx,ny,ns) + !print *,"running mask_MAPL_1d() ..." + allocate(msk1d(nt_ocn)) + call mask_MAPL_1d(msk1d,t2loni,t2lati,nt_ocn,res_MAPL,nx_MAPL,ny_MAPL) + deallocate(t2loni,t2lati) + !print *,"running mask_MAPL_2d() ..." + allocate(msk2d(nx,ny)) + call mask_MAPL_2d(rst_ocn,msk1d,msk2d,nt_ocn,nx,ny) + deallocate(rst_ocn,msk1d) + !print *,"running mask_MAPL_bcs() ..." + allocate(mask(nx,ny)) + call mask_MAPL_bcs(rst_ocn_lnd,msk2d,mask,nx,ny,nl_ocn_lnd,nt_ocn_lnd) + deallocate(msk2d,rst_ocn_lnd) + !print *,"running ocean_boundary() ..." + allocate(boundary(nx,ny)) + call ocean_boundary(mask,boundary,nx,ny) + !print *,"running ocean_boundary_num() ..." + call ocean_boundary_num(boundary,nx,ny,nsh) + !print *,"ocean boundary point num is ",nsh + allocate(lonsh(nsh),latsh(nsh)) + !print *,"running ocean_boundary_points() ..." + call ocean_boundary_points(boundary,lon30s,lat30s,lonsh,latsh,nx,ny,nsh) + deallocate(boundary) + !print *,"running move_to_ocean() ..." + call move_to_ocean(loni_lnd,lati_lnd,lon_lnd,lat_lnd,mask,lonsh,latsh,lons_adj,lats_adj,ns,nx,ny,nsh) + deallocate(mask,lonsh,latsh) + !print *,"running sinkxy_ocean() ..." + call sinkxy_ocean(lons_adj,lats_adj,lon30s,lat30s,loni_ocn,lati_ocn,ns,nx,ny) + !print *,"running update_outlets() ..." + call update_outlets(loni_ocn,lati_ocn,ns_map,lons,lats,nx,ny,ns) + + deallocate(loni_lnd,lati_lnd,lons_adj,lats_adj,loni_ocn,lati_ocn) + deallocate(lon30s,lat30s) + deallocate(ns_map,lon_lnd,lat_lnd) + + end subroutine outlets_to_ocean + + !------------------------------------------------------------------------- + ! This subroutine gets the name of 'Gridname_ocn' and 'Gridname_ocn_lnd' from the input name 'Gridname'. + ! It also gets the resolution of the ocean domain. + subroutine get_domain_name(Gridname,Gridname_ocn,Gridname_ocn_lnd,res_MAPL,nx_MAPL,ny_MAPL) + + character(len=*), intent(in) :: Gridname !input domain name, eg: CF0180x6C_M6TP0072x0036-Pfafstetter + character(len=*), intent(out) :: Gridname_ocn,Gridname_ocn_lnd,res_MAPL + !Gridname_ocn: ocean domain name, eg: M6TP0072x0036 + !Gridname_ocn_lnd: ocean-land domain name, eg: M6TP0072x0036-Pfafstetter + !res_MAPL: ocean resolution name, eg: M6TP0072x0036 + integer, intent(out) :: nx_MAPL,ny_MAPL + !nx_MAPL: number of lon of ocean domain, eg: 72 + !ny_MAPL: number of lat of ocean domain, eg: 36 + + character*100 :: nx_str,ny_str,fileT + integer :: px,plats,plate,plons,plone,plonss,pocns,pocne + integer :: nstr1,nstr2 + integer :: i,length + + Gridname_ocn="" + Gridname_ocn_lnd="" + res_MAPL="" + nx_MAPL=-9999 + ny_MAPL=-9999 + fileT = "til/"//trim(Gridname)//".til" + open(10,file=fileT, form="formatted", status="old") + do i=1,5 + read(10,*) + enddo + read(10,*)Gridname_ocn_lnd + do i=100,1,-1 + if(Gridname_ocn_lnd(i:i).eq."-")then + Gridname_ocn(1:i-1)=Gridname_ocn_lnd(1:i-1) + exit + endif + enddo + read(10,*)nx_MAPL + read(10,*)ny_MAPL + nx_str="" + ny_str="" + write(nx_str,*)nx_MAPL + write(ny_str,*)ny_MAPL + res_MAPL="" + length = len( trim(adjustl(nx_str))//"x"//trim(adjustl(ny_str)) ) + res_MAPL(1:length)=trim(adjustl(nx_str))//"x"//trim(adjustl(ny_str)) + + + end subroutine get_domain_name + + !------------------------------------------------------------------------- + ! This subroutine reads rst and til files + + subroutine read_rst_til_files(nx,ny,Gridname_ocn,Gridname_ocn_lnd,rst_ocn,rst_ocn_lnd,t2loni,t2lati,nt_ocn_lnd,nl_ocn_lnd,nt_ocn,lon30s,lat30s) + + integer, intent(in) :: nx,ny !number of lon and lat, eg: 43200 and 21600 in 30s resolution + character(len=*), intent(in) :: Gridname_ocn,Gridname_ocn_lnd ! input filename, eg: M6TP0072x0036 and M6TP0072x0036-Pfafstetter + + integer, intent(out) :: rst_ocn(nx,ny),rst_ocn_lnd(nx,ny) !data from rst files + integer,pointer,intent(out), dimension(:) :: t2loni,t2lati !relationship between ocn tile idx in M6TP0072x0036.rst and lat/lon idx in MAPL_Tripolar.nc + integer, intent(out) :: nt_ocn_lnd,nl_ocn_lnd,nt_ocn + !nt_ocn_lnd: number of total tiles in the M6TP0072x0036-Pfafstetter + !nl_ocn_lnd: number of land tiles in the M6TP0072x0036-Pfafstetter + !nt_ocn: number of total tiles in the M6TP0072x0036 + real*8, intent(out) :: lon30s(nx),lat30s(ny)!lon and lat value arrays of the 30s map + + character*100 :: fileT_ocn, fileR_ocn + character*100 :: fileT_ocn_lnd, fileR_ocn_lnd + integer :: i,j,type,np,k,l + real :: num1,num2,num3,num4 + real*8 :: dx,dy + + fileT_ocn = "til/"//trim(Gridname_ocn)//".til" ! input + fileR_ocn = "rst/"//trim(Gridname_ocn)//".rst" ! input + fileT_ocn_lnd = "til/"//trim(Gridname_ocn_lnd)//".til" ! input + fileR_ocn_lnd = "rst/"//trim(Gridname_ocn_lnd)//".rst" ! input + + !print *, "Reading rst file "//trim(fileR_ocn) + open(20,file=fileR_ocn,form="unformatted",status="old") + do j=1,ny + read(20) rst_ocn(:,j) + enddo + close(20) + + !print *, "Reading rst file "//trim(fileR_ocn_lnd) + open(21,file=fileR_ocn_lnd,form="unformatted",status="old") + do j=1,ny + read(21) rst_ocn_lnd(:,j) + enddo + close(21) + + open(10,file=fileT_ocn, form="formatted", status="old") + read(10,*) nt_ocn + allocate(t2lati(nt_ocn),t2loni(nt_ocn)) + do i=1,4 + read(10,*) + enddo + do l=1,nt_ocn + read(10,*)type,num1,num2,num3,t2loni(l),t2lati(l) + enddo + close(10) + + open(10,file=fileT_ocn_lnd, form="formatted", status="old") + read(10,*) np + do i=1,4 + read(10,*) + enddo + k=0 + do l=1,np + read(10,*)type,num1,num2,num3,num4 + if(type/=100)exit + k=k+1 + enddo + close(10) + nt_ocn_lnd=np + nl_ocn_lnd=k + + dx=360.d0/nx + dy=180.d0/ny + do i=1,nx + lon30s(i)=-180.d0+dx/2.d0+dx*(i-1) + enddo + do j=1,ny + lat30s(j)=-90.d0+dy/2.d0+dy*(j-1) + enddo + + end subroutine read_rst_til_files + + !------------------------------------------------------------------------- + ! This subroutine counts the number of outlet points from input outlets lon/lat idx map. + + subroutine outlets_num(rst_ocn_lnd,nl,nt,lons,lats,nx,ny,ns) + + integer, intent(in) :: nx,ny,nl,nt + !nx: number of lon of 30s map, 43200 + !ny: number of lat of 30s map, 21600 + !nl: number of land tiles + !nt: number of total tiles + integer, intent(inout) :: lons(nx,ny),lats(nx,ny) !map of lon/lat idx of outlets for each cells on the 30s map + integer, intent(in) :: rst_ocn_lnd(nx,ny) !map of tile idx from ocean-land rst map, eg: from M6TP0072x0036-Pfafstetter.rst + integer, intent(out) :: ns !number of the outlets + + integer, allocatable, dimension(:) :: lonp,latp + integer, allocatable, dimension(:,:) :: acc,np_map + integer :: i,j,k,l,lonc,latc,flag,maxbak,status,num + + allocate(acc(nx,ny)) + + !It first masks out the outlets map with the ocean-land rst map... + !If a cell is not defined as land or glacier in the ocean-land rst map, we ignore it. + do i=1,nx + do j=1,ny + if(rst_ocn_lnd(i,j)>nl.and.rst_ocn_lnd(i,j)/=nt)then + lons(i,j)=-999 + lats(i,j)=-999 + endif + enddo + enddo + + !Counting the outlets... + acc=0 + k=0 + do i=1,nx + do j=1,ny + if(lons(i,j)/=-999.and.lats(i,j)/=-999)then + lonc=lons(i,j) + latc=lats(i,j) + if(acc(lonc,latc)==0)then + k=k+1 + acc(lonc,latc)=1 + endif + endif + enddo + enddo + ns=k + deallocate(acc) + + end subroutine outlets_num + + !------------------------------------------------------------------------ + ! This subroutine retrieves the outlet locations from the outlets map (nx,ny) to a list (ns). + ! It also stores the outlets idx on the 30s map. + + subroutine retrieve_outlets(lons,lats,lon30s,lat30s,lonp,latp,lon_lnd,lat_lnd,ns_map,nx,ny,ns) + + integer, intent(in) :: nx,ny,ns + !nx: number of lon of 30s map, 43200 + !ny: number of lat of 30s map, 21600 + !ns: number of the outlets + integer, intent(in) :: lons(nx,ny),lats(nx,ny) !input map of outlets idx + real*8, intent(in) :: lon30s(nx),lat30s(ny) !lon and lat value arrays of the 30s map + integer, intent(out) :: lonp(ns),latp(ns) !list of lon and lat idx for the ns outlets + real*8, intent(out) :: lon_lnd(ns),lat_lnd(ns) !list of lon and lat value for the ns outlets + integer, intent(out) :: ns_map(nx,ny) !It stores the outlets idx on the 30s map + + integer, allocatable,dimension(:,:) :: acc + integer :: i,j,k,l,lonc,latc + + allocate(acc(nx,ny)) + ns_map=-9999 + acc=0 + k=0 + do i=1,nx + do j=1,ny + if(lons(i,j)/=-999.and.lats(i,j)/=-999)then + lonc=lons(i,j) + latc=lats(i,j) + if(acc(lonc,latc)==0)then + k=k+1 + acc(lonc,latc)=1 + lonp(k)=lonc + latp(k)=latc + ns_map(lonc,latc)=k + endif + endif + enddo + enddo + + do i=1,ns + lon_lnd(i)=lon30s(lonp(i)) + lat_lnd(i)=lat30s(latp(i)) + enddo + + deallocate(acc) + + end subroutine retrieve_outlets + + !------------------------------------------------------------------------ + ! convert the ocean mask in MAPL_Tripolar.nc to 1d list for each ocean tile defined by ocn rst, eg: M6TP0072x0036.rst + + subroutine mask_MAPL_1d(msk_tile,t2loni,t2lati,nt,res_MAPL,nlon,nlat) + + integer, intent(in) :: nt,nlon,nlat + !nt: number of ocean tiles + !nlon: number of lon in MAPL_Tripolar.nc, eg 72 + !nlat: number of lat in MAPL_Tripolar.nc, eg 36 + integer, intent(in) :: t2loni(nt),t2lati(nt) !relationship between ocn tile idx in M6TP0072x0036.rst and lat/lon idx in MAPL_Tripolar.nc + character(len=*), intent(in) :: res_MAPL !name of the ocean resolution, eg M6TP0072x0036 + integer, intent(out) :: msk_tile(nt) !1d list for the ocean msk for each ocean tile + + real, allocatable, dimension(:,:) :: msk_MAPL + integer :: i + + allocate(msk_MAPL(nlon,nlat)) + call read_oceanModel_mapl(res_MAPL,msk_MAPL,nlon,nlat) + + do i=1,nt + msk_tile(i)=int(msk_MAPL(t2loni(i),t2lati(i))) + enddo + + deallocate(msk_MAPL) + + end subroutine mask_MAPL_1d + + !------------------------------------------------------------------------ + ! read ocean mask from "MAPL_Tripolar.nc" + + subroutine read_oceanModel_mapl(res_MAPL,wetMask,nx,ny) + + + character(len=*), intent(in) :: res_MAPL !ocn resolution name + integer, intent(in) :: nx, ny !ocn resolution numbers + real, intent(out) :: wetMask(nx,ny) !ocn mask from MAPL_Tripolar.nc + + integer :: ncid, varid, ret + character(len=4) :: subname="read" + + ! try MOM6 first + + ret=nf90_open( trim(MAKE_BCS_INPUT_DIR) // "/ocean/MOM6/" // trim(res_MAPL) // "/MAPL_Tripolar.nc", 0, ncid ) + + ! if MOM6 did not work, try MOM5, + + if(ret /= NF_NOERR)then + call check_ret( nf90_open( trim(MAKE_BCS_INPUT_DIR) // "/ocean/MOM5/" // trim(res_MAPL) // "/MAPL_Tripolar.nc", 0, ncid ), subname) + endif + + ! read "mask" from netcdf file into "wetMask" + + call check_ret( nf90_inq_varid(ncid,"mask",varid), subname) + call check_ret( nf90_get_var( ncid,varid,wetMask), subname) + call check_ret( nf90_close( ncid), subname) + + end subroutine read_oceanModel_mapl + + !------------------------------------------------------------------------ + + subroutine check_ret(ret, calling) + + integer, intent(in) :: ret + character(len=*) :: calling + + if (ret /= NF_NOERR) then + write(6,*)'netcdf error from ',trim(calling) + call endrun(nf_strerror(ret)) + end if + + end subroutine check_ret + + !----------------------------------------------------------------------- + + subroutine endrun(msg,subname) + + character(len=*), intent(in), optional :: msg ! string to be printed + character(len=*), intent(in), optional :: subname ! subname + + if (present(subname)) then + write(6,*) 'ERROR in subroutine :', trim(subname) + end if + + if (present(msg)) then + write(6,*)'ENDRUN: ', msg + else + write(6,*) 'ENDRUN: called without a message string' + end if + + stop + + end subroutine endrun + + !------------------------------------------------------------------------ + ! convert 1d list of ocn mask to a 30s map + + subroutine mask_MAPL_2d(ocean,mask1d,msk2d,nt,nlon,nlat) + + integer, intent(in) :: nt,nlon,nlat + !nt: number of ocean tiles + !nlon: number of lon in 30s map, eg 43200 + !nlat: number of lat in 30s map, eg 21600 + integer, intent(in) :: ocean(nlon,nlat) !tile idx from ocn rst file, eg: from M6TP0072x0036.rst + integer, intent(in) :: mask1d(nt) !1d list of ocn mask got from subroutine mask_MAPL_1d + integer, intent(out) :: msk2d(nlon,nlat) !output of the ocn mask on the 30s map + + real*8, allocatable,dimension(:) :: lon,lat + integer :: i,j,xi,yi,tid + + do i=1,nlon + do j=1,nlat + tid=ocean(i,j) + msk2d(i,j)=mask1d(tid) + enddo + enddo + + end subroutine mask_MAPL_2d + + !------------------------------------------------------------------------ + ! further mask the ocn mask from MAPL_Tripolar.nc based on the land-ocean rst eg: M6TP0072x0036-Pfafstetter.rst + + subroutine mask_MAPL_bcs(rst_ocn_lnd,mask_mapl,mask,nlon,nlat,nl,nt) + + integer,intent(in) :: nlon,nlat,nl,nt + !nlon: number of lon in 30s map, eg 43200 + !nlat: number of lat in 30s map, eg 21600 + !nl: number of lnd tile + !nt: number of total tile + integer,intent(in) :: rst_ocn_lnd(nlon,nlat) !tile idx from land-ocean rst eg: M6TP0072x0036-Pfafstetter.rst + integer,intent(in) :: mask_mapl(nlon,nlat) !ocn mask map from subroutine mask_MAPL_2d + integer,intent(out) :: mask(nlon,nlat) !ocn mask map masked further by land-ocean rst + + mask=0 + !tile idx<=nl are land tiles + !tile idx==nt are glacier tiles + !tile idx==nt-1 are lake tiles + !rest are all ocean tiles + where(rst_ocn_lnd>nl.and.rst_ocn_lnd/=nt.and.rst_ocn_lnd/=nt-1.and.mask_mapl==1)mask=1 + + end subroutine mask_MAPL_bcs + + !------------------------------------------------------------------------ + ! find the ocean boundary cells (that are next to non-ocean cell) on the 30s map based on the ocn mask from mask_MAPL_bcs + + subroutine ocean_boundary(mask,boundary,nlon,nlat) + + integer, intent(in) :: nlon,nlat + !nlon: number of lon in 30s map, eg 43200 + !nlat: number of lat in 30s map, eg 21600 + integer, intent(in) :: mask(nlon,nlat) !ocn mask from subroutine mask_MAPL_bcs + integer, intent(out) :: boundary(nlon,nlat) !map of ocean boundary cells + + real*8, allocatable :: lon(:),lat(:) + integer :: xi,yi,id + integer :: xp1,xm1,yp1,ym1 + + boundary=mask + boundary=-9999 + + do xi=2,nlon-1 + do yi=2,nlat-1 + id=mask(xi,yi) + if(id==1)then + boundary(xi,yi)=0 + if(mask(xi+1,yi)==1.and.& + mask(xi+1,yi-1)==1.and.& + mask(xi ,yi-1)==1.and.& + mask(xi-1,yi-1)==1.and.& + mask(xi-1,yi)==1.and.& + mask(xi-1,yi+1)==1.and.& + mask(xi ,yi+1)==1.and.& + mask(xi+1,yi+1)==1)then + boundary(xi,yi)=-9999 + endif + endif + enddo + enddo + + end subroutine ocean_boundary + + !------------------------------------------------------------------------ + ! count the number of ocean boundary cells + + subroutine ocean_boundary_num(boundary,nlon,nlat,nsh) + + integer, intent(in) :: nlon,nlat + !nlon: number of lon in 30s map, eg 43200 + !nlat: number of lat in 30s map, eg 21600 + integer, intent(in) :: boundary(nlon,nlat) !map of ocean boundary cells + integer, intent(out) :: nsh !number of ocean boundary cells + + integer :: i,xi,yi,k + + k=0 + do xi=1,nlon + do yi=1,nlat + if(boundary(xi,yi)==0)then + k=k+1 + endif + enddo + enddo + nsh=k + + end subroutine ocean_boundary_num + + !------------------------------------------------------------------------ + ! list the lat and lon of ocean boundary cells + + subroutine ocean_boundary_points(boundary,lon30s,lat30s,lonsh,latsh,nlon,nlat,nsh) + + integer,intent(in) :: nlon,nlat,nsh + !nlon: number of lon in 30s map, eg 43200 + !nlat: number of lat in 30s map, eg 21600 + !nsh: number of ocean boundary cells + integer,intent(in) :: boundary(nlon,nlat) !map of ocean boundary cells + real*8,intent(in) :: lon30s(nlon),lat30s(nlat) !lon and lat value arrays of the 30s map + real*8,intent(out) :: lonsh(nsh),latsh(nsh) !lists of the lat and lon values of the ocean boundary cells + integer i,xi,yi,k + + k=0 + do xi=1,nlon + do yi=1,nlat + if(boundary(xi,yi)==0)then + k=k+1 + lonsh(k)=lon30s(xi) + latsh(k)=lat30s(yi) + endif + enddo + enddo + end subroutine ocean_boundary_points + + !------------------------------------------------------------------------ + ! move the outlet locations to the nearest ocean boundary cell + + subroutine move_to_ocean(lonsi,latsi,lons,lats,mask,lonsh,latsh,lons_adj,lats_adj,ns,nlon,nlat,nsh) + + integer, intent(in) :: ns,nlon,nlat,nsh + !ns: number of the outlets + !nlon: number of lon in 30s map, eg 43200 + !nlat: number of lat in 30s map, eg 21600 + !nsh: number of ocean boundary cells + integer, intent(in) :: lonsi(ns),latsi(ns) !lon and lat idx of the outlets before moving + real*8, intent(in) :: lons(ns),lats(ns) !lon and lat values of the outlets before moving + integer, intent(in) :: mask(nlon,nlat) !!ocn mask from subroutine mask_MAPL_bcs + real*8, intent(in) :: lonsh(nsh),latsh(nsh) !lists of the lat and lon values of the ocean boundary cells + real*8, intent(out) :: lons_adj(ns),lats_adj(ns) !lon and lat values of the outlets after moving + + real,allocatable :: dist(:) + + integer :: i,j + real :: dy,dy2,dx,dx2,dxA,dxB,dist_temp + + allocate(dist(ns)) + do i=1,ns + IF(mask(lonsi(i),latsi(i))==0)THEN + dist(i)=1.e12 + do j=1,nsh + dy=abs(lats(i)-latsh(j)) + dy2=dy*dy + dxA=abs(lons(i)-lonsh(j)) + dxB=360.-dxA + dx=min(dxA,dxB) + dx2=dx*dx + dist_temp=sqrt(dx2+dy2) + if(dist_tempns)then + print *,"ns_map is Incorrect, ind=",ind + stop + endif + lons(i,j)=loni_ocn(ind) + lats(i,j)=lati_ocn(ind) + + endif + enddo + enddo + + end subroutine update_outlets + + !------------------------------------------------------------------------ + +end program mk_runofftbl + +! ============================ EOF ===================================================== 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 ad373f92d..d7867004d 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 @@ -52,6 +52,7 @@ module rmTinyCatchParaMod character*6, public, save :: SOILBCS = 'UNDEF' character*6, public, save :: MODALB = 'UNDEF' character*10, public, save :: SNOWALB = 'UNDEF' + character*5, public, save :: OUTLETV = 'UNDEF' REAL, public, save :: GNU = MAPL_UNDEF character*512 :: MAKE_BCS_INPUT_DIR @@ -91,7 +92,13 @@ SUBROUTINE init_bcs_config (LBCSV) ! NGDC : Soil parameters from Reynolds et al. 2000, doi:10.1029/2000WR900130 (MERRA-2, Fortuna, Ganymed, Icarus) ! HWSD : Merged HWSDv1.21-STATSGO2 soil properties on 43200x21600 with Woesten et al. (1999) parameters ! HWSD_b : As in HWSD but with surgical fix of Argentina peatland issue (38S,60W) - + ! + ! OUTLETV: Definition of outlet locations. DEFAULT : N/A + ! N/A : No information (do not create routing "TRN" files). + ! v1 : Outlet locations file produced manually by Randy Koster. + ! v2 : Outlet locations file produced by run_routing_raster.py using routing information encoded + ! in SRTM-based Pfafstetter catchments and Greenland outlets info provided by Lauren Andrews. + implicit none character(*), intent (in) :: LBCSV ! land BCs version @@ -103,6 +110,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'NGDC' MODALB = 'MODIS1' SNOWALB = 'LUT' + OUTLETV = "N/A" GNU = 2.17 use_PEATMAP = .false. jpl_height = .false. @@ -112,6 +120,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'NGDC' MODALB = 'MODIS2' SNOWALB = 'LUT' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. @@ -121,6 +130,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'LUT' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. @@ -130,6 +140,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'LUT' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .false. jpl_height = .true. @@ -139,6 +150,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'LUT' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. @@ -148,6 +160,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'MODC061' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. @@ -157,6 +170,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'LUT' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .true. jpl_height = .false. @@ -166,6 +180,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'MODC061' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .false. jpl_height = .false. @@ -175,6 +190,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'MODC061' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .true. jpl_height = .false. @@ -184,6 +200,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'MODC061v2' + OUTLETV = "N/A" GNU = 1.0 use_PEATMAP = .true. jpl_height = .false. @@ -193,6 +210,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD' MODALB = 'MODIS2' SNOWALB = 'MODC061v2' + OUTLETV = "v1" GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. @@ -202,6 +220,7 @@ SUBROUTINE init_bcs_config (LBCSV) SOILBCS = 'HWSD_b' MODALB = 'MODIS2' SNOWALB = 'MODC061v2' + OUTLETV = "v2" GNU = 1.0 use_PEATMAP = .true. jpl_height = .true. diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt index 15405d696..d4eca3cfc 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/CMakeLists.txt @@ -1 +1 @@ -esma_add_subdirectories (soil) +esma_add_subdirectories (soil routing) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt new file mode 100644 index 000000000..46bdc0dc2 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/CMakeLists.txt @@ -0,0 +1,30 @@ +esma_set_this () + +set(srcs + routing_constant.f90 +) + +set (exe_srcs + get_finalID_msk.f90 + get_landocean_Greenland_real.f90 + get_outlets_land_allcat.f90 + get_sinkxy_land.f90 + get_outlets_catchindex.f90 + get_outlets_land.f90 + Pfaf_to_2d_30s_land.f90 +) + +esma_add_library (${this} + SRCS ${srcs} +) +foreach (src ${exe_srcs}) + string (REGEX REPLACE ".f90" ".x" exe ${src}) + ecbuild_add_executable ( + TARGET ${exe} + SOURCES ${src} + NOINSTALL + LIBS ${this} NetCDF::NetCDF_Fortran MPI::MPI_Fortran) +endforeach () + +# copy to the build directory +configure_file(run_routing_raster.py run_routing_raster.py) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/Pfaf_to_2d_30s_land.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/Pfaf_to_2d_30s_land.f90 new file mode 100644 index 000000000..3935fbb3d --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/Pfaf_to_2d_30s_land.f90 @@ -0,0 +1,71 @@ +program main + + use routing_constant,only : nc,ng,nlon,nlat + implicit none + + character(len=100) :: var1="outlet_sinky_allcat" + character(len=100) :: var2="outlet_sinkx_allcat" + character(len=100) :: map="Pfafstetter_Greenland_real" + + real*8,allocatable :: lon(:),lat(:) + integer,allocatable :: catchind(:,:) + integer,allocatable :: lons(:,:),lats(:,:) + integer,allocatable :: data_Pfaf(:) + + integer :: i,j,xi,yi,id,nall + +! Transform the 1d list to the unformatted Fortran binary file "Outlet_latlon.43200x21600" that can be read directly by "mk_runofftbl.F90" of makebcs. +! nc= number of land catchments (excluding Greenland), 291284 in our case. +! ng= number of Greenland catchments, 525 in our case +! nall =number of the total catchments (including Greenland) +! catchind = catchment index: 1-291284 for land catchments, and 291285-291809 for Greenland catchments +! data_Pfaf= lat (or lon) index (on a 30s map) of the final sink point for each catchment, +! lats= a map (with a resolution of 30s) of lat values. For each pixel, the value is the latitude of its final outlet point +! lons= a map (with a resolution of 30s) of lon values. For each pixel, the value is the longitude of its final outlet point + + nall=nc+ng + allocate(catchind(nlon,nlat),lons(nlon,nlat),lats(nlon,nlat)) + allocate(lon(nlon),lat(nlat)) + +! Read the raster array of catchment indices + open(30,file="outputs/"//trim(map),form="unformatted") + do j = 1,nlat + read (30) catchind(:,j) + end do + + allocate(data_Pfaf(nall)) +! Read in the latitudes associated with each catchment + open(77,file="outputs/"//trim(var1)//".txt") + read(77,*)data_Pfaf + lats=-999 +! For each raster point, find the catchment index and use that, in conjunction with data_Pfaf, to compute a raster array of latitudes + do xi=1,nlon + do yi=1,nlat + if(catchind(xi,yi)>=1.and.catchind(xi,yi)<=nall)then + id=catchind(xi,yi) + lats(xi,yi)=data_Pfaf(id) + endif + enddo + enddo + +! Read in the longitudes associated with each catchment + open(77,file="outputs/"//trim(var2)//".txt") + read(77,*)data_Pfaf + lons=-999 + ! For each raster point, find the catchment index and use that, in conjunction with data_Pfaf, to compute a raster array of longitudes + do xi=1,nlon + do yi=1,nlat + if(catchind(xi,yi)>=1.and.catchind(xi,yi)<=nall)then + id=catchind(xi,yi) + lons(xi,yi)=data_Pfaf(id) + endif + enddo + enddo + + open(30,file="Outlet_latlon.43200x21600",form="unformatted") + do j = 1, nlat + write (30) lons(:,j) + write (30) lats(:,j) + end do + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_finalID_msk.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_finalID_msk.f90 new file mode 100644 index 000000000..9f65bfa9b --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_finalID_msk.f90 @@ -0,0 +1,176 @@ +program main + + use routing_constant,only : nc + implicit none + + integer,allocatable,dimension(:) :: downid,finalid + real*8,allocatable,dimension(:) :: pfaf + integer,allocatable,dimension(:,:) :: pfaf_digit + integer*8,allocatable,dimension(:) :: res + integer,allocatable,dimension(:) :: pfaf_last,pfaf_msk,code,behind + integer,allocatable,dimension(:) :: first,last + + integer :: i,j,jj,k,p,down,cur,idx,num,ok,samed + integer :: fulli(12),fullj(12) + real :: val(9) + + character(len=100) :: file_path !/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat + +! Get downstream catchment and final destination ID for each catchment, and determine whether it directs to an ocean or inland lake. +! downid=Pfafstetter index of catchment just downstream +! finalid=Pfafstetter index of catchment at outlet point +! pfaf= Pfafstetter number for catchment +! pfaf_digit= The 12 digits in a Pfafstetter number, separated +! pfaf_last= The index of the last nonzero digit in a Pfafstetter number (counting from the left) +! pfaf_msk =1 for non-sink catchments, 2 for sink catchments with endpoints in ocean, =3 for sink catchments with endpoints in interior lake +! last= The index of the last digit in a Pfafstetter number after removing any 11..000 tail. +! first= The index of the last zero (but not the zero at the very end). However, if there are no zeroes until the end, first =2 (the second index, since the first index indicates the continent). + + if (command_argument_count() /= 1) then + print *, "no found" + stop + endif + call get_command_argument(1, file_path) + + open(77,file=file_path, form="formatted", status="old") + read(77,*)num + + allocate(downid(nc),finalid(nc),pfaf(nc),pfaf_digit(nc,12),res(nc),pfaf_last(nc),pfaf_msk(nc)) + allocate(first(nc),last(nc)) + + do i=1,nc + read(77,*)idx,pfaf(i) + enddo + +! Separate Pfafstetter number into individual digits + res=int8(pfaf) + pfaf_digit(:,1)=res/(int8(10)**int8(11)) + do i=2,12 + res=res-int8(10)**int8(13-i)*int8(pfaf_digit(:,i-1)) + pfaf_digit(:,i)=res/(int8(10)**int8(12-i)) + enddo + +! Determine positions of last nonzero digit (pfaf_last) and the last digit that’s neither 0 nor 1 (at the end) + first=2 + last=2 + do i=1,nc + do j=12,1,-1 + if(pfaf_digit(i,j)/=0)then + pfaf_last(i)=j + do k=0,j-1 + if(pfaf_digit(i,j-k)/=1)then + last(i)=j-k + exit + endif + enddo + exit + endif + enddo + enddo + do i=1,nc + if(last(i)<=1) last(i)=2 + enddo + +! Determine position of final zero that has some nonzero digits after it + do i=1,nc + do j=last(i),2,-1 + if(pfaf_digit(i,j)==0)then + first(i)=j + exit + endif + enddo + enddo + + do i=1,nc + + if(first(i)>last(i)-1)then + downid(i)=-1 + else + + allocate(code(1:last(i)-first(i))) + code=pfaf_digit(i,first(i):last(i)-1) + if(any(code==2).or.any(code==4).or.any(code==6).or.any(code==8))then + ! If all digits (after the first) are odd, the Pfafstetter logic implies that the catchment will be on the coast. + fulli=pfaf_digit(i,:) + do j=i-1,1,-1 ! Test each catchment to see if it lies just downstream of catchment i + ok=1 + fullj=pfaf_digit(j,:) + samed=0 + do k=1,min(pfaf_last(i),pfaf_last(j)) ! Determine the index (samed) up to which the Pfaf numbers of catchment I and j match + if(fulli(k)==fullj(k))then + samed=samed+1 + else + exit + endif + enddo ! end k loop + if(samed+1<=pfaf_last(j))then + ! Check that none of catchment j’s indices (after samed) are even, which would imply a downstream branching off from the river on which catchment i lies. + allocate(behind(1:pfaf_last(j)-samed)) + behind=fullj(samed+1:pfaf_last(j)) + if(any(mod(behind,2)==0)) ok=0 + deallocate(behind) + else + ok=0 + endif + if(ok==1)then + downid(i)=j + exit + endif + enddo ! end j loop + else + downid(i)=-1 + endif + deallocate(code) + + endif ! end i loop + + + enddo + + + open(88,file="outputs/Pfaf_downid.txt") + do i=1,nc + write(88,*)downid(i) + enddo + +! Keep “moving downstream” until you find the catchment with no downstream catchment: + do i=1,nc + cur=i + down=downid(i) + do while(down/=-1) + cur=down + down=downid(cur) + enddo + finalid(i)=cur + enddo + + open(88,file="outputs/Pfaf_finalID.txt") + do i=1,nc + write(88,*)finalid(i) + enddo + + + ! Set masks: 1 = has downstream catchment; 2 = drains to ocean; 3 = drains to inland lake + do i=1,nc + if(downid(i)/=-1)then + pfaf_msk(i)=1 + else + allocate(code(1:pfaf_last(i))) + code=pfaf_digit(i,1:pfaf_last(i)) + if(any(code==0))then + pfaf_msk(i)=3 + else + pfaf_msk(i)=2 + endif + deallocate(code) + end if + enddo + + open(88,file="outputs/Pfaf_msk.txt") + do i=1,nc + write(88,*)pfaf_msk(i) + enddo + + + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_landocean_Greenland_real.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_landocean_Greenland_real.f90 new file mode 100644 index 000000000..63da40671 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_landocean_Greenland_real.f90 @@ -0,0 +1,76 @@ +program main + + use routing_constant, only : nc,nlon,nlat,nlon_G,nlat_G,loni_min,loni_max,lati_min,lati_max,& + nlon1m,nlat1m + + implicit none + include 'netcdf.inc' + + real*8,allocatable,dimension(:) :: lon_G,lat_G + integer,allocatable,dimension(:,:) :: catchind,landocean,Greenland + integer,allocatable,dimension(:) :: Pfaf_real, countc + + integer :: i,j,ret,ncid,varid,ntile + real :: val(4) + + character(len=100) :: file_path1 !/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/topo/v1/SRTM-TopoData/SRTM_PfafData.nc + character(len=100) :: file_path2 !/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/land/routing/GreenlandID_30s.nc + +! Insert the Greenland index map into the catchment index map. + + if (command_argument_count() /= 2) then + print *, "no found" + stop + endif + call get_command_argument(1, file_path1) + call get_command_argument(2, file_path2) + + allocate(catchind(nlon1m,nlat1m)) + ret=nf_open(file_path1,0,ncid) + ret=nf_inq_varid(ncid,"CatchIndex",varid) + ret=nf_get_var_int(ncid,varid,catchind) + ret=nf_close(ncid) + allocate(landocean(nlon,nlat)) + landocean=-9999 + do i=1,nlon1m + do j=1,nlat1m + if(catchind(i,j)/=-9999)then + landocean(2*i-1:2*i,2*j-1:2*j)=catchind(i,j) + endif + enddo + enddo + + + allocate(Greenland(nlon_G,nlat_G)) + allocate(lon_G(nlon_G),lat_G(nlat_G)) + ret=nf_open(file_path2,0,ncid) + ret=nf_inq_varid(ncid,"lon",varid) + ret=nf_get_var_double(ncid,varid,lon_G) + ret=nf_close(ncid) + ret=nf_open(file_path2,0,ncid) + ret=nf_inq_varid(ncid,"lat",varid) + ret=nf_get_var_double(ncid,varid,lat_G) + ret=nf_close(ncid) + ret=nf_open(file_path2,0,ncid) + ret=nf_inq_varid(ncid,"data",varid) + ret=nf_get_var_int(ncid,varid,Greenland) + ret=nf_close(ncid) + + + where(Greenland/=-9999) landocean(loni_min:loni_max,lati_min:lati_max)=Greenland + + + do i=1,nlon + do j=1,nlat + if(landocean(i,j)>=700000000)then + landocean(i,j)=landocean(i,j)-700000000+nc + endif + enddo + enddo + + open(30,file="outputs/Pfafstetter_Greenland_real",form="unformatted") + do j = 1,nlat + write (30) landocean(:,j) + end do + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_catchindex.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_catchindex.f90 new file mode 100644 index 000000000..d2f69a12e --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_catchindex.f90 @@ -0,0 +1,60 @@ +program main + + use routing_constant,only : nc,nl,ng + implicit none + + integer,allocatable,dimension(:) :: msk,outid,mskall,final,finalall + + integer :: k,i,ntot,ns + +! Get a list of sink catchment IDs (both land and Greenland) +! nc = number of the land catchments (excluding Greenland), 291284 in our case. +! ng = number of the Greenland catchments, 525 in our case. +! nl = number of outlets to ocean in land (not including Greenland outlets) +! ns = number of the total outlets (including Greenland); note that all the Greenland catchments are sink catchments (ie, no downstream catchment). +! ntot = number of the total catchments (including Greenland) +! outid = a list of indices for sink catchments (to ocean) + + ntot=nc+ng + ns=nl+ng + allocate(msk(nc),outid(ns),mskall(ntot),final(nc),finalall(ntot)) + open(77,file="outputs/Pfaf_msk.txt") !!! + read(77,*)msk + k=0 + do i=1,nc + if(msk(i).eq.2)then ! msk=2 is for a catchment that drains directly into the ocean + k=k+1 + outid(k)=i + end if + end do + + ! Add Greenland catchments to list of outlet catchments; write outlet catchments to file + do i=k+1,ns + outid(i)=nc+i-k + end do + open(88,file="outputs/outlet_catchindex.txt") + do i=1,ns + write(88,*)outid(i) + enddo + +! Append msk values for Greenland to original msk array; write out + mskall(1:nc)=msk + mskall(nc+1:)=2 + open(88,file="outputs/Pfaf_msk_all.txt") + do i=1,ntot + write(88,*)mskall(i) + enddo + +! Append final values for Greenland to original “final” array; write out + open(77,file="outputs/Pfaf_finalID.txt") + read(77,*)final + finalall(1:nc)=final + do i=nc+1,ntot + finalall(i)=i + end do + open(88,file="outputs/Pfaf_finalID_all.txt") + do i=1,ntot + write(88,*)finalall(i) + enddo + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land.f90 new file mode 100644 index 000000000..a8db67847 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land.f90 @@ -0,0 +1,120 @@ +program main + + use routing_constant,only : nc,nl,ng,nlon=>nlon1m,nlat=>nlat1m + implicit none + include 'netcdf.inc' + + real*8,allocatable :: lon(:),lat(:),long(:),latg(:),lons(:),lats(:) + integer,allocatable :: catchind(:,:) + real,allocatable :: acah(:,:) + integer,allocatable :: down(:),sx(:),sy(:),msk(:) + real,allocatable :: acas(:) + + integer :: id,xi,yi,i,k,xis,yis,ntot,ncid,ret,varid + + character(len=100) :: file_path1 !/discover/nobackup/projects/gmao/bcs_shared/make_bcs_inputs/land/topo/v1/SRTM-TopoData/SRTM_PfafData.nc + character(len=100) :: file_path2 !/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/land/routing/HydroSHEDS_drainage_area.nc + character(len=100) :: file_path3 !/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/land/routing/Greenland_outlets_lat.txt + character(len=100) :: file_path4 !/discover/nobackup/projects/gmao/bcs_shared/preprocessing_bcs_inputs/land/routing/Greenland_outlets_lon.txt + + if (command_argument_count() /= 4) then + print *, "no found" + stop + endif + call get_command_argument(1, file_path1) + call get_command_argument(2, file_path2) + call get_command_argument(3, file_path3) + call get_command_argument(4, file_path4) + +! Get sink points on land or in Greenland (from Lauren Andrews) by picking the point (i.e., 1 minute grid cell) within each sink catchment that has the largest drainage area per the HydroSHEDS (https://www.hydrosheds.org/) dataset. +! acah = map of drainage area with a resolution of 1m from HydroSHEDS. +! down = catchment index of the downstream catchment. +! sx = the lon index (on a 1m map) of the outlet point (>0 only for sink catchments). +! sy = the lat index (on a 1m map) of the outlet point (>0 only for sink catchments). +! msk: 1 = has downstream catchment; 2 = drains to ocean; 3 = drains to inland lake +! acas = maximum drainage area, defined (>0) only for sink catchments +! ntot = number of the total outlets (including Greenland) +! nl = number of outlets to ocean in land (not including Greenland) +! ng = number of outlets to ocean in Greenland + + ntot=nl+ng + allocate(catchind(nlon,nlat),acah(nlon,nlat)) + allocate(lon(nlon),lat(nlat)) + allocate(sx(nc),sy(nc),acas(nc),down(nc),msk(nc)) + allocate(long(ng),latg(ng),lons(ntot),lats(ntot)) + + + ret=nf_open(file_path1,0,ncid) + ret=nf_inq_varid(ncid,"longitude",varid) + ret=nf_get_var_double(ncid,varid,lon) + ret=nf_close(ncid) + ret=nf_open(file_path1,0,ncid) + ret=nf_inq_varid(ncid,"latitude",varid) + ret=nf_get_var_double(ncid,varid,lat) + ret=nf_close(ncid) + + ret=nf_open(file_path1,0,ncid) + ret=nf_inq_varid(ncid,"CatchIndex",varid) + ret=nf_get_var_int(ncid,varid,catchind) + ret=nf_close(ncid) + + ret=nf_open(file_path2,0,ncid) + ret=nf_inq_varid(ncid,"data",varid) + ret=nf_get_var_real(ncid,varid,acah) + ret=nf_close(ncid) + + open(77,file="outputs/Pfaf_downid.txt") + read(77,*)down + open(77,file="outputs/Pfaf_msk.txt") + read(77,*)msk + +! For each long/lat location, determine if the catchment holding it is an outlet catchment to the ocean, +! and if so, determine if this point has the maximum drainage area + acas=-9999. + sx=0 + sy=0 + do xi=1,nlon + do yi=1,nlat + if(catchind(xi,yi)>=1)then + id=catchind(xi,yi) + if(down(id)==-1.and.acah(xi,yi)>=acas(id))then + acas(id)=acah(xi,yi) + sx(id)=xi + sy(id)=yi + endif + endif + enddo + enddo + +! Construct arrays of longitudes and latitudes of sink points. + where(down/=-1)sx=-1 + where(down/=-1)sy=-1 + k=0 + do i=1,nc + if(msk(i)==2)then + k=k+1 + lons(k)=lon(sx(i)) + lats(k)=lat(sy(i)) + endif + enddo + +! Append Greenland values to the longitude and latitude arrays. + open(77,file=file_path3) + read(77,*)latg + open(77,file=file_path4) + read(77,*)long + + lons(k+1:ntot)=long + lats(k+1:ntot)=latg + +! Write out arrays of sink point longitudes and latitudes + open(88,file="outputs/outlet_sinklat.txt") + do i=1,ntot + write(88,*)lats(i) + enddo + open(88,file="outputs/outlet_sinklon.txt") + do i=1,ntot + write(88,*)lons(i) + enddo + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land_allcat.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land_allcat.f90 new file mode 100644 index 000000000..49b6207f7 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land_allcat.f90 @@ -0,0 +1,66 @@ +program main + + use routing_constant,only : nc,nl,ng + implicit none + + integer, allocatable, dimension(:) :: id_final,id_outlet,msk + integer,allocatable,dimension(:) :: lati_outlet,loni_outlet + integer,allocatable,dimension(:) :: lati_full,loni_full + + integer :: i,j,nall,ns + +! Assign outlet locations to all upstream catchments to create a 1d list showing the final x and y indexes for each catchment. +! id_final = index of the final sink catchment for each catchment +! id_outlet =catchment index for each outlet +! msk = 1 = has downstream catchment; 2 = drains to ocean; 3 = drains to inland lake +! lati_outlet = lat index for each outlet +! loni_outlet =lon index for each outlet +! lati_full = lat index (on a 30s map) of its final sink point for each catchment +! loni_full = lon index (on a 30s map) of its final sink point for each catchment + + nall=nc+ng + ns=nl+ng + allocate(id_final(nall),id_outlet(ns),msk(nall),& + lati_outlet(ns),loni_outlet(ns),lati_full(nall),loni_full(nall)) + + open(77,file="outputs/Pfaf_finalID_all.txt") + read(77,*)id_final + open(77,file="outputs/outlet_catchindex.txt") + read(77,*)id_outlet + open(77,file="outputs/outlet_sinky.txt") + read(77,*)lati_outlet + open(77,file="outputs/outlet_sinkx.txt") + read(77,*)loni_outlet + open(77,file="outputs/Pfaf_msk_all.txt") + read(77,*)msk + + lati_full=-999 + loni_full=-999 + + do i=1,nall +! For each catchment, check to see if its final sink catchment actually drains to the ocean. + if(msk(id_final(i)).eq.2)then + ! For the catchment being tested, loop over the indices of all sink catchments and find the one that matches its own sink catchment. + ! Assign lat/longs of outlet point to the catchment being tested. + do j=1,ns + if(id_outlet(j).eq.id_final(i))then + lati_full(i)=lati_outlet(j) + loni_full(i)=loni_outlet(j) + end if + enddo + else if(msk(id_final(i)).eq.3)then + lati_full(i)=-999 + loni_full(i)=-999 + endif + end do + + open(88,file="outputs/outlet_sinky_allcat.txt") + do i=1,nall + write(88,*)lati_full(i) + enddo + open(88,file="outputs/outlet_sinkx_allcat.txt") + do i=1,nall + write(88,*)loni_full(i) + enddo + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_sinkxy_land.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_sinkxy_land.f90 new file mode 100644 index 000000000..41d830888 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_sinkxy_land.f90 @@ -0,0 +1,64 @@ +program main + + use routing_constant,only : nl,ng,nlon,nlat + implicit none + + real*8,allocatable,dimension(:) :: lats,lons,lat30s,lon30s,lat_dis,lon_dis + integer,allocatable,dimension(:) :: lati,loni + + integer :: i,temp(1),ns + real*8 :: dlat,dlon + +! Convert outlet locations in degree lat/lon to indices on the 30 arc-sec raster grid. +! lats = lat degree of the outlets +! lons = lon degree of the outlets +! lati = lat index of the outlets on the 30s map +! loni =lon index of the outlets on the 30s map +! lat30s = lat coordinate of the 30s map +! lon30s =lon coordinate of the 30s map +! lat_dis = distance of lat between the center of each 30s pixel and the outlet point. +! lon_dis = distance of lon between the center of each 30s pixel and the outlet point. +! nlat = number of latitude indices. For 30s map, nlat = 21600 +! nlon = number of longitude indices. For 30s map, nlon = 43200 + + ns=nl+ng + allocate(lats(ns),lons(ns),lati(ns),loni(ns)) + allocate(lat30s(nlat),lon30s(nlon),lat_dis(nlat),lon_dis(nlon)) + open(77,file="outputs/outlet_sinklat.txt") + read(77,*)lats + open(77,file="outputs/outlet_sinklon.txt") + read(77,*)lons + + dlat=180.D0/nlat + dlon=360.D0/nlon + lat30s(1)=-90.D0+dlat/2.D0 + lon30s(1)=-180.D0+dlon/2.D0 + do i=2,nlat + lat30s(i)=lat30s(i-1)+dlat + enddo + do i=2,nlon + lon30s(i)=lon30s(i-1)+dlon + enddo + +! For each sink catchment, find the 30s-latitude and 30s-longitude closest to its outlet point. + do i=1,ns + lat_dis=abs(lat30s-lats(i)) + temp=minloc(lat_dis) + lati(i)=temp(1) + enddo + do i=1,ns + lon_dis=abs(lon30s-lons(i)) + temp=minloc(lon_dis) + loni(i)=temp(1) + enddo + + open(88,file="outputs/outlet_sinky.txt") + do i=1,ns + write(88,*)lati(i) + enddo + open(88,file="outputs/outlet_sinkx.txt") + do i=1,ns + write(88,*)loni(i) + enddo + +end program main diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/readme.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/readme.txt new file mode 100644 index 000000000..c47c3b3dc --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/readme.txt @@ -0,0 +1,35 @@ +v1 11/08/2023, Yujin Zeng + +The "preproc/routing" package is used for creating a 30-arcsec raster file with the locations of river outlets to the ocean. + +The output from this package is a binary file "Outlet_latlon.43200x21600". + +The river outlets are located in land or landice tiles as defined in the file "SRTM_PfafData.nc". + +The "Outlet_latlon.43200x21600" file is the input for "mk_runofftbl.F90" in the makebcs package, which further adjusts the outlet locations to be consistent with the ocean model resolution and domain ("mk_runofftbl.F90"). + +If on NCCS/Discover, the package can be run using the script "run_routing_raster.py". Users should source g5_modules before run to get the necessary env. If not on Discover, please contact yujin.zeng@nasa.gov. + +The tasks completed by each f90 program are briefly described as follows: + +1. get_finalID_msk.f90: +Get downstream catchment and final destination ID for each catchment, and determine whether it directs to ocean or inland lake. + +2. get_outlets_catchindex.f90: +Get a list of sink catchment IDs. + +3. get_outlets_land.f90: +Get sink points on land or in Greenland (from Lauren Andrews) by picking the point (i.e., 15-arcsec grid cell) within each sink catchment that has the largest drainage area per the HydroSHEDS (https://www.hydrosheds.org/) dataset. + +4. get_sinkxy_land.f90: +Convert outlet locations in degree lat/lon to indices on the 30 arc-sec raster grid. + +5. get_outlets_land_allcat.f90: +Assign outlet locations to all upstream catchments to create a 1d list showing the final x and y indexes for each catchment. + +6. get_landocean_Greenland_real.f90: +Insert the Greenland index map into the catchment index map. + +7. Pfaf_to_2d_30s_land.f90: +Transform the 1d list above to the unformatted Fortran binary file "Outlet_latlon.43200x21600" that can be read directly by "mk_runofftbl.F90" of makebcs. + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 new file mode 100644 index 000000000..c4fbfa3cb --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/routing_constant.f90 @@ -0,0 +1,22 @@ +module routing_constant + + ! hardwired constants for GEOS river routing scheme + + implicit none + public + + integer,parameter :: nlon = 43200 ! number of lat of world grid with 30 sec resolution + integer,parameter :: nlat = 21600 ! number of lon of world grid with 30 sec resolution + integer,parameter :: nlon1m = 21600 ! number of lat of world grid with 1 min resolution + integer,parameter :: nlat1m = 10800 ! number of lon of world grid with 1 min resolution + integer,parameter :: nlon_G = 8400 ! number of lat of Greenland grid (30 sec resolution) + integer,parameter :: nlat_G = 4800 ! number of lon of Greenland grid (30 sec resolution) + integer,parameter :: loni_min = 12001 ! ind of lon start of Greenland grid in world grid (30 sec resolution) + integer,parameter :: loni_max = 20400 ! = loni_min + nlon_G - 1, ind of lon end of Greenland grid in world grid (30 sec resolution) + integer,parameter :: lati_min = 16801 ! ind of lat start of Greenland grid in world grid (30 sec resolution) + integer,parameter :: lati_max = 21600 ! = lati_min + nlat_G - 1, ind of lat end of Greenland grid in world grid (30 sec resolution) + integer,parameter :: nc = 291284 ! number of catchments in land + integer,parameter :: ng = 525 ! number of catchments in Greenland + integer,parameter :: nl = 22116 ! number of outlets to ocean in land (not including Greenland) + +end module routing_constant diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/run_routing_raster.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/run_routing_raster.py new file mode 100755 index 000000000..c76a06062 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/run_routing_raster.py @@ -0,0 +1,62 @@ +#!/usr/bin/env python3 + +#source g5_modules before run to get the necessary env + +import os +import subprocess + +# Path to "bcs_shared" directory: + +file_path = "/discover/nobackup/projects/gmao/bcs_shared/" # NCCS Discover +#file_path = "/nobackup/gmao_SIteam/ModelData/bcs_shared/" # NAS + +# Input files + +file_Pfafcatch = file_path + "/make_bcs_inputs/land/topo/v1/SRTM-TopoData/Pfafcatch-routing.dat" +file_SRTMPfaf = file_path + "/make_bcs_inputs/land/topo/v1/SRTM-TopoData/SRTM_PfafData.nc" +file_Drainage = file_path + "/preprocessing_bcs_inputs/land/routing/HydroSHEDS_drainage_area.nc" +file_GrnLat = file_path + "/preprocessing_bcs_inputs/land/routing/Greenland_outlets_lat.txt" +file_GrnLon = file_path + "/preprocessing_bcs_inputs/land/routing/Greenland_outlets_lon.txt" +file_GrnMap = file_path + "/preprocessing_bcs_inputs/land/routing/GreenlandID_30s.nc" + +name_Pfafcatch = os.path.basename(file_Pfafcatch) +name_SRTMPfaf = os.path.basename(file_SRTMPfaf) +name_Drainage = os.path.basename(file_Drainage) +name_GrnLat = os.path.basename(file_GrnLat) +name_GrnLon = os.path.basename(file_GrnLon) +name_GrnMap = os.path.basename(file_GrnMap) + +files=[file_Pfafcatch,file_SRTMPfaf,file_Drainage,file_GrnLat,file_GrnLon,file_GrnMap] +# Remove files and directories +os.system("rm -rf inputs >& /dev/null") +os.system("rm -rf outputs >& /dev/null") +os.system("rm -f Outlet_latlon.43200x21600 >& /dev/null") + +# Create directories and symbolic links +os.makedirs("inputs", exist_ok=True) +os.makedirs("outputs", exist_ok=True) +for file in files: + file_name = os.path.basename(file) + os.symlink(file, os.path.join("inputs", file_name)) + +program_inputs = { + "./get_finalID_msk.x": [f"inputs/{name_Pfafcatch}"], + "./get_outlets_catchindex.x": [], + "./get_outlets_land.x": [f"inputs/{name_SRTMPfaf}", f"inputs/{name_Drainage}", f"inputs/{name_GrnLat}", f"inputs/{name_GrnLon}"], + "./get_sinkxy_land.x": [], + "./get_outlets_land_allcat.x": [], + "./get_landocean_Greenland_real.x": [f"inputs/{name_SRTMPfaf}", f"inputs/{name_GrnMap}"], + "./Pfaf_to_2d_30s_land.x": [] +} + +for program, input_files in program_inputs.items(): + print(f"running {program}") + command = [program] + input_files + subprocess.run(command) + +print("Outlet_latlon.43200x21600 created!") + +# Clean up +print("Removing temporary output files ...") +os.system("rm -rf outputs") +os.system("rm -rf inputs")