Routing Outlets in Boundary Conditions

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/Pfaf_to_2d_30s_land.f90

@@ -0,0 +1,59 @@
+program main
+
+use omp_lib
+use rwncfile
+implicit none
+
+character(len=100) :: var1="outlet_sinky_allcat"
+character(len=100) :: var2="outlet_sinkx_allcat"
+character(len=100) :: map="Pfafstetter_Greenland_real.nc"
+integer,parameter :: nc=291809
+integer,parameter :: nlon=43200
+integer,parameter :: nlat=21600
+
+real*8,allocatable :: lon(:),lat(:)
+integer,allocatable :: catchind(:,:)
+integer,allocatable :: data2d(:,:)
+integer,allocatable :: data_Pfaf(:)
+
+integer :: xi,yi,id
+
+
+allocate(catchind(nlon,nlat),data2d(nlon,nlat))
+allocate(lon(nlon),lat(nlat))
+call read_ncfile_double1d("outputs/"//trim(map),"lon",lon,nlon)
+call read_ncfile_double1d("outputs/"//trim(map),"lat",lat,nlat)
+call read_ncfile_int2d("outputs/"//trim(map),"data",catchind,nlon,nlat)
+
+allocate(data_Pfaf(nc))
+
+open(77,file="outputs/"//trim(var1)//".txt")
+read(77,*)data_Pfaf
+data2d=-999
+do xi=1,nlon
+  do yi=1,nlat
+    if(catchind(xi,yi)>=1.and.catchind(xi,yi)<=nc)then
+      id=catchind(xi,yi)
+      data2d(xi,yi)=data_Pfaf(id)
+    endif 
+  enddo
+enddo
+call create_ncfile_int2d_fill("outputs/"//trim(var1)//"_2d.nc","data",data2d,lon,lat,nlon,nlat,-999.)
+
+open(77,file="outputs/"//trim(var2)//".txt")
+read(77,*)data_Pfaf
+data2d=-999
+do xi=1,nlon
+  do yi=1,nlat
+    if(catchind(xi,yi)>=1.and.catchind(xi,yi)<=nc)then
+      id=catchind(xi,yi)
+      data2d(xi,yi)=data_Pfaf(id)
+    endif 
+  enddo
+enddo
+call create_ncfile_int2d_fill("outputs/"//trim(var2)//"_2d.nc","data",data2d,lon,lat,nlon,nlat,-999.)
+
+
+
+
+end

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/build

@@ -0,0 +1,19 @@
+#!/bin/bash
+
+if [ $# -lt 1 ]; then
+        echo "no f90 specified"
+        exit
+fi
+
+string=$1
+array=(${string//./ })  
+
+FILENAME=${array[0]}
+
+
+NETCDF_PATH=/discover/nobackup/yzeng3/apps/netcdf-4.2.1.1
+
+ifort -qopenmp ncdioMod.f90 rwncMod.f90 ${FILENAME}.f90 -I$NETCDF_PATH/include -L$NETCDF_PATH/lib -L/usr/local/intel/oneapi/2021/compiler/2021.4.0/linux/lib -lnetcdf -lnetcdff -o ${FILENAME}.out
+
+
+

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_landocean_Greenland_real.f90

@@ -0,0 +1,66 @@
+program main
+
+use omp_lib
+use rwncfile
+implicit none
+
+integer,parameter :: nc=291809
+integer,parameter :: nlon=43200
+integer,parameter :: nlat=21600
+integer,parameter :: nlon_G=8400
+integer,parameter :: nlat_G=4800
+integer,parameter :: loni_min=12001
+integer,parameter :: loni_max=20400
+integer,parameter :: lati_min=16801
+integer,parameter :: lati_max=21600
+
+integer,parameter :: id_glac=290191
+integer,parameter :: id_lake=290190
+integer,parameter :: id_landend=290188
+
+real*8,allocatable,dimension(:) :: lon,lat,lon_G,lat_G
+integer,allocatable,dimension(:,:) :: landocean,Greenland
+integer,allocatable,dimension(:) :: Pfaf_real, countc
+
+integer :: i,j
+
+allocate(landocean(nlon,nlat))
+allocate(lon(nlon),lat(nlat))
+call read_ncfile_double1d("inputs/Pfafstetter.nc","lon",lon,nlon)
+call read_ncfile_double1d("inputs/Pfafstetter.nc","lat",lat,nlat)
+call read_ncfile_int2d("inputs/Pfafstetter.nc","data",landocean,nlon,nlat)
+
+allocate(Greenland(nlon_G,nlat_G))
+allocate(lon_G(nlon_G),lat_G(nlat_G))
+call read_ncfile_double1d("inputs/GreenlandID_30s.nc","lon",lon_G,nlon_G)
+call read_ncfile_double1d("inputs/GreenlandID_30s.nc","lat",lat_G,nlat_G)
+call read_ncfile_int2d("inputs/GreenlandID_30s.nc","data",Greenland,nlon_G,nlat_G)
+
+where(Greenland/=-9999.and.(landocean(loni_min:loni_max,lati_min:lati_max)<=id_landend.or.&
+	landocean(loni_min:loni_max,lati_min:lati_max)==id_glac ))&
+ landocean(loni_min:loni_max,lati_min:lati_max)=Greenland
+
+
+where(landocean>id_landend.and.landocean<id_lake) landocean=-9999
+where(landocean==id_lake.or.landocean==id_glac) landocean=0
+
+allocate(Pfaf_real(id_landend))
+open(77,file="inputs/Pfaf_real.txt")
+read(77,*)Pfaf_real
+
+do i=1,nlon
+  do j=1,nlat
+    if(landocean(i,j)<=id_landend.and.landocean(i,j)>=1)then
+      landocean(i,j)=Pfaf_real(landocean(i,j))
+    else if(landocean(i,j)>=700000000)then
+      landocean(i,j)=landocean(i,j)-700000000+291284
+    endif
+  enddo
+enddo
+
+
+
+call create_ncfile_int2d("outputs/Pfafstetter_Greenland_real.nc","data",landocean,lon,lat,nlon,nlat)
+
+
+end

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_catchindex.f90

@@ -0,0 +1,52 @@
+program main
+
+implicit none
+
+integer,parameter :: nc=291284
+integer,parameter :: ns=22612
+integer,parameter :: ng=525
+
+integer,allocatable,dimension(:) :: msk,outid,mskall,final,finalall
+
+integer :: k,i,ntot
+
+ntot=nc+ng
+allocate(msk(nc),outid(ns),mskall(ntot),final(nc),finalall(ntot))
+open(77,file="inputs/Pfaf_msk.txt")
+read(77,*)msk
+k=0
+do i=1,nc
+  if(msk(i).eq.2)then
+    k=k+1  	
+    outid(k)=i
+  end if
+end do
+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
+
+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
+
+open(77,file="inputs/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

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land.f90

@@ -0,0 +1,87 @@
+program main
+
+use omp_lib
+use rwncfile
+implicit none
+
+integer,parameter :: nc=291284
+integer,parameter :: nl=22087
+integer,parameter :: ng=525
+integer,parameter :: nlon=21600
+integer,parameter :: nlat=10800
+
+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
+
+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))
+
+call read_ncfile_double1d("inputs/CatchIndex.nc","lon",lon,nlon)
+call read_ncfile_double1d("inputs/CatchIndex.nc","lat",lat,nlat)
+call read_ncfile_int2d("inputs/CatchIndex.nc","data",catchind,nlon,nlat)
+call read_ncfile_real2d("inputs/HydroSHEDS_drainage_area.nc","data",acah,nlon,nlat)
+
+
+open(77,file="inputs/downstream_1D_new_noadj.txt")
+read(77,*)down
+open(77,file="inputs/Pfaf_msk.txt")
+read(77,*)msk
+
+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
+
+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
+!print *,k
+
+open(77,file="inputs/Greenland_outlets_lat.txt")
+read(77,*)latg
+open(77,file="inputs/Greenland_outlets_lon.txt")
+read(77,*)long
+
+lons(k+1:ntot)=long
+lats(k+1:ntot)=latg
+
+
+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

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_outlets_land_allcat.f90

@@ -0,0 +1,59 @@
+program main
+
+use omp_lib
+
+implicit none
+
+integer,parameter :: nall=291809
+integer,parameter :: nc=22612
+
+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
+
+allocate(id_final(nall),id_outlet(nc),msk(nall),&
+	lati_outlet(nc),loni_outlet(nc),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
+  !if(mod(i,1000)==0) print *,i
+  if(msk(id_final(i)).eq.2)then
+    do j=1,nc
+      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

GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/routing/get_sinkxy_land.f90

@@ -0,0 +1,46 @@
+program main
+
+use rwncfile
+implicit none
+
+integer,parameter :: ns=22612
+integer,parameter :: nlat=21600
+integer,parameter :: nlon=43200
+
+real*8,allocatable,dimension(:) :: lats,lons,lat30s,lon30s,lat_dis,lon_dis
+integer,allocatable,dimension(:) :: lati,loni
+
+integer :: i,temp(1)
+
+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
+open(77,file="inputs/lat_30s.txt")
+read(77,*)lat30s
+open(77,file="inputs/lon_30s.txt")
+read(77,*)lon30s
+
+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