small update for levee test scripts

etc/levee_test/s03-analysis.sh → etc/levee_test/s02-analysis.sh

@@ -2,30 +2,43 @@
 
 # link output directory
 
-ln -sf ../../out/levee_test       out_lev
+ln -sf ../../out/levee_sample     out_lev
 ln -sf ../../out/test1-glb_15min  out_ori
+ln -sf ../../map/glb_15min        map
 
 mkdir -p data
 mkdir -p fig
 
 # some analysis on river depth & flood fraction
 
-./src/analysis_output
+#./src/analysis_output
 
 ###########
 
 # Downscale flood depth for Mississippi
 
-WEST=-93
+WEST=-94
 EAST=-86
-SOUTH=29
-NORTH=35
-HIRES="3sec"
-#HIRES="1min"
+SOUTH=34
+NORTH=40
+#HIRES="3sec"
+HIRES="1min"
+
+LEVFRC="map/lev_frc_global.bin"
+
+############
+
+rm -f dph_lev.bin
+rm -f dph_ori.bin
+rm -f flood_lev.bin
+rm -f flood_ori.bin
+rm -f hand.bin
+rm -f protect_pix.bin
 
 # define levee protected pixels 
-echo "./src/define_protect_pix  $WEST $EAST $SOUTH $NORTH $HIRES"
-./src/define_protect_pix  $WEST $EAST $SOUTH $NORTH $HIRES
+echo "./src/define_protect_pix  $WEST $EAST $SOUTH $NORTH $HIRES $LEVFRC"
+./src/define_protect_pix  $WEST $EAST $SOUTH $NORTH $HIRES  $LEVFRC
+
 
 # downscale flood depth with levee
 

etc/levee_test/src/define_protect_pix.F90

@@ -27,11 +27,12 @@ program downscale_flddph
 !
       integer*4,allocatable    ::  nextXX(:,:)                    !! downstream (jXX,jYY)
       integer*2,allocatable    ::  catmXX(:,:), catmYY(:,:)        !! catchment (iXX,iYY) of pixel (ix,iy)
+      integer*1,allocatable    ::  catmZZ(:,:)                    !! catchment Z layer (100 levels)
+
       real,allocatable         ::  flddif(:,:), rivwth(:,:), hand(:,:)      !! height above channel [m]
       real,allocatable         ::  lon(:), lat(:)
 
-      real,allocatable         ::  levbas(:,:)                   !! flood depth [m] (coarse resolution)
-      real,allocatable         ::  levhgt(:,:)                   !! flood depth [m] (coarse resolution)
+      real,allocatable         ::  levfrc(:,:)                   !! levee unprotected fraction
 !
       real,allocatable         ::  protect(:,:)                    !! downscaled flood depth [m]
       real,allocatable         ::  flddif2(:,:),hand2(:,:)                    !! downscaled flood depth [m]
@@ -40,7 +41,7 @@ program downscale_flddph
       character*256            ::  mapdir, hires
       parameter                   (mapdir='./map/')              !! map directory (please make a symbolic link)
       character*256            ::  fnextxy, fflood, rfile
-      character*256            ::  flevbas, flevhgt
+      character*256            ::  flevfrc
       character*256            ::  buf
       integer                  ::  ios
 ! ===============================================
@@ -56,6 +57,8 @@ program downscale_flddph
 
       call getarg(5,hires)       !! downscale resolution
 
+      call getarg(6,flevfrc)       !! downscale resolution
+
       param=trim(mapdir)//'params.txt'
 
       open(11,file=param,form='formatted')
@@ -102,25 +105,19 @@ program downscale_flddph
 
 ! ==========
 
-      allocate(nextXX(nXX,nYY),levbas(nXX,nYY),levhgt(nXX,nYY))
+      allocate(nextXX(nXX,nYY),levfrc(nXX,nYY))
       allocate(protect(mx,my),flddif2(mx,my),hand2(mx,my))
       protect(:,:)=-9999
 
 ! ===============================================
       fnextxy=trim(mapdir)//'nextxy.bin'
-      flevbas=trim(mapdir)//'levbas.bin'
-      flevhgt=trim(mapdir)//'levhgt.bin'
 
       open(11, file=fnextxy, form='unformatted', access='direct', recl=4*nXX*nYY)
       read(11,rec=1) nextXX
       close(11)
 
-      open(12, file=flevbas, form='unformatted', access='direct', recl=4*nXX*nYY)
-      read(12,rec=1) levbas
-      close(12)
-
-      open(12, file=flevhgt, form='unformatted', access='direct', recl=4*nXX*nYY)
-      read(12,rec=1) levhgt
+      open(12, file=flevfrc, form='unformatted', access='direct', recl=4*nXX*nYY)
+      read(12,rec=1) levfrc
       close(12)
 
 
@@ -134,7 +131,7 @@ program downscale_flddph
         read(11,*) buf, area, lon_ori, lon_end, lat_end, lat_ori, nx, ny, csize
 
         if( lon_end<west2 .or. lon_ori>east2 .or.lat_ori<south2 .or. lat_end>north2 ) goto 1090  !! out of domain
-        allocate(catmXX(nx,ny),catmYY(nx,ny),flddif(nx,ny),rivwth(nx,ny),hand(nx,ny))
+        allocate(catmXX(nx,ny),catmYY(nx,ny),catmZZ(nx,ny),flddif(nx,ny),rivwth(nx,ny),hand(nx,ny))
         allocate(lon(nx),lat(ny))
 
         rfile=trim(mapdir)//trim(hires)//'/'//trim(area)//'.catmxy.bin'
@@ -149,7 +146,19 @@ program downscale_flddph
           print *, 'no data: ', rfile
           stop
         endif
-
+
+        rfile=trim(mapdir)//trim(hires)//'/'//trim(area)//'.catmz100.bin'
+        print *, rfile
+        open(21,file=rfile,form='unformatted',access='direct',recl=1*nx*ny,status='old',iostat=ios)
+        if( ios==0 )then
+          read(21,rec=1) catmZZ
+          close(21)
+        else
+          print *, '*******************'
+          print *, 'no data: ', rfile
+          stop
+        endif
+
         rfile=trim(mapdir)//trim(hires)//'/'//trim(area)//'.flddif.bin'
         open(21,file=rfile,form='unformatted',access='direct',recl=4*nx*ny,status='old',iostat=ios)
         if( ios==0 )then
@@ -209,12 +218,11 @@ program downscale_flddph
               endif
 
               if( catmXX(ix,iy)>0 )then
-                protect(jx,jy)=0
+                protect(jx,jy)=0                                                    !! 
                 iXX=catmXX(ix,iy)
                 iYY=catmYY(ix,iy)
-                if( flddif(ix,iy)>levbas(iXX,iYY) .and. levbas(iXX,iYY)>0 )then
+                if( catmZZ(ix,iy)>levfrc(iXX,iYY)*100 )then  !! protected pixel
                   protect(jx,jy)=5
-                  if( flddif(ix,iy)>levhgt(iXX,iYY) ) protect(jx,jy)=6
                 endif
 
                 if( rivwth(ix,iy)/=-9999 .and. rivwth(ix,iy)/=0 )then !! permanent water
@@ -227,7 +235,7 @@ program downscale_flddph
           end do
         end do
 
-        deallocate(catmXX,catmYY,flddif,rivwth,hand)
+        deallocate(catmXX,catmYY,catmZZ,flddif,rivwth,hand)
         deallocate(lon,lat)
 
  1090   continue