diff --git a/etc/levee_test/s01-simulation_sample_param.sh b/etc/levee_test/s01-simulation_sample_param.sh new file mode 100755 index 0000000..ce369b6 --- /dev/null +++ b/etc/levee_test/s01-simulation_sample_param.sh @@ -0,0 +1,522 @@ +#!/bin/sh +#========================================================== +# CaMa-Flood sample go script: levee scheme test +# -- Multi 1-year simulations (2000 spinup -> 2000 -> 2001) +# -- Daily runoff forcing (plain binary) at 1deg resolution +# +# (C) D.Yamazaki & E. Dutra (U-Tokyo/FCUL) Aug 2019 +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# You may not use this file except in compliance with the License. +# You may obtain a copy of the License at: http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software distributed under the License is +# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and limitations under the License. +#========================================================== + +#*** PBS setting when needed +#PBS -q E20 +#PBS -l select=1:ncpus=20:mem=10gb +#PBS -j oe +#PBS -m ea +#PBS -V + +#================================================ +# (0) Basic Setting (for workstation) + +#*** 0a. Set CaMa-Flood base directory +BASE=`pwd`/../../ +# BASE="/home/yamadai/work/CaMa_v404/cmf_v404_pkg" # setting for PBS in cluster + +echo $BASE + +#*** 0b. Set dynamic library if needed +export IFORTLIB="/opt/intel/lib:/opt/intel/mkl/lib" +export HDF5LIB="/opt/local/hdf5-1.10.5/lib" +export DYLD_LIBRARY_PATH="${HDF5LIB}:${IFORTLIB}:${DYLD_LIBRARY_PATH}" + +#*** 0c. OpenMP thread number +export OMP_NUM_THREADS=16 # OpenMP cpu num + +#================================================ +# (1) Experiment setting +# -- some non-default options can be modified in NAMELIST section + +#============================ +#*** 1a. Experiment directory setting +EXP="levee_sample" # experiment name (output directory name) +RDIR=${BASE}/out/${EXP} # directory to run CaMa-Flood +EXE="MAIN_cmf" # Execute file name +PROG=${BASE}/src/${EXE} # location of Fortran main program +NMLIST="./input_cmf.nam" # standard namelist +LOGOUT="./log_CaMa.txt" # standard log output + + +#============================ +#*** 1b. Model physics option +DT=86400 # base DT (modified in physics loop by LADPSTP) +LADPSTP=".TRUE." # .TRUE. for adaptive time step + +LFPLAIN=".TRUE." # .TRUE. to activate floodplain storage +LKINE=".FALSE." # .TRUE. to use kinematic wave equation +LFLDOUT=".TRUE." # .TRUE. to activate floodplain discharge +LPTHOUT=".TRUE." # .TRUE. to activate bifurcation flow, mainly for delta simulation +LDAMOUT=".FALSE." # .TRUE. to activate reservoir operation (under development) + +#*** opt. Levee Scheme +LLEVEE=".TRUE." # .TRUE. to activate levee scheme + +#============================ +#*** 1c. simulation time +YSTA=2000 # start year ( from YSTA / Jan 1st _ 00:00) +YEND=2001 # end year (until YEND / Dec 31st _ 24:00) +SPINUP=0 # [0]: zero-storage start, [1]: from restart file +NSP=1 # spinup repeat time + + +#============================ +#*** 1d. spinup setting + +#* input restart file +LRESTART="" # see (3) set each year # TRUE. to use restart initial condition +CRESTSTO="" # see (3) set each year # input restart FIle +LSTOONLY=".FALSE." # .TRUE. for storage only restart (for assimilation) + +#* output restart file +CRESTDIR="./" # output restart file directory +CVNREST="restart" # output restart file prefix +LRESTCDF=".FALSE." # .TRUE. to use netCDF restart file +LRESTDBL=".FALSE." # .TRUE. for binary restart double precision +IFRQ_RST="0" # output restat frequency. + # [0]: only at last time, [1,2,3,...,24] hourly restart, [30]: monthly restart + + +#============================ +#*** 1e. forcing setting +IFRQ_INP="24" # input forcing frequency: [1,2,3,...,24] hour +DROFUNIT="86400000" # [mm/day->m/s] # runoff unit conversion + +#----- for plain binary runoff forcing +LINPCDF=".FALSE." # true for netCDF runoff +LINTERP=".TRUE." # .TRUE. to interporlate with input matrix +LINTERPCDF=".FALSE." # .TRUE. to use netCDF input matrix +CROFDIR="${BASE}/inp/test_1deg/runoff/" # runoff directory +CROFPRE="Roff____" # runoff prefix/suffix +CROFSUF=".one" # $(CROFPRE)YYYYMMDD$(CROFSUF) + +###** sub-surface runoff scheme (not available with plain binary runoff) +LROSPLIT=".FALSE." # .TRUE. for sub-surface runoff +###CSUBDIR="NONE" # sub-surface runoff directory +###CSUBPRE="NONE" # sub-surface runoff prefix/suffix +###CSUBSUF="NONE" # $(PREFIX)YYYYMMDD$(SUFFIX) + +#----- for netCDF runoff forcing ### +###LINPCDF=".TRUE." # true for netCDF runoff +###LINTERP=".TRUE." # .TRUE. to interporlate with input matrix +###LINTERPCDF=".FALSE." # .TRUE. to use netCDF input matrix +###CROFDIR="${BASE}/inp/test_15min_nc/" # runoff directory +###CROFPRE="e2o_ecmwf_wrr2_glob15_day_Runoff_" # runoff prefix/suffix +###CROFCDF="" # see (3) set each year # netCDF runoff file +###CVNROF="Runoff" # netCDF runoff variable name +###CVNSUB="" # netCDF runoffsub variable name +###SYEARIN="" # see (3) set each year # netCDF runoff file, start date +###SMONIN="" # see (3) set each year +###SDAYIN="" # see (3) set each year +###SHOURIN="" # see (3) set each year + + +#============================ +#*** 1f. river map & topography +FMAP="${BASE}/map/glb_15min" # map directory +CDIMINFO="${FMAP}/diminfo_test-1deg.txt" # dimention information file +CINPMAT=${FMAP}/inpmat_test-1deg.bin # runoff input matrix for interporlation +#CDIMINFO="${FMAP}/diminfo_test-15min_nc.txt" # dimention information file +#CINPMAT=${FMAP}/inpmat_test-15min_nc.bin # runoff input matrix for interporlation + +#----- for plain binary map input +#** basic topography +LMAPCDF=".FALSE." # .TRUE. for netCDF map +CNEXTXY="${FMAP}/nextxy.bin" # downstream xy (river network map) +CGRAREA="${FMAP}/ctmare.bin" # unit-catchment area [m2] +CELEVTN="${FMAP}/elevtn.bin" # channel top elevation [m] +CNXTDST="${FMAP}/nxtdst.bin" # downstream distance [m] +CRIVLEN="${FMAP}/rivlen.bin" # channel length [m] +CFLDHGT="${FMAP}/fldhgt.bin" # floodplain elevation profile (height above 'elevtn') [m] + +#** channel parameter +###CRIVWTH=${FMAP}/rivwth.bin" # channel width [m] (empirical power-low) +CRIVWTH="${FMAP}/rivwth_gwdlr.bin" # channel width [m] (GWD-LR + filled with empirical) +CRIVHGT="${FMAP}/rivhgt.bin" # channel depth [m] (empirical power-low) +CRIVMAN="${FMAP}/rivman.bin" # manning coefficient river (The one in flood plain is a global parameter; set $PMANFLD below.) + +CLEVFRC="${FMAP}/lev_frc_global.bin" # levee unprotected fraction [0-1] +CLEVHGT="${FMAP}/lev_hgt_global_vic_up4_new.bin" # levee height [m] + + + +#** bifurcation channel info +CPTHOUT="${FMAP}/bifprm.txt" # bifurcation channel list + +###** groundwater delay (not available in plain binary runoff/map) +LGDWDLY=".FALSE." # .TRUE. to actuvate groundwater delay +#CGDWDLY="" # ground water delay map + +###** mean sea level +LMEANSL=".FALSE." # .TRUE. to use mean sea level data +#CMEANSL="" # mean sea level map + +#----- for netCDF map input +###LMAPCDF=".TRUE." # .TRUE. for netCDF map +###CRIVCLINC="" # netCDF topography map +###CRIVPARNC="" # netCDF river parameters +###CMEANSLNC="" # netCDF mean sea level + + +#============================ +#*** 1g. Dynamic Boundary Sea Level (not default) +LSEALEV=".FALSE." # .TRUE. to activate dynamic sea level +###LSEALEVCDF=".FALSE." +###CSEALEVDIR="NONE" # Sea level boundary DIRECTORY +###CSEALEVPRE="NONE" # Sea level boundary PREFIX +###CSEALEVSUF="NONE" # Sea level boundary SUFFIX +###CSEALEVCDF="NONE" # * Sea level netCDF file name +###CVNSEALEV="sealev" # * Sea Level netCDF variable name +###SYEARSL=1 # * netCDF sea level start year +###SMONSL=1 # * netCDF sea level start year +###SDAYSL=1 # * netCDF sea level start year +###SHOURSL=0 # * netCDF sea level start year +###NSTATIONS=1 # sea level data points +###CSLMAP="NONE" # sea level sta->XY conversion table + +#============================ +#*** 1h. Output Settings +LOUTPUT=".TRUE." # .TRUE. to use CaMa-Flood standard output +IFRQ_OUT=24 # output frequency: [1,2,3,...,24] hour + +LOUTCDF=".FALSE." # .TRUE. netCDF output, .FALSE. plain binary output +COUTDIR="./" # output directory +#CVARSOUT="outflw,storge,fldfrc,maxdph,flddph" # list output variable (comma separated) +CVARSOUT="rivout,rivsto,rivdph,rivvel,fldout,fldsto,flddph,fldfrc,fldare,sfcelv,outflw,storge,pthflw,pthout,maxsto,maxflw,maxdph,levsto,levdph" # list output variable (comma separated) +COUTTAG="" # see (3) set each year # output tag $(COUTDIR)/$(VARNAME)$(OUTTAG).bin + +##### Model Parameters ################ +PMANRIV="0.03D0" # manning coefficient river +PMANFLD="0.10D0" # manning coefficient floodplain +PCADP="0.7" # satety coefficient for CFL condition +PDSTMTH="10000.D0" # downstream distance at river mouth [m] + +#================================================ +# (2) Initial setting + +#*** 2a. create running dir +mkdir -p ${RDIR} +cd ${RDIR} + +#*** 2b. for new simulation, remove old files in running directory + +if [ ${SPINUP} -eq 0 ]; then + rm -rf ${RDIR}/????-sp* + rm -rf ${RDIR}/*.bin + rm -rf ${RDIR}/*.pth + rm -rf ${RDIR}/*.vec + rm -rf ${RDIR}/*.nc + rm -rf ${RDIR}/*.log + rm -rf ${RDIR}/*.txt + rm -rf ${RDIR}/restart* +else + NSP=0 # restart, no spinup +fi + + + +#================================================ +# (3) For each simulation year, modify setting +#-- loop 1-year simulation from $YSTART to $YEND + +ISP=1 ## spinup count +IYR=${YSTA} ## curent year +while [ ${IYR} -le ${YEND} ]; +do + CYR=`printf %04d ${IYR}` ## update file name, bugfix in v3.96a + + #*** 3a. modify restart setting + if [ ${SPINUP} -eq 0 ];then + LRESTART=".FALSE." ## from zero storage + CRESTSTO="" + else + LRESTART=".TRUE." + CRESTSTO="${CVNREST}${CYR}010100.bin" ## from restart file + if [ ${LRESTCDF} = ".TRUE." ]; then + CRESTSTO="${CVNREST}${CYR}010100.nc" ## from restart file + fi + fi + + #*** 3b. update start-end year + SYEAR=$IYR + SMON=1 + SDAY=1 + SHOUR=0 + + EYEAR=`expr $SYEAR + 1` + EMON=1 + EDAY=1 + EHOUR=0 + + ln -sf $PROG $EXE + + #*** 3c. update input / output file data + CSYEAR=`printf %04d ${SYEAR}` + COUTTAG=${CSYEAR} # output file tag + + #CROFCDF="${CROFDIR}/${CROFPRE}${CSYEAR}.nc" # input netCDF runoff file + #SYEARIN=$IYR + #SMONIN=1 + #SDAYIN=1 + #SHOURIN=0 + + + +#================================================ +# (4) Create NAMELIST for simulation year +# it is OK to remove optional variables (set to default in CaMa-Flood) + +rm -f ${NMLIST} + +#*** 0. config +cat >> ${NMLIST} << EOF +&NRUNVER +LADPSTP = ${LADPSTP} ! true: use adaptive time step +LFPLAIN = ${LFPLAIN} ! true: consider floodplain (false: only river channel) +LKINE = ${LKINE} ! true: use kinematic wave +LFLDOUT = ${LFLDOUT} ! true: floodplain flow (high-water channel flow) active +LPTHOUT = ${LPTHOUT} ! true: activate bifurcation scheme +LDAMOUT = ${LDAMOUT} ! true: activate dam operation (under development) +LLEVEE = ${LLEVEE} ! true: activate levee scheme (under development) +LROSPLIT = ${LROSPLIT} ! true: input if surface (Qs) and sub-surface (Qsb) runoff +LGDWDLY = ${LGDWDLY} ! true: Activate ground water reservoir and delay +LSLPMIX = .FALSE. ! true: activate mixed kinematic and local inertia based on slope +LMEANSL = ${LMEANSL} ! true: boundary condition for mean sea level +LSEALEV = ${LSEALEV} ! true: boundary condition for variable sea level +LRESTART = ${LRESTART} ! true: initial condition from restart file +LSTOONLY = ${LSTOONLY} ! true: storage only restart (mainly for data assimilation) +LOUTPUT = ${LOUTPUT} ! true: use standard output (to file) +LGRIDMAP = .TRUE. ! true: for standard XY gridded 2D map +LLEAPYR = .TRUE. ! true: neglect leap year (Feb29 skipped) +LMAPEND = .FALSE. ! true: for map data endian conversion +LSTG_ES = .FALSE. ! true: for Vector Processor optimization (CMF_OPT_FLDSTG_ES) +/ +&NDIMTIME +CDIMINFO = "${CDIMINFO}" ! text file for dimention information +DT = ${DT} ! time step length (sec) +IFRQ_INP = ${IFRQ_INP} ! input forcing update frequency (hour) +/ +&NPARAM +PMANRIV = ${PMANRIV} ! manning coefficient river +PMANFLD = ${PMANFLD} ! manning coefficient floodplain +PGRV = 9.8D0 ! gravity accerelation +PDSTMTH = ${PDSTMTH} ! downstream distance at river mouth [m] +PCADP = ${PCADP} ! CFL coefficient +PMINSLP = 1.D-5 ! minimum slope (kinematic wave) +IMIS = -9999 ! missing value for integer +RMIS = 1.E20 ! missing value for real*4 +DMIS = 1.E20 ! missing value for real*8 +CSUFBIN = '.bin' ! file suffix for plain binary 2D map +CSUFVEC = '.vec' ! file suffix for plain binary 1D vector +CSUFPTH = '.pth' ! file suffix for plain binary bifurcation channel +CSUFCDF = '.nc' ! file suffix for netCDF +/ +EOF + +#*** 1. time +cat >> ${NMLIST} << EOF +&NSIMTIME +SYEAR = ${SYEAR} ! start year +SMON = ${SMON} ! month +SDAY = ${SDAY} ! day +SHOUR = ${SHOUR} ! houe +EYEAR = ${EYEAR} ! end year +EMON = ${EMON} ! month +EDAY = ${EDAY} ! day +EHOUR = ${EHOUR} ! hour +/ +EOF + +#*** 2. map +cat >> ${NMLIST} << EOF +&NMAP +LMAPCDF = ${LMAPCDF} ! * true for netCDF map input +CNEXTXY = "${CNEXTXY}" ! river network nextxy +CGRAREA = "${CGRAREA}" ! catchment area +CELEVTN = "${CELEVTN}" ! bank top elevation +CNXTDST = "${CNXTDST}" ! distance to next outlet +CRIVLEN = "${CRIVLEN}" ! river channel length +CFLDHGT = "${CFLDHGT}" ! floodplain elevation profile +CRIVWTH = "${CRIVWTH}" ! channel width +CRIVHGT = "${CRIVHGT}" ! channel depth +CRIVMAN = "${CRIVMAN}" ! river manning coefficient +CPTHOUT = "${CPTHOUT}" ! bifurcation channel table +CGDWDLY = "${CGDWDLY}" ! Groundwater Delay Parameter +CMEANSL = "${CMEANSL}" ! mean sea level +CRIVCLINC = "${CRIVCLINC}" ! * river map netcdf +CRIVPARNC = "${CRIVPARNC}" ! * river parameter netcdf (width, height, manning, ground water delay) +CMEANSLNC = "${CMEANSLNC}" ! * mean sea level netCDF +/ +EOF + +#*** 3. restart +cat >> ${NMLIST} << EOF +&NRESTART +CRESTSTO = "${CRESTSTO}" ! restart file +CRESTDIR = "${CRESTDIR}" ! restart directory +CVNREST = "${CVNREST}" ! restart variable name +LRESTCDF = ${LRESTCDF} ! * true for netCDF restart file (double precision) +LRESTDBL = ${LRESTDBL} ! * true for double precision binary restart +IFRQ_RST = ${IFRQ_RST} ! restart write frequency (1-24: hour, 0:end of run) +/ +EOF + +#*** 4. forcing +if [ ${LINPCDF} = ".FALSE." ]; then +cat >> ${NMLIST} << EOF +&NFORCE +LINPCDF = ${LINPCDF} ! true for netCDF runoff +LINTERP = ${LINTERP} ! true for runoff interpolation using input matrix +LINPEND = .FALSE. ! true for runoff endian conversion +CINPMAT = "${CINPMAT}" ! input matrix file name +DROFUNIT = ${DROFUNIT} ! runoff unit conversion +CROFDIR = "${CROFDIR}" ! runoff input directory +CROFPRE = "${CROFPRE}" ! runoff input prefix +CROFSUF = "${CROFSUF}" ! runoff input suffix +CSUBDIR = "${CSUBDIR}" ! sub-surface runoff input directory +CSUBPRE = "${CSUBPRE}" ! sub-surface runoff input prefix +CSUBSUF = "${CSUBSUF}" ! sub-surface runoff input suffix +/ +EOF + +elif [ ${LINPCDF} = ".TRUE." ]; then +cat >> ${NMLIST} << EOF +&NFORCE +LINPCDF = ${LINPCDF} ! true for netCDF runoff +LINTERP = ${LINTERP} ! true for runoff interpolation using input matrix +LINPEND = .FALSE. ! true for runoff endian conversion +LITRPCDF = ${LINTERPCDF} ! * true for netCDF input matrix +CINPMAT = "${CINPMAT}" ! input matrix file name +DROFUNIT = ${DROFUNIT} ! runoff unit conversion +CROFCDF = "${CROFCDF}" ! * netCDF input runoff file name +CVNROF = "${CVNROF}" ! * netCDF input runoff variable name +CVNSUB = "${CVNSUB}" ! * netCDF input runoff variable name +SYEARIN = ${SYEARIN} ! * netCDF input start year +SMONIN = ${SMONIN} ! * netCDF input start year +SDAYIN = ${SDAYIN} ! * netCDF input start year +SHOURIN = ${SHOURIN} ! * netCDF input start year +/ +EOF + +fi # (if LINPCDF) + +#*** 5. outputs +cat >> ${NMLIST} << EOF +&NOUTPUT +COUTDIR = "${COUTDIR}" ! OUTPUT DIRECTORY +CVARSOUT = "${CVARSOUT}" ! Comma-separated list of output variables to save +COUTTAG = "${COUTTAG}" ! Output Tag Name for each experiment +LOUTVEC = .FALSE ! TRUE FOR VECTORIAL OUTPUT, FALSE FOR NX,NY OUTPUT +LOUTCDF = ${LOUTCDF} ! * true for netcdf outptu false for binary +NDLEVEL = 0 ! * NETCDF DEFLATION LEVEL +IFRQ_OUT = ${IFRQ_OUT} ! output data write frequency (hour) +/ +EOF + +#### 6. sea level (optional) +#cat >> ${NMLIST} << EOF +#&NBOUND +#LSEALEVCDF = ${LSEALEVCDF} ! * true : netCDF sea level boundary +#CSEALEVDIR = "${CSEALEVDIR}" ! Sea level boundary DIRECTORY +#CSEALEVPRE = "${CSEALEVPRE}" ! Sea level boundary PREFIX +#CSEALEVSUF = "${CSEALEVSUF}" ! Sea level boundary SUFFIX +#CSEALEVCDF = "${CSEALEVCDF}" ! * Sea level netCDF file name +#CVNSEALEV = "${CVNSEALEV}" ! * Sea Level netCDF variable name +#SYEARSL = ${SYEARSL} ! * netCDF sea level start year +#SMONSL = ${SMONSL} ! * netCDF sea level start year +#SDAYSL = ${SDAYSL} ! * netCDF sea level start year +#SHOURSL = ${SHOURSL} ! * netCDF sea level start year +#NSTATIONS = ${NSTATIONS} ! sea level data points +#CSLMAP = "${CSLMAP} ! station to XY conversion table +#IFRQ_SL = ${IFRQ_SL} ! sea level boundary update frequency (min) +#/ +#EOF + +#*** 5. outputs +cat >> ${NMLIST} << EOF +&NLEVEE +CLEVFRC = "${CLEVFRC}" ! Levee unprotectdd fraction +CLEVHGT = "${CLEVHGT}" ! Levee Height +/ +EOF + +#================================================ +# (5) Execute main program + +echo "start: ${SYEAR}" `date` >> log.txt +time ./${EXE} >> log.txt +echo "end: ${SYEAR}" `date` >> log.txt + +mv ${LOGOUT} log_CaMa-${CYR}.txt + +#================================================ +# (6) manage spin up + +# if curent spinup time $ISP < required spinup time $NSP +# copy the restart file restart$(IYR+1) to restart${IYR} +# copy the outputs to directory "${IYR}-sp1" + +SPINUP=1 +if [ ${IYR} -eq ${YSTA} ]; +then + if [ ${ISP} -le ${NSP} ]; + then + IYR1=`expr ${IYR} + 1` + CYR1=`printf %04d ${IYR1}` + cp -f ${CVNREST}${CYR1}010100.bin ${CVNREST}${CYR}010100.bin-sp${ISP} 2> /dev/null + mv -f ${CVNREST}${CYR1}010100.bin ${CVNREST}${CYR}010100.bin 2> /dev/null + + cp -f ${CVNREST}${CYR1}010100.bin.pth ${CVNREST}${CYR}010100.bin.pth-sp${ISP} 2> /dev/null + mv -f ${CVNREST}${CYR1}010100.bin.pth ${CVNREST}${CYR}010100.bin.pth 2> /dev/null + + cp -f ${CVNREST}${CYR1}010100.nc ${CVNREST}${CYR}010100.nc-sp${ISP} 2> /dev/null + mv -f ${CVNREST}${CYR1}010100.nc ${CVNREST}${CYR}010100.nc 2> /dev/null + + mkdir -p ${CYR}-sp${ISP} + mv -f ./${CVNREST}${CYR}010100.bin-sp${ISP} ${CYR}-sp${ISP} 2> /dev/null + mv -f ./${CVNREST}${CYR}010100.bin.pth-sp${ISP} ${CYR}-sp${ISP} 2> /dev/null + mv -f ./${CVNREST}${CYR}010100.nc-sp${ISP} ${CYR}-sp${ISP} 2> /dev/null + mv -f ./*${CYR}.bin ${CYR}-sp${ISP} 2> /dev/null + mv -f ./*${CYR}.pth ${CYR}-sp${ISP} 2> /dev/null + mv -f ./o_*${CYR}.nc ${CYR}-sp${ISP} 2> /dev/null + mv -f ./*${CYR}.log ${CYR}-sp${ISP} 2> /dev/null + mv -f ./log_CaMa-${CYR}.txt ${CYR}-sp${ISP} 2> /dev/null + + ISP=`expr ${ISP} + 1` + else + ISP=0 + IYR=`expr ${IYR} + 1` + fi +else + IYR=`expr ${IYR} + 1` +fi + +#================================================ +# (7) End of each year loop. Back to (3) + +done # loop to next year simulation + + +# Simulation without levee +echo "please execute gosh/test1-glb_15min.sh as reference simulation" + + +exit 0 + + + diff --git a/etc/levee_test/s03-analysis.sh b/etc/levee_test/s02-analysis.sh similarity index 81% rename from etc/levee_test/s03-analysis.sh rename to etc/levee_test/s02-analysis.sh index 558fbf2..5c22dc2 100755 --- a/etc/levee_test/s03-analysis.sh +++ b/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 diff --git a/etc/levee_test/s02-simulation.sh b/etc/levee_test/s02-simulation.sh deleted file mode 100755 index 8b4368e..0000000 --- a/etc/levee_test/s02-simulation.sh +++ /dev/null @@ -1,8 +0,0 @@ -#!/bin/sh - -# Execute simulaton with levee -./test_levee.sh - -# Simulation without levee -echo "please execute gosh/test1-glb_15min.sh as reference simulation" - diff --git a/etc/levee_test/src/define_protect_pix.F90 b/etc/levee_test/src/define_protect_pix.F90 index bb5c092..07308bc 100755 --- a/etc/levee_test/src/define_protect_pix.F90 +++ b/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_endeast2 .or.lat_orinorth2 ) 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 diff --git a/etc/levee_test/src/define_protect_pix_old.F90 b/etc/levee_test/src/define_protect_pix_old.F90 new file mode 100755 index 0000000..bb5c092 --- /dev/null +++ b/etc/levee_test/src/define_protect_pix_old.F90 @@ -0,0 +1,251 @@ + program downscale_flddph +! =============================================== + implicit none +! CaMa-Flood parameters + character*256 :: param !! river map parameters + integer :: iXX, iYY + integer :: nXX, nYY !! grid number (river network map) + integer :: nflp !! floodplain layers + real*8 :: gsize !! grid size [deg] + real*8 :: west, east, north, south !! domain (river network map) + +! HydroSHEDS parameters !! (from location.txt) + integer :: i, narea !! area ID + character*256 :: area !! area code + integer :: ix, iy, jx, jy + integer :: nx, ny !! grid number (hires data) + real*8 :: csize !! size of pixel [deg] + real*8 :: lon_ori !! west edge + real*8 :: lon_end !! east edge + real*8 :: lat_ori !! north edge + real*8 :: lat_end !! south edge + + real*8 :: west2, east2, north2, south2 !! output domain + integer :: mx, my + + character*256 :: list_loc +! + integer*4,allocatable :: nextXX(:,:) !! downstream (jXX,jYY) + integer*2,allocatable :: catmXX(:,:), catmYY(:,:) !! catchment (iXX,iYY) of pixel (ix,iy) + 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 :: protect(:,:) !! downscaled flood depth [m] + real,allocatable :: flddif2(:,:),hand2(:,:) !! downscaled flood depth [m] + +! + 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 :: buf + integer :: ios +! =============================================== +! downscale target domain + call getarg(1,buf) + read(buf,*) west2 + call getarg(2,buf) + read(buf,*) east2 + call getarg(3,buf) + read(buf,*) south2 + call getarg(4,buf) + read(buf,*) north2 + + call getarg(5,hires) !! downscale resolution + + param=trim(mapdir)//'params.txt' + + open(11,file=param,form='formatted') + read(11,*) nXX + read(11,*) nYY + read(11,*) nflp + read(11,*) gsize + read(11,*) west + read(11,*) east + read(11,*) south + read(11,*) north + close(11) + +! CaMa-Flood simulation domain + east =west +real(nXX)*gsize + south=north-real(nXX)*gsize + + if( trim(hires)=='1sec' )then + csize=1./3600. + elseif( trim(hires)=='3sec' )then + csize=1./1200. + elseif( trim(hires)=='5sec' )then + csize=1./720. + elseif( trim(hires)=='15sec' )then + csize=1./240. + elseif( trim(hires)=='30sec' )then + csize=1./120. + elseif( trim(hires)=='1min' )then + csize=1./60. + else + stop + endif + +! calculate number of cell size + mx=nint( (east2 -west2 )/csize ) !! downscale domain x*y + my=nint( (north2-south2)/csize ) + + if( mx*my>525000000 )then + print *, 'downscale domain too large: mx*my*4>integer limit' + stop + endif + + print *, 'domain:', west2, east2, south2, north2, mx, my + +! ========== + + allocate(nextXX(nXX,nYY),levbas(nXX,nYY),levhgt(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 + close(12) + + +! open hires files + list_loc=trim(mapdir)//trim(hires)//'/location.txt' + open(11,file=list_loc,form='formatted') + read(11,*) narea + read(11,*) + + do i=1, narea + read(11,*) buf, area, lon_ori, lon_end, lat_end, lat_ori, nx, ny, csize + + if( lon_endeast2 .or.lat_orinorth2 ) goto 1090 !! out of domain + allocate(catmXX(nx,ny),catmYY(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' + print *, rfile + open(21,file=rfile,form='unformatted',access='direct',recl=2*nx*ny,status='old',iostat=ios) + if( ios==0 )then + read(21,rec=1) catmXX + read(21,rec=2) catmYY + 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 + read(21,rec=1) flddif + close(21) + else + print *, '*******************' + print *, 'no data: ', rfile + stop + endif + + + rfile=trim(mapdir)//trim(hires)//'/'//trim(area)//'.rivwth.bin' + open(21,file=rfile,form='unformatted',access='direct',recl=4*nx*ny,status='old',iostat=ios) + if( ios==0 )then + read(21,rec=1) rivwth + close(21) + else + print *, '*******************' + print *, 'no data: ', rfile + stop + endif + + rfile=trim(mapdir)//trim(hires)//'/'//trim(area)//'.hand.bin' + open(21,file=rfile,form='unformatted',access='direct',recl=4*nx*ny,status='old',iostat=ios) + if( ios==0 )then + read(21,rec=1) hand + close(21) + else + print *, '*******************' + print *, 'no data: ', rfile + stop + endif + + do ix=1, nx + lon(ix)=lon_ori+(real(ix)-0.5)*csize + if( lon(ix)>=180. ) lon(ix)=lon(ix)-360. + if( lon(ix)<-180. ) lon(ix)=lon(ix)+360. + end do + do iy=1, ny + lat(iy) =lat_ori-(real(iy)-0.5)*csize + end do + + do iy=1, ny + do ix=1, nx + if( lon(ix)>west2 .and. lon(ix)south2 .and. lat(iy)0 )then + hand2(jx,jy) =hand(ix,iy)**0.3 + hand2(jx,jy) =min(hand2(jx,jy),5.) + else + hand2(jx,jy)=-1 + endif + + if( catmXX(ix,iy)>0 )then + 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 + 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 + if( protect(jx,jy)==0 ) protect(jx,jy)=1 + endif + elseif( catmXX(ix,iy)/=-9999 )then + protect(jx,jy)=-1 + endif + endif + end do + end do + + deallocate(catmXX,catmYY,flddif,rivwth,hand) + deallocate(lon,lat) + + 1090 continue + enddo + + hand2(1,1)= 5 + hand2(1,2)=-1 + + fflood='./data/protect_pix.bin' + open(11, file=fflood, form='unformatted', access='direct', recl=4*mx*my) + write(11,rec=1) protect + write(11,rec=2) flddif2 + close(11) + + fflood='./data/hand.bin' + open(11, file=fflood, form='unformatted', access='direct', recl=4*mx*my) + write(11,rec=1) hand2 + close(11) + + end program downscale_flddph + diff --git a/etc/levee_test/s01-set_levee_param.sh b/etc/levee_test/x01-levee_param_simple.sh similarity index 100% rename from etc/levee_test/s01-set_levee_param.sh rename to etc/levee_test/x01-levee_param_simple.sh diff --git a/etc/levee_test/test_levee.sh b/etc/levee_test/x02-simulation_simple.sh similarity index 99% rename from etc/levee_test/test_levee.sh rename to etc/levee_test/x02-simulation_simple.sh index e4b2273..1386b70 100755 --- a/etc/levee_test/test_levee.sh +++ b/etc/levee_test/x02-simulation_simple.sh @@ -511,4 +511,12 @@ fi done # loop to next year simulation + +# Simulation without levee +echo "please execute gosh/test1-glb_15min.sh as reference simulation" + + exit 0 + + + diff --git a/etc/levee_test/t03-downscale.sh b/etc/levee_test/x03-downscale.sh similarity index 100% rename from etc/levee_test/t03-downscale.sh rename to etc/levee_test/x03-downscale.sh