diff --git a/driver_scripts/driver_grid.cray.sh b/driver_scripts/driver_grid.cray.sh index 2abe26f72..28e68c4c7 100755 --- a/driver_scripts/driver_grid.cray.sh +++ b/driver_scripts/driver_grid.cray.sh @@ -36,21 +36,23 @@ # "regional_gfdl" - stand-alone gfdl regional grid # "regional_esg" - stand-alone extended Schmidt gnomonic # (esg) regional grid -# 3) For "stretch" and "nest" grids, set the stretching factor - +# 3) For "uniform" grids - to include lake fraction and +# depth, set "add_lake" to true, and the "lake_cutoff" value. +# 4) For "stretch" and "nest" grids, set the stretching factor - # "stretch_fac", and center lat/lon of highest resolution # tile - "target_lat" and "target_lon". -# 4) For "nest" grids, set the refinement ratio - "refine_ratio", +# 5) For "nest" grids, set the refinement ratio - "refine_ratio", # the starting/ending i/j index location within the parent # tile - "istart_nest", "jstart_nest", "iend_nest", "jend_nest" -# 5) For "regional_gfdl" grids, set the "halo". Default is three +# 6) For "regional_gfdl" grids, set the "halo". Default is three # rows/columns. -# 6) For "regional_esg" grids, set center lat/lon of grid, +# 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 7) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TMPDIR - and path to the repository # clone - home_dir. -# 8) Submit script: "cat $script | bsub". -# 9) All files will be placed in "out_dir". +# 9) Submit script: "cat $script | bsub". +# 10) All files will be placed in "out_dir". # #----------------------------------------------------------------------- @@ -62,11 +64,14 @@ module list # Set grid specs here. #----------------------------------------------------------------------- -export gtype=uniform # 'uniform', 'stretch', 'nest' +export gtype=uniform # 'uniform', 'stretch', 'nest' # 'regional_gfdl', 'regional_esg' if [ $gtype = uniform ]; then export res=96 + export add_lake=false # Add lake frac and depth to orography data. + # Uniform grids only. + export lake_cutoff=0.20 # lake frac less than lake_cutoff is ignored elif [ $gtype = stretch ]; then export res=96 export stretch_fac=1.5 # Stretching factor for the grid @@ -109,6 +114,10 @@ export home_dir=$LS_SUBCWD/.. export TMPDIR=/gpfs/hps3/stmp/$LOGNAME/fv3_grid.$gtype export out_dir=/gpfs/hps3/stmp/$LOGNAME/my_grids +#----------------------------------------------------------------------- +# Should not need to change anything below here. +#----------------------------------------------------------------------- + export NODES=1 export APRUN="aprun -n 1 -N 1 -j 1 -d 1 -cc depth" export APRUN_SFC="aprun -j 1 -n 24 -N 24" @@ -117,6 +126,7 @@ export OMP_NUM_THREADS=6 export OMP_STACKSIZE=2048m export KMP_AFFINITY=disabled export machine=WCOSS_C +export NCDUMP=/gpfs/hps/usrx/local/prod/NetCDF/4.2/intel/sandybridge/bin/ncdump ulimit -a ulimit -s unlimited diff --git a/driver_scripts/driver_grid.dell.sh b/driver_scripts/driver_grid.dell.sh index a1d891fe7..c6da158fc 100755 --- a/driver_scripts/driver_grid.dell.sh +++ b/driver_scripts/driver_grid.dell.sh @@ -38,21 +38,23 @@ # "regional_gfdl" - stand-alone gfdl regional grid # "regional_esg" - stand-alone extended Schmidt gnomonic # (esg) regional grid -# 3) For "stretch" and "nest" grids, set the stretching factor - +# 3) For "uniform" grids - to include lake fraction and +# depth, set "add_lake" to true, and the "lake_cutoff" value. +# 4) For "stretch" and "nest" grids, set the stretching factor - # "stretch_fac", and center lat/lon of highest resolution # tile - "target_lat" and "target_lon". -# 4) For "nest" grids, set the refinement ratio - "refine_ratio", +# 5) For "nest" grids, set the refinement ratio - "refine_ratio", # the starting/ending i/j index location within the parent # tile - "istart_nest", "jstart_nest", "iend_nest", "jend_nest" -# 5) For "regional_gfdl" grids, set the "halo". Default is three +# 6) For "regional_gfdl" grids, set the "halo". Default is three # rows/columns. -# 6) For "regional_esg" grids, set center lat/lon of grid, +# 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 7) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TMPDIR - and path to the repository # clone - home_dir. -# 8) Submit script: "cat $script | bsub". -# 9) All files will be placed in "out_dir". +# 9) Submit script: "cat $script | bsub". +# 10) All files will be placed in "out_dir". # #----------------------------------------------------------------------- @@ -65,10 +67,13 @@ module list #----------------------------------------------------------------------- export gtype=uniform # 'uniform', 'stretch', 'nest', - # 'regional_gfdl', 'regional_esg' + # 'regional_gfdl', 'regional_esg' if [ $gtype = uniform ]; then export res=96 + export add_lake=false # Add lake frac and depth to orography data. + # Uniform grids only. + export lake_cutoff=0.20 # lake frac less than lake_cutoff is ignored elif [ $gtype = stretch ]; then export res=96 export stretch_fac=1.5 # Stretching factor for the grid diff --git a/driver_scripts/driver_grid.hera.sh b/driver_scripts/driver_grid.hera.sh index 50f13fbee..75c986aba 100755 --- a/driver_scripts/driver_grid.hera.sh +++ b/driver_scripts/driver_grid.hera.sh @@ -36,21 +36,23 @@ # "regional_gfdl" - stand-alone gfdl regional grid # "regional_esg" - stand-alone extended Schmidt gnomonic # (esg) regional grid -# 3) For "stretch" and "nest" grids, set the stretching factor - +# 3) For "uniform" grids - to include lake fraction and +# depth, set "add_lake" to true, and the "lake_cutoff" value. +# 4) For "stretch" and "nest" grids, set the stretching factor - # "stretch_fac", and center lat/lon of highest resolution # tile - "target_lat" and "target_lon". -# 4) For "nest" grids, set the refinement ratio - "refine_ratio", +# 5) For "nest" grids, set the refinement ratio - "refine_ratio", # the starting/ending i/j index location within the parent # tile - "istart_nest", "jstart_nest", "iend_nest", "jend_nest" -# 5) For "regional_gfdl" grids, set the "halo". Default is three +# 6) For "regional_gfdl" grids, set the "halo". Default is three # rows/columns. -# 6) For "regional_esg" grids, set center lat/lon of grid, +# 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 7) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TMPDIR - and path to the repository # clone - home_dir. -# 8) Submit script: "sbatch $script". -# 9) All files will be placed in "out_dir". +# 9) Submit script: "sbatch $script". +# 10) All files will be placed in "out_dir". # #----------------------------------------------------------------------- @@ -69,6 +71,9 @@ export gtype=uniform # 'uniform', 'stretch', 'nest' if [ $gtype = uniform ]; then export res=96 + export add_lake=false # Add lake frac and depth to orography data. + # Uniform grids only. + export lake_cutoff=0.20 # lake frac less than lake_cutoff is ignored elif [ $gtype = stretch ]; then export res=96 export stretch_fac=1.5 # Stretching factor for the grid diff --git a/driver_scripts/driver_grid.jet.sh b/driver_scripts/driver_grid.jet.sh index 179864abd..650b39820 100755 --- a/driver_scripts/driver_grid.jet.sh +++ b/driver_scripts/driver_grid.jet.sh @@ -37,21 +37,23 @@ # "regional_gfdl" - stand-alone gfdl regional grid # "regional_esg" - stand-alone extended Schmidt gnomonic # (esg) regional grid -# 3) For "stretch" and "nest" grids, set the stretching factor - +# 3) For "uniform" grids - to include lake fraction and +# depth, set "add_lake" to true, and the "lake_cutoff" value. +# 4) For "stretch" and "nest" grids, set the stretching factor - # "stretch_fac", and center lat/lon of highest resolution # tile - "target_lat" and "target_lon". -# 4) For "nest" grids, set the refinement ratio - "refine_ratio", +# 5) For "nest" grids, set the refinement ratio - "refine_ratio", # the starting/ending i/j index location within the parent # tile - "istart_nest", "jstart_nest", "iend_nest", "jend_nest" -# 5) For "regional_gfdl" grids, set the "halo". Default is three +# 6) For "regional_gfdl" grids, set the "halo". Default is three # rows/columns. -# 6) For "regional_esg" grids, set center lat/lon of grid, +# 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 7) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TMPDIR - and path to the repository # clone - home_dir. -# 8) Submit script: "sbatch $script". -# 9) All files will be placed in "out_dir". +# 9) Submit script: "sbatch $script". +# 10) All files will be placed in "out_dir". # #----------------------------------------------------------------------- @@ -70,6 +72,9 @@ export gtype=uniform # 'uniform', 'stretch', 'nest' if [ $gtype = uniform ]; then export res=96 + export add_lake=false # Add lake frac and depth to orography data. + # Uniform grids only. + export lake_cutoff=0.20 # lake frac less than lake_cutoff is ignored elif [ $gtype = stretch ]; then export res=96 export stretch_fac=1.5 # Stretching factor for the grid diff --git a/driver_scripts/driver_grid.orion.sh b/driver_scripts/driver_grid.orion.sh index 66a2955d5..208596da3 100755 --- a/driver_scripts/driver_grid.orion.sh +++ b/driver_scripts/driver_grid.orion.sh @@ -36,21 +36,23 @@ # "regional_gfdl" - stand-alone gfdl regional grid # "regional_esg" - stand-alone extended Schmidt # gnomonic (esg) regional grid -# 3) For "stretch" and "nest" grids, set the stretching factor - +# 3) For "uniform" grids - to include lake fraction and +# depth, set "add_lake" to true, and the "lake_cutoff" value. +# 4) For "stretch" and "nest" grids, set the stretching factor - # "stretch_fac", and center lat/lon of highest resolution # tile - "target_lat" and "target_lon". -# 4) For "nest" grids, set the refinement ratio - "refine_ratio", +# 5) For "nest" grids, set the refinement ratio - "refine_ratio", # the starting/ending i/j index location within the parent # tile - "istart_nest", "jstart_nest", "iend_nest", "jend_nest" -# 5) For "regional_gfdl" grids, set the "halo". Default is three +# 6) For "regional_gfdl" grids, set the "halo". Default is three # rows/columns. -# 6) For "regional_esg" grids, set center lat/lon of grid, +# 7) For "regional_esg" grids, set center lat/lon of grid, # - "target_lat/lon" - the i/j dimensions - "i/jdim", the # x/y grid spacing - "delx/y", and halo. -# 7) Set working directory - TMPDIR - and path to the repository +# 8) Set working directory - TMPDIR - and path to the repository # clone - home_dir. -# 8) Submit script: "sbatch $script". -# 9) All files will be placed in "out_dir". +# 9) Submit script: "sbatch $script". +# 10) All files will be placed in "out_dir". # #----------------------------------------------------------------------- @@ -69,6 +71,9 @@ export gtype=uniform # 'uniform', 'stretch', 'nest', if [ $gtype = uniform ]; then export res=96 + export add_lake=false # Add lake frac and depth to orography data. + # Uniform grids only. + export lake_cutoff=0.20 # lake frac less than lake_cutoff is ignored elif [ $gtype = stretch ]; then export res=96 export stretch_fac=1.5 # Stretching factor for the grid diff --git a/sorc/CMakeLists.txt b/sorc/CMakeLists.txt index fbf3ce95c..6daacc548 100644 --- a/sorc/CMakeLists.txt +++ b/sorc/CMakeLists.txt @@ -15,5 +15,7 @@ add_subdirectory(mkgfsnemsioctl.fd) add_subdirectory(fre-nctools.fd) add_subdirectory(grid_tools.fd) add_subdirectory(chgres_cube.fd) -add_subdirectory(orog.fd) +add_subdirectory(orog_mask_tools.fd/orog.fd) +add_subdirectory(orog_mask_tools.fd/lake.fd) +add_subdirectory(orog_mask_tools.fd/inland.fd) add_subdirectory(sfc_climo_gen.fd) diff --git a/sorc/orog_mask_tools.fd/inland.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/inland.fd/CMakeLists.txt new file mode 100644 index 000000000..e2078e76f --- /dev/null +++ b/sorc/orog_mask_tools.fd/inland.fd/CMakeLists.txt @@ -0,0 +1,19 @@ +set(fortran_src + nb.F90 + inland.F90) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -assume byterecl") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU|Clang|AppleClang)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fno-range-check") +endif() + +set(exe_name inland) + +add_executable(${exe_name} ${fortran_src}) + +target_link_libraries( + ${exe_name} + NetCDF::NetCDF_Fortran) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir}) diff --git a/sorc/orog_mask_tools.fd/inland.fd/inland.F90 b/sorc/orog_mask_tools.fd/inland.fd/inland.F90 new file mode 100644 index 000000000..cdeb5a303 --- /dev/null +++ b/sorc/orog_mask_tools.fd/inland.fd/inland.F90 @@ -0,0 +1,225 @@ +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +! This program creates the inland mask and writes it to the orography data files. +! +! Ning Wang, July 1, 2020, original version. +! +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +PROGRAM inland_mask + USE cs_nb + IMPLICIT NONE + + INTEGER :: tile, i, j + TYPE(nb_gp_idx) :: nbs + + REAL, ALLOCATABLE :: inland(:,:,:) + REAL, ALLOCATABLE :: land_frac(:,:,:) + INTEGER :: i_ctr, j_ctr, tile_beg, tile_end + INTEGER :: cs_res + CHARACTER(len=32) :: arg + INTEGER :: stat + INTEGER :: max_rd + REAL :: cutoff + + LOGICAL, ALLOCATABLE :: done(:,:,:) + +! CALL getarg(0, arg) ! get the program name +! IF (iargc() /= 3) THEN +! PRINT*, 'Usage: ', trim(arg), ' [resolution (48,96, ...)] [non-land cutoff] [max recursive depth]' +! STOP +! ENDIF + + CALL getarg(1, arg) + READ(arg,*,iostat=stat) cs_res + CALL getarg(2, arg) + READ(arg,*,iostat=stat) cutoff + CALL getarg(3, arg) + READ(arg,*,iostat=stat) max_rd + + ALLOCATE(done(cs_res,cs_res,6)) + ALLOCATE(inland(cs_res,cs_res,6)) + ALLOCATE(land_frac(cs_res,cs_res,6)) + + tile_beg = 1; tile_end = 6 + +! init inter-panel neighbor index + CALL idx_init(cs_res) + +! read in orography data + CALL read_orog(cs_res) + +! create a inland mask + CALL mark_global_inland(cs_res) + +! write back to the orography data files + CALL write_inland(cs_res) + +CONTAINS + +SUBROUTINE mark_global_inland(cs_res) + INTEGER, INTENT(IN) :: cs_res + + done = .false. + inland = 1.0 + i_ctr = cs_res/2; j_ctr = cs_res/2 + + CALL mark_global_inland_rec_d(i_ctr, j_ctr, 2, 0) + +END SUBROUTINE mark_global_inland + +RECURSIVE SUBROUTINE mark_global_inland_rec_d(i, j, t, rd) + INTEGER, INTENT(IN) :: i, j, t, rd + + TYPE(nb_gp_idx) :: nbs + INTEGER :: k, nrd + + IF (land_frac(i,j,t) <= 0.15) THEN + nrd = 1 + ELSE + nrd = rd + 1 + ENDIF + + IF (nrd > max_rd) RETURN + + IF (done(i,j,t)) RETURN + IF (land_frac(i,j,t) < cutoff) THEN + done(i,j,t) = .true. + inland(i,j,t) = 0.0 + CALL neighbors(t, i, j, nbs) + ! recursively go through k neighbors + DO k = 1, 4 + CALL mark_global_inland_rec_d(nbs%ijt(1,k),nbs%ijt(2,k),nbs%ijt(3,k),nrd) + ENDDO + ENDIF + +END SUBROUTINE mark_global_inland_rec_d + +RECURSIVE SUBROUTINE mark_global_inland_rec(i, j, t) + INTEGER, INTENT(IN) :: i, j, t + + TYPE(nb_gp_idx) :: nbs + INTEGER :: k + + IF (done(i,j,t)) RETURN + IF (land_frac(i,j,t) < 0.9) THEN + done(i,j,t) = .true. + inland(i,j,t) = 0.0 + CALL neighbors(t, i, j, nbs) + ! recursively go through k neighbors + DO k = 1, 4 + CALL mark_global_inland_rec(nbs%ijt(1,k), nbs%ijt(2,k), nbs%ijt(3,k)) + ENDDO + ENDIF + +END SUBROUTINE mark_global_inland_rec + +SUBROUTINE read_orog(cs_res) + USE netcdf + INTEGER, INTENT(IN) :: cs_res + + INTEGER :: tile_sz, tile_num + INTEGER :: stat, ncid, x_dimid, y_dimid, varid + INTEGER :: land_frac_id, slmsk_id, geolon_id, geolat_id + CHARACTER(len=256) :: filename,string + CHARACTER(len=1) :: ich + CHARACTER(len=4) res_ch + + INTEGER :: i, j + REAL, ALLOCATABLE :: var_tmp(:,:) + + tile_sz = cs_res*cs_res + ALLOCATE(var_tmp(cs_res,cs_res)) + + WRITE(res_ch,'(I4)') cs_res + DO tile_num = tile_beg, tile_end + WRITE(ich, '(I1)') tile_num + filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc" + print *,'Read, update, and write ',trim(filename) + stat = nf90_open(filename, NF90_NOWRITE, ncid) + CALL nc_opchk(stat, "nf90_open oro_data.nc") +! original orodata netcdf file uses (y, x) order, so we made change to match it. + stat = nf90_inq_varid(ncid, "land_frac", land_frac_id) + CALL nc_opchk(stat, "nf90_inq_varid: land_frac") + stat = nf90_get_var(ncid, land_frac_id, var_tmp, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_get_var: land_frac") + land_frac(:,:,tile_num) = var_tmp(:,:) + stat = nf90_close(ncid) + CALL nc_opchk(stat, "nf90_close oro_data.nc") + ENDDO + + DEALLOCATE(var_tmp) + +END SUBROUTINE read_orog + +SUBROUTINE write_inland(cs_res) + USE netcdf + INTEGER, INTENT(IN) :: cs_res + + CHARACTER(len=256) :: filename + CHARACTER(len=1) :: ich + CHARACTER(len=4) res_ch + + INTEGER :: tile_num + INTEGER :: stat, ncid, x_dimid, y_dimid, inland_id, dimids(2) + REAL, ALLOCATABLE :: var_tmp(:,:) + + ALLOCATE(var_tmp(cs_res,cs_res)) + + WRITE(res_ch,'(I4)') cs_res + DO tile_num = tile_beg, tile_end + WRITE(ich, '(I1)') tile_num + filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc" + print *,'write inland to ',trim(filename) + stat = nf90_open(filename, NF90_WRITE, ncid) + CALL nc_opchk(stat, "nf90_open oro_data.nc") + stat = nf90_inq_dimid(ncid, "lon", x_dimid) + CALL nc_opchk(stat, "nf90_inq_dim: x") + stat = nf90_inq_dimid(ncid, "lat", y_dimid) + CALL nc_opchk(stat, "nf90_inq_dim: y") + +! original orodata netcdf file uses (y, x) order, so we made change to match it. + dimids = (/ x_dimid, y_dimid /) + +! define a new variables + stat = nf90_redef(ncid) + CALL nc_opchk(stat, "nf90_redef") + stat = nf90_def_var(ncid,"inland",NF90_FLOAT,dimids,inland_id) + CALL nc_opchk(stat, "nf90_def_var: inland") + stat = nf90_put_att(ncid, inland_id,'coordinates','geolon geolat') + CALL nc_opchk(stat, "nf90_put_att: inland:coordinates") + stat = nf90_put_att(ncid, inland_id,'description', & + 'inland = 1 indicates grid cells away from coast') + CALL nc_opchk(stat, "nf90_put_att: inland:description") + + stat = nf90_enddef(ncid) + CALL nc_opchk(stat, "nf90_enddef") + + var_tmp(:,:) = inland(:,:,tile_num) + stat = nf90_put_var(ncid, inland_id, var_tmp, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_put_var: inland") + + stat = nf90_close(ncid) + CALL nc_opchk(stat, "nf90_close oro_data.nc") + ENDDO + DEALLOCATE(var_tmp) + +END SUBROUTINE write_inland + +SUBROUTINE nc_opchk(stat,opname) + USE netcdf + IMPLICIT NONE + INTEGER stat + CHARACTER(len=*) opname + CHARACTER(64) msg + + IF (stat .NE.0) THEN + msg = trim(opname) // ' Error, status code and message:' + PRINT*,trim(msg), stat, nf90_strerror(stat) + STOP + END IF + +END SUBROUTINE nc_opchk + +END PROGRAM inland_mask + diff --git a/sorc/orog_mask_tools.fd/inland.fd/nb.F90 b/sorc/orog_mask_tools.fd/inland.fd/nb.F90 new file mode 100644 index 000000000..32b77abdf --- /dev/null +++ b/sorc/orog_mask_tools.fd/inland.fd/nb.F90 @@ -0,0 +1,335 @@ +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +! This module contains the subroutines that find the adjacent neighbors of a +! given cell, in a cubed sphere grid. Each neighbor is in the form of (i,j,tile). +! +! Ning Wang, July 1, 2020, original version. +! +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + +MODULE cs_nb + IMPLICIT NONE + + TYPE nb_tile_idx + INTEGER :: nb_tile_num + CHARACTER (1) :: nb_tile_bndry + END TYPE nb_tile_idx + + TYPE nb_gp_idx + INTEGER :: gp_type + INTEGER :: ijt(3,8) + END TYPE nb_gp_idx + + TYPE(nb_tile_idx):: nb_tile(4,6) + TYPE(nb_gp_idx):: nb_index + + INTEGER :: cres + + +CONTAINS + +! _______1_______ +! | | 1-upper, 2-bottom, 3-left, 4-right +! | | +! | | +! 3| |4 +! | | +! | | +! |_______________| +! 2 +! Figure 1. Boundary numbers +! + SUBROUTINE idx_init(cres_in) + INTEGER :: cres_in + + nb_tile(1,1)%nb_tile_num = 3; nb_tile(1,1)%nb_tile_bndry ='l' + nb_tile(2,1)%nb_tile_num = 6; nb_tile(2,1)%nb_tile_bndry ='t' + nb_tile(3,1)%nb_tile_num = 5; nb_tile(3,1)%nb_tile_bndry ='t' + nb_tile(4,1)%nb_tile_num = 2; nb_tile(4,1)%nb_tile_bndry ='l' + + nb_tile(1,2)%nb_tile_num = 3; nb_tile(1,2)%nb_tile_bndry ='b' + nb_tile(2,2)%nb_tile_num = 6; nb_tile(2,2)%nb_tile_bndry ='r' + nb_tile(3,2)%nb_tile_num = 1; nb_tile(3,2)%nb_tile_bndry ='r' + nb_tile(4,2)%nb_tile_num = 4; nb_tile(4,2)%nb_tile_bndry ='b' + + nb_tile(1,3)%nb_tile_num = 5; nb_tile(1,3)%nb_tile_bndry ='l' + nb_tile(2,3)%nb_tile_num = 2; nb_tile(2,3)%nb_tile_bndry ='t' + nb_tile(3,3)%nb_tile_num = 1; nb_tile(3,3)%nb_tile_bndry ='t' + nb_tile(4,3)%nb_tile_num = 4; nb_tile(4,3)%nb_tile_bndry ='l' + + nb_tile(1,4)%nb_tile_num = 5; nb_tile(1,4)%nb_tile_bndry ='b' + nb_tile(2,4)%nb_tile_num = 2; nb_tile(2,4)%nb_tile_bndry ='r' + nb_tile(3,4)%nb_tile_num = 3; nb_tile(3,4)%nb_tile_bndry ='r' + nb_tile(4,4)%nb_tile_num = 6; nb_tile(4,4)%nb_tile_bndry ='b' + + nb_tile(1,5)%nb_tile_num = 1; nb_tile(1,5)%nb_tile_bndry ='l' + nb_tile(2,5)%nb_tile_num = 4; nb_tile(2,5)%nb_tile_bndry ='t' + nb_tile(3,5)%nb_tile_num = 3; nb_tile(3,5)%nb_tile_bndry ='t' + nb_tile(4,5)%nb_tile_num = 6; nb_tile(4,5)%nb_tile_bndry ='l' + + nb_tile(1,6)%nb_tile_num = 1; nb_tile(1,6)%nb_tile_bndry ='b' + nb_tile(2,6)%nb_tile_num = 4; nb_tile(2,6)%nb_tile_bndry ='r' + nb_tile(3,6)%nb_tile_num = 5; nb_tile(3,6)%nb_tile_bndry ='r' + nb_tile(4,6)%nb_tile_num = 2; nb_tile(4,6)%nb_tile_bndry ='b' + + cres = cres_in + + END SUBROUTINE idx_init + + INTEGER FUNCTION bndry(i, j) + INTEGER :: i,j + + bndry = 0 ! no boundary + + IF (j == cres) THEN ! upper boundary + bndry = 1 + IF (i == 1) THEN + bndry = 13 + ELSE IF (i == cres) THEN + bndry = 14 + ENDIF + ELSE IF (j == 1) THEN ! bottom boundary + bndry = 2 + IF (i == 1) THEN + bndry = 23 + ELSE IF (i == cres) THEN + bndry = 24 + ENDIF + ELSE IF (i == 1) THEN ! left boundary + bndry = 3 + ELSE IF (i == cres) THEN ! right boundary + bndry = 4 + ENDIF + + END FUNCTION bndry + +! ______________ +! | | | | ________ +! | 5 | 1 | 6 | /\ 1 \ 6 \ +! |____|____|____| / \___\___\ +! | | | | /\2 / c / 3 / +! | 2 | c | 3 | / \/___/___/ +! |____|____|____| \7 / 4 / 8 / +! | | | | \/___/___/ +! | 7 | 4 | 8 | +! |____|____|____| +! +! Figure 2. Eight neighbors of cell 'c' and special cases at upper left +! cornner of the tile +! + SUBROUTINE neighbors(tile, i, j, nb) + INTEGER :: tile, i, j + TYPE(nb_gp_idx) :: nb + + INTEGER :: bd, nb_t_num + + nb%gp_type = bndry(i,j) + IF (nb%gp_type == 0) THEN ! interior (non-boundary) cell + ! top, bottom, left, and right + nb%ijt(1,1) = i; nb%ijt(2,1) = j+1; nb%ijt(3,1) = tile + nb%ijt(1,2) = i-1; nb%ijt(2,2) = j; nb%ijt(3,2) = tile + nb%ijt(1,3) = i+1; nb%ijt(2,3) = j; nb%ijt(3,3) = tile + nb%ijt(1,4) = i; nb%ijt(2,4) = j-1; nb%ijt(3,4) = tile + ! top left, top right, bottom left, and bottom right + nb%ijt(1,5) = i-1; nb%ijt(2,5) = j+1; nb%ijt(3,5) = tile + nb%ijt(1,6) = i+1; nb%ijt(2,6) = j+1; nb%ijt(3,6) = tile + nb%ijt(1,7) = i-1; nb%ijt(2,7) = j-1; nb%ijt(3,7) = tile + nb%ijt(1,8) = i+1; nb%ijt(2,8) = j-1; nb%ijt(3,8) = tile + ELSEIF (nb%gp_type == 1) THEN ! top boundary cell + bd = 1 + nb_t_num = nb_tile(nb%gp_type,tile)%nb_tile_num + nb%ijt(3,1)=nb_t_num; nb%ijt(3,5)=nb_t_num; nb%ijt(3,6)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'l') THEN + nb%ijt(1,1) = 1; nb%ijt(2,1) = cres+1-i; + nb%ijt(1,5) = 1; nb%ijt(2,5) = cres+1-(i-1); + nb%ijt(1,6) = 1; nb%ijt(2,6) = cres+1-(i+1); + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 'b') THEN + nb%ijt(1,1) = i; nb%ijt(2,1) = 1 + nb%ijt(1,5) = i-1; nb%ijt(2,5) = 1 + nb%ijt(1,6) = i+1; nb%ijt(2,6) = 1 + ENDIF + nb%ijt(1,2) = i-1; nb%ijt(2,2) = j; nb%ijt(3,2) = tile + nb%ijt(1,3) = i+1; nb%ijt(2,3) = j; nb%ijt(3,3) = tile + nb%ijt(1,4) = i; nb%ijt(2,4) = j-1; nb%ijt(3,4) = tile + nb%ijt(1,7) = i-1; nb%ijt(2,7) = j-1; nb%ijt(3,7) = tile + nb%ijt(1,8) = i+1; nb%ijt(2,8) = j-1; nb%ijt(3,8) = tile + ELSEIF (nb%gp_type == 2) THEN ! bottom boundary cell + bd = 2 + nb_t_num = nb_tile(nb%gp_type,tile)%nb_tile_num + nb%ijt(3,4)=nb_t_num; nb%ijt(3,7)=nb_t_num; nb%ijt(3,8)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'r') THEN + nb%ijt(1,4) = cres; nb%ijt(2,4) = cres+1-i; + nb%ijt(1,7) = cres; nb%ijt(2,7) = cres+1-(i-1); + nb%ijt(1,8) = cres; nb%ijt(2,8) = cres+1-(i+1); + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 't') THEN + nb%ijt(1,4) = i; nb%ijt(2,4) = cres + nb%ijt(1,7) = i-1; nb%ijt(2,7) = cres + nb%ijt(1,8) = i+1; nb%ijt(2,8) = cres + ENDIF + nb%ijt(1,1) = i; nb%ijt(2,1) = j+1; nb%ijt(3,1) = tile + nb%ijt(1,2) = i-1; nb%ijt(2,2) = j; nb%ijt(3,2) = tile + nb%ijt(1,3) = i+1; nb%ijt(2,3) = j; nb%ijt(3,3) = tile + nb%ijt(1,5) = i-1; nb%ijt(2,5) = j+1; nb%ijt(3,5) = tile + nb%ijt(1,6) = i+1; nb%ijt(2,6) = j+1; nb%ijt(3,6) = tile + ELSEIF (nb%gp_type == 3) THEN ! left boundary cell + bd = 3 + nb_t_num = nb_tile(nb%gp_type,tile)%nb_tile_num + nb%ijt(3,2)=nb_t_num; nb%ijt(3,5)=nb_t_num; nb%ijt(3,7)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'r') THEN + nb%ijt(1,2) = cres; nb%ijt(2,2) = j; + nb%ijt(1,5) = cres; nb%ijt(2,5) = j+1; + nb%ijt(1,7) = cres; nb%ijt(2,7) = j-1; + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 't') THEN + nb%ijt(1,2) = cres+1-j; nb%ijt(2,2) = cres + nb%ijt(1,5) = cres+1-(j+1); nb%ijt(2,5) = cres + nb%ijt(1,7) = cres+1-(j-1); nb%ijt(2,7) = cres + ENDIF + nb%ijt(1,1) = i; nb%ijt(2,1) = j+1; nb%ijt(3,1) = tile + nb%ijt(1,3) = i+1; nb%ijt(2,3) = j; nb%ijt(3,3) = tile + nb%ijt(1,4) = i; nb%ijt(2,4) = j-1; nb%ijt(3,4) = tile + nb%ijt(1,6) = i+1; nb%ijt(2,6) = j+1; nb%ijt(3,6) = tile + nb%ijt(1,8) = i+1; nb%ijt(2,8) = j-1; nb%ijt(3,8) = tile + ELSEIF (nb%gp_type == 4) THEN ! right boundary cell + bd = 4 + nb_t_num = nb_tile(nb%gp_type,tile)%nb_tile_num + nb%ijt(3,3)=nb_t_num; nb%ijt(3,6)=nb_t_num; nb%ijt(3,8)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'l') THEN + nb%ijt(1,3) = 1; nb%ijt(2,3) = j; + nb%ijt(1,6) = 1; nb%ijt(2,6) = j+1; + nb%ijt(1,8) = 1; nb%ijt(2,8) = j-1; + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 'b') THEN + nb%ijt(1,3) = cres+1-j; nb%ijt(2,3) = 1 + nb%ijt(1,6) = cres+1-(j+1); nb%ijt(2,6) = 1 + nb%ijt(1,8) = cres+1-(j-1); nb%ijt(2,8) = 1 + ENDIF + nb%ijt(1,1) = i; nb%ijt(2,1) = j+1; nb%ijt(3,1) = tile + nb%ijt(1,2) = i-1; nb%ijt(2,2) = j; nb%ijt(3,2) = tile + nb%ijt(1,4) = i; nb%ijt(2,4) = j-1; nb%ijt(3,4) = tile + nb%ijt(1,5) = i-1; nb%ijt(2,5) = j+1; nb%ijt(3,5) = tile + nb%ijt(1,7) = i-1; nb%ijt(2,7) = j-1; nb%ijt(3,7) = tile + ELSEIF (nb%gp_type == 13) THEN ! upper left coner + bd = 1 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,1)=nb_t_num; nb%ijt(3,6)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'l') THEN + nb%ijt(1,1) = 1; nb%ijt(2,1) = cres+1-i + nb%ijt(1,6) = 1; nb%ijt(2,6) = cres+1-(i+1) + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 'b') THEN + nb%ijt(1,1) = i; nb%ijt(2,1) = 1 + nb%ijt(1,6) = i+1; nb%ijt(2,6) = 1 + ENDIF + bd = 3 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,2)=nb_t_num; nb%ijt(3,7)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'r') THEN + nb%ijt(1,2) = cres; nb%ijt(2,2) = j + nb%ijt(1,7) = cres; nb%ijt(2,7) = j-1 + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 't') THEN + nb%ijt(1,2) = cres+1-j; nb%ijt(2,2) = cres + nb%ijt(1,7) = cres+1-(j-1); nb%ijt(2,7) = cres + ENDIF + nb%ijt(3,5)=0 + nb%ijt(1,3) = i+1; nb%ijt(2,3) = j; nb%ijt(3,3) = tile + nb%ijt(1,4) = i; nb%ijt(2,4) = j-1; nb%ijt(3,4) = tile + nb%ijt(1,8) = i+1; nb%ijt(2,8) = j-1; nb%ijt(3,8) = tile + ELSEIF (nb%gp_type == 14) THEN ! upper left coner + bd = 1 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,1)=nb_t_num; nb%ijt(3,5)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'l') THEN + nb%ijt(1,1) = 1; nb%ijt(2,1) = cres+1-i + nb%ijt(1,5) = 1; nb%ijt(2,5) = cres+1-(i-1) + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 'b') THEN + nb%ijt(1,1) = i; nb%ijt(2,1) = 1 + nb%ijt(1,5) = i-1; nb%ijt(2,5) = 1 + ENDIF + bd = 4 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,3)=nb_t_num; nb%ijt(3,8)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'l') THEN + nb%ijt(1,3) = 1; nb%ijt(2,3) = j + nb%ijt(1,8) = 1; nb%ijt(2,8) = j-1 + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 'b') THEN + nb%ijt(1,3) = cres+1-j; nb%ijt(2,3) = 1 + nb%ijt(1,8) = cres+1-(j-1); nb%ijt(2,8) = 1 + ENDIF + nb%ijt(3,6)=0 + nb%ijt(1,2) = i-1; nb%ijt(2,2) = j; nb%ijt(3,2) = tile + nb%ijt(1,4) = i; nb%ijt(2,4) = j-1; nb%ijt(3,4) = tile + nb%ijt(1,7) = i-1; nb%ijt(2,7) = j-1; nb%ijt(3,7) = tile + ELSEIF (nb%gp_type == 23) THEN ! upper left coner + bd = 2 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,4)=nb_t_num; nb%ijt(3,8)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'r') THEN + nb%ijt(1,4) = cres; nb%ijt(2,4) = cres+1-i + nb%ijt(1,8) = cres; nb%ijt(2,8) = cres+1-(i+1) + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 't') THEN + nb%ijt(1,4) = i; nb%ijt(2,4) = cres + nb%ijt(1,8) = i+1; nb%ijt(2,8) = cres + ENDIF + bd = 3 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,2)=nb_t_num; nb%ijt(3,5)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'r') THEN + nb%ijt(1,2) = cres; nb%ijt(2,2) = j + nb%ijt(1,5) = cres; nb%ijt(2,5) = j+1 + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 't') THEN + nb%ijt(1,2) = cres+1-j; nb%ijt(2,2) = cres + nb%ijt(1,5) = cres+1-(j+1); nb%ijt(2,5) = cres + ENDIF + nb%ijt(3,7)=0 + nb%ijt(1,1) = i; nb%ijt(2,1) = j+1; nb%ijt(3,1) = tile + nb%ijt(1,3) = i+1; nb%ijt(2,3) = j; nb%ijt(3,3) = tile + nb%ijt(1,6) = i+1; nb%ijt(2,6) = j+1; nb%ijt(3,6) = tile + ELSEIF (nb%gp_type == 24) THEN ! upper left coner + bd = 2 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,4)=nb_t_num; nb%ijt(3,7)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'r') THEN + nb%ijt(1,4) = cres; nb%ijt(2,4) = cres+1-i + nb%ijt(1,7) = cres; nb%ijt(2,7) = cres+1-(i-1) + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 't') THEN + nb%ijt(1,4) = i; nb%ijt(2,4) = cres + nb%ijt(1,7) = i-1; nb%ijt(2,7) = cres + ENDIF + bd = 4 + nb_t_num = nb_tile(bd,tile)%nb_tile_num + nb%ijt(3,3)=nb_t_num; nb%ijt(3,6)=nb_t_num + IF (nb_tile(bd,tile)%nb_tile_bndry == 'l') THEN + nb%ijt(1,3) = 1; nb%ijt(2,3) = j + nb%ijt(1,6) = 1; nb%ijt(2,6) = j+1 + ELSEIF (nb_tile(bd,tile)%nb_tile_bndry == 'b') THEN + nb%ijt(1,3) = cres+1-j; nb%ijt(2,3) = 1 + nb%ijt(1,6) = cres+1-(j+1); nb%ijt(2,6) = 1 + ENDIF + nb%ijt(3,8)=0 + nb%ijt(1,1) = i; nb%ijt(2,1) = j+1; nb%ijt(3,1) = tile + nb%ijt(1,2) = i-1; nb%ijt(2,2) = j; nb%ijt(3,2) = tile + nb%ijt(1,5) = i-1; nb%ijt(2,5) = j+1; nb%ijt(3,5) = tile + + ENDIF + END SUBROUTINE neighbors + +END MODULE cs_nb + +#ifdef TEST_CS_NB +PROGRAM test_nb + USE cs_nb + + INTEGER :: tile, i, j + TYPE(nb_gp_idx) :: nbs + + INTEGER, PARAMETER :: res = 96 + + CALL idx_init(res) +! tile = 1; i = 1; j = 5 +! tile = 1; i = 5; j = 1 +! tile = 1; i = 96; j = 10 + tile = 2; i = 96; j = 96 + CALL neighbors(tile, i, j, nbs) + print*, 'tile = ', tile, ' i = ',i, ' j = ', j + print*, 'nbs%type: ', nbs%gp_type + print*, 'nbs%ijt: ' + print*, nbs%ijt +END PROGRAM test_nb +#endif diff --git a/sorc/orog_mask_tools.fd/lake.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/lake.fd/CMakeLists.txt new file mode 100644 index 000000000..b33ddfd6d --- /dev/null +++ b/sorc/orog_mask_tools.fd/lake.fd/CMakeLists.txt @@ -0,0 +1,20 @@ +set(fortran_src + enclosure_cnvx.F90 + find_limit.F90 + lakefrac.F90) + +if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -assume byterecl") +elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU|Clang|AppleClang)$") + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fno-range-check") +endif() + +set(exe_name lakefrac) + +add_executable(${exe_name} ${fortran_src}) + +target_link_libraries( + ${exe_name} + NetCDF::NetCDF_Fortran) + +install(TARGETS ${exe_name} RUNTIME DESTINATION ${exec_dir}) diff --git a/sorc/orog_mask_tools.fd/lake.fd/enclosure_cnvx.F90 b/sorc/orog_mask_tools.fd/lake.fd/enclosure_cnvx.F90 new file mode 100644 index 000000000..35bd766ee --- /dev/null +++ b/sorc/orog_mask_tools.fd/lake.fd/enclosure_cnvx.F90 @@ -0,0 +1,82 @@ +#ifdef INCLUDE_TEST_DRIVER + PROGRAM testenc + IMPLICIT NONE + REAL*8 :: v(2,4) + REAL*8 :: p(2) + + REAL*8 :: d2r + LOGICAL:: enclosure_cnvx, inside + INTEGER :: co_gc + + d2r = acos(-1.0)/180.0D0 + + v(1,1) = 10.0D0*d2r; v(2,1) = 20.0D0*d2r + v(1,2) = 15.0D0*d2r; v(2,2) = 30.0D0*d2r + v(1,3) = 17.7D0*d2r; v(2,3) = 25.0D0*d2r + v(1,4) = 20.0D0*d2r; v(2,4) = 20.0D0*d2r + +! p(1) = 15.0D0*d2r; p(2) = 30.00000001D0*d2r +! p(1) = 20.00000000D0*d2r; p(2) = 20.0D0*d2r +! p(1) = 9.999999999D0*d2r; p(2) = 20.0D0*d2r +! p(1) = 10.00000000*d2r; p(2) = 20.000000001D0*d2r + p(1) = 17.7D0*d2r; p(2) = 25.000000001D0*d2r + + inside = enclosure_cnvx(v,4,p,co_gc) + IF (inside) THEN + PRINT*, 'inside ', co_gc + ELSE + PRINT*, 'outside ', co_gc + ENDIF + + END PROGRAM +#endif + + +LOGICAL FUNCTION enclosure_cnvx(v, n, p, co_gc) +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +! Function enclosure tests whether a given point p(2) is inside a convex spherical polygon +! defined with a sequence of n vertices v(2,n). Both the test point and the polygon are +! apecified in latitude and longitude radians. +! +! N. Wang, Jan 2007 - Initial version +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> + IMPLICIT NONE + REAL*8, INTENT(IN) :: v(2,n), p(2) + INTEGER, INTENT(IN) :: n + INTEGER, INTENT(OUT) :: co_gc + + REAL*8 v_xy(2, n) + REAL*8 cp_z(n), cos_d2c, eps + + INTEGER :: i, ip1 + + + eps = 0.000000000000001D0 + co_gc = 0 + DO i = 1, n + cos_d2c = sin(p(1))*sin(v(1,i)) + cos(p(1))*cos(v(1,i))*cos(v(2,i)-p(2)) + v_xy(1,i) = (cos(v(1,i))*sin(v(2,i)-p(2)))/cos_d2c + v_xy(2,i) = (cos(p(1))*sin(v(1,i))-sin(p(1))*cos(v(1,i))*cos(v(2,i)-p(2)))/cos_d2c + + ENDDO + + DO i = 1, n + ip1 = mod(i,n)+1 + cp_z(i) = v_xy(1,i)*v_xy(2,ip1)-v_xy(2,i)*v_xy(1,ip1) + IF (abs(cp_z(i)) < eps) co_gc = i + ENDDO + + DO i = 1, n-1 + ip1 = mod(i,n)+1 + IF (cp_z(i)*cp_z(ip1) .LT. -eps) THEN + enclosure_cnvx = .false. + RETURN + ENDIF + ENDDO + + enclosure_cnvx = .true. + RETURN + +END FUNCTION enclosure_cnvx + +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> diff --git a/sorc/orog_mask_tools.fd/lake.fd/find_limit.F90 b/sorc/orog_mask_tools.fd/lake.fd/find_limit.F90 new file mode 100644 index 000000000..cc6964f59 --- /dev/null +++ b/sorc/orog_mask_tools.fd/lake.fd/find_limit.F90 @@ -0,0 +1,100 @@ +!#define DIAG +SUBROUTINE find_limit (p1_in, p2_in, latmin, latmax) + REAL*8, INTENT(IN) :: p1_in(2), p2_in(2) + REAL*8, INTENT(OUT) :: latmin, latmax + + REAL*8 :: p1(2),p2(2), pm(2) + REAL*8 :: r2d = 180.0/acos(-1.0) + + p1 = p1_in/r2d; p2 = p2_in/r2d + latmin = min(p1(1), p2(1)) + latmax = max(p1(1), p2(1)) + + CALL middle(p1, p2, pm) +#ifdef DIAG + PRINT*, 'before loop', p1(1)*r2d,p2(1)*r2d,pm(1)*r2d +#endif + + DO WHILE (abs(p1(1)-p2(1)) > 0.00001 .AND. & + abs(p1(2)-p2(2)) > 0.00001 ) + IF (abs(p1(1)-pm(1)) < abs(p2(1)-pm(1))) THEN + p2 = pm + ELSE + p1 = pm + ENDIF + CALL middle(p1, p2, pm) +#ifdef DIAG + PRINT*, 'in loop', p1(1)*r2d,p2(1)*r2d, pm(1)*r2d +#endif + ENDDO + + latmin = min(latmin, pm(1)) + latmax = max(latmax, pm(1)) + + latmin = latmin *r2d + latmax = latmax *r2d + +END SUBROUTINE find_limit + +!====================================================== +! This subroutine computes the latitude and longitude +! of the middle point between two given ponits. +! +! There are two formulae available to compute it. +! +! One derived from a more general m-sect formula: +! +! xyz = sin((1-f)*theta) / sin(theta) * xyz1 + +! sin(f*theta) /sin(theta) * xyz2 ; +! where theta is the angle of xyz1, and xyz2. +! +! xyz = 0.5 / sqrt[(1+dot(xyz1,xyz2))/2] * (xyz1+xyz2) +! +! and the other one is the normalized middle point of +! the two end points: +! +! xyz = 0.5 * (xyz1+xyz2), xyz = xyz / sqrt(dot(xyz,xyz)) +! +! Author: Ning Wang, March, 2006 +!====================================================== + +SUBROUTINE middle(p1,p2,p) + IMPLICIT NONE + + ! Two given points in lat/lon: + REAL*8, INTENT(IN) :: p1(2),p2(2) + REAL*8, INTENT(OUT) :: p(2) + REAL*8 :: pi + + REAL*8 :: xyz1(3),xyz2(3),xyz(3) + + pi = acos(-1.0) + + ! Convert them into Cardesian coor: + xyz1(1) = cos(p1(1)) * cos(p1(2)) + xyz1(2) = cos(p1(1)) * sin(p1(2)) + xyz1(3) = sin(p1(1)) + + xyz2(1) = cos(p2(1)) * cos(p2(2)) + xyz2(2) = cos(p2(1)) * sin(p2(2)) + xyz2(3) = sin(p2(1)) + + ! middle point: + +! coeff = 0.5 / sqrt((1.0 + dot_product(xyz1,xyz2)) / 2) +! xyz = coeff * (xyz1 + xyz2) + + xyz = 0.5 * (xyz1 + xyz2) + + xyz = xyz / sqrt(dot_product(xyz,xyz)) + + ! Convert the middle point to lat/lon coor: + p(1) = atan2(xyz(3), sqrt(xyz(1) * xyz(1) + xyz(2) * xyz(2))) + p(2) = atan2(xyz(2), xyz(1)) + + IF (p(2) < -pi / 2.0) THEN + p(2) = p(2) + 2 * pi + END IF + +END SUBROUTINE middle + diff --git a/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 new file mode 100644 index 000000000..9673dacfb --- /dev/null +++ b/sorc/orog_mask_tools.fd/lake.fd/lakefrac.F90 @@ -0,0 +1,896 @@ +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +! This program computes lake fraction and depth numbers for FV3 cubed sphere +! grid cells, from a high resolution lat/lon data set. +! +! Ning Wang, July 2018: Original version +! +! Shan Sun, Aug. 2018: Added Caspian Sea and Aral Sea to the lake fraction +! and lake depth fields. +! Shan Sun, Dec. 2018: Added round up and round down with respect to a +! numerical minimum value and a cut-off value, for lake +! fraction number. +! Ning Wang, Apr. 2019: Extended the program to process the same lake data +! for FV3 stand-alone regional (SAR) model. +! +!<><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><><> +!#define DIAG_N_VERBOSE +#define ADD_ATT_FOR_NEW_VAR +PROGRAM lake_frac + USE netcdf + IMPLICIT NONE + + CHARACTER(len=256) :: sfcdata_path + INTEGER :: cs_res, ncsmp, ncscp, i + INTEGER :: res_x, res_y + + INTEGER*1, ALLOCATABLE :: lakestatus(:) + INTEGER*2, ALLOCATABLE :: lakedepth(:) + REAL, ALLOCATABLE :: cs_grid(:,:) + REAL, ALLOCATABLE :: cs_lakestatus(:), cs_lakedepth(:) + REAL, ALLOCATABLE :: src_grid_lon(:), src_grid_lat(:) + + INTEGER :: tile_req, tile_beg, tile_end + REAL :: lake_cutoff + + INTEGER, PARAMETER :: nlat = 21600, nlon = 43200 + REAL, PARAMETER :: d2r = acos(-1.0) / 180.0 + REAL, PARAMETER :: r2d = 180.0 /acos(-1.0) + REAL, PARAMETER :: pi = acos(-1.0) + REAL*8, PARAMETER :: gppdeg = 119.99444445 + REAL*8, PARAMETER :: delta = 1.0 / 119.99444445 + + CHARACTER(len=32) :: arg + CHARACTER(len=256) :: lakedata_path + INTEGER :: stat + + CALL getarg(0, arg) ! get the program name + IF (iargc() /= 3 .AND. iargc() /= 4) THEN + PRINT*, 'Usage: ', trim(arg), & + ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [path to lake data file]' + PRINT*, 'Or: ', trim(arg), & + ' [tile_num (0:all tiles, 7:regional)] [resolution (48,96, ...)] [path to lake data file] [lake_cutoff]' + STOP + ENDIF + CALL getarg(1, arg) + READ(arg,*,iostat=stat) tile_req + CALL getarg(2, arg) + READ(arg,*,iostat=stat) cs_res + CALL getarg(3, lakedata_path) + + IF (iargc() == 3) THEN + lake_cutoff = 0.20 + ELSE + CALL getarg(4, arg) + READ(arg,*,iostat=stat) lake_cutoff + ENDIF + + IF (tile_req == 0) THEN + tile_beg = 1; tile_end = 6 + PRINT*, 'Process tile 1 - 6 at resolution C',cs_res + ELSE IF (tile_req /= 7) THEN + tile_beg = tile_req; tile_end = tile_req + PRINT*, 'Process tile',tile_req, 'at resolution C',cs_res + ELSE + tile_beg = 1; tile_end = 1 + PRINT*, 'Process regional tile (tile', tile_req, ') at resolution C',cs_res + ENDIF + +! cs_res = 96 ! tile_beg = 3; tile_end = 3 + + ! read in grid spec data for each tile and concatenate them together + + ncsmp = (2*cs_res+1)*(2*cs_res+1)*6 + IF (tile_req /= 7) THEN + PRINT*, 'Read in cubed sphere grid information ... ',ncsmp,'pairs of lat/lons' + ENDIF + + IF (tile_req /= 7) THEN + ALLOCATE(cs_grid(ncsmp, 2)) + CALL read_cubed_sphere_grid(cs_res, cs_grid) + ELSE + CALL read_cubed_sphere_reg_grid(cs_res, cs_grid, 3, res_x, res_y) + ENDIF + ! compute source grid + ALLOCATE(src_grid_lon(nlon), src_grid_lat(nlat)) + DO i = 1, nlon + src_grid_lon(i) = -180.0 + delta * (i-1) + ENDDO + DO i = 1, nlat + src_grid_lat(i) = 90.0 - delta * (i-1) + ENDDO + + ! read in lake data file +! sfcdata_path = '/scratch1/NCEPDEV/global/glopara/fix/fix_orog/' + lakedata_path = trim(lakedata_path) // "/" + ALLOCATE(lakestatus(nlon*nlat),lakedepth(nlon*nlat)) + PRINT*, 'Read in lake data file ...' + CALL read_lakedata(lakedata_path,lakestatus,lakedepth,nlat,nlon) + + ! calculate fraction numbers for all cs-cells + ncscp = cs_res*cs_res*6 + ALLOCATE(cs_lakestatus(ncscp)) + ALLOCATE(cs_lakedepth(ncscp)) + + PRINT*, 'Calculate lake fraction and average depth for cubed-sphere cells ...' + CALL cal_lake_frac_depth(lakestatus,cs_lakestatus,lakedepth,cs_lakedepth) + + ! write lake status (in terms of fraction) and lake depth(average) + ! to a netcdf file + IF (tile_req /= 7) THEN + PRINT*, 'Write lake fraction/depth on cubed sphere grid cells to netCDF files ...' + CALL write_lakedata_to_orodata(cs_res, cs_lakestatus, cs_lakedepth) + ELSE + PRINT*, 'Write lake fraction/depth on regional FV3 grid cells to a netCDF file ...' + CALL write_reg_lakedata_to_orodata(cs_res, res_x, res_y, cs_lakestatus, cs_lakedepth) + ENDIF + + DEALLOCATE(cs_lakestatus,cs_lakedepth) + DEALLOCATE(cs_grid) + DEALLOCATE(lakestatus,lakedepth) + STOP +CONTAINS + +SUBROUTINE cal_lake_frac_depth(lakestat,cs_lakestat,lakedpth,cs_lakedpth) + INTEGER*1, INTENT(IN) :: lakestat(:) + INTEGER*2, INTENT(IN) :: lakedpth(:) + REAL, INTENT(OUT) :: cs_lakestat(:), cs_lakedpth(:) + + REAL*8 lolf(2), lort(2), uplf(2), uprt(2), sd_ltmn(4), sd_ltmx(4) + REAL*8 :: v(2,4), p(2) + REAL :: latmin1, latmax1 + REAL :: latmin, latmax, lonmin, lonmax, lontmp, lat_sz_max, lon_sz_max + INTEGER :: tile_num, i, j, gp, row, col, cs_grid_idx, cs_data_idx + INTEGER :: sidex_res, sidey_res, sidex_sz, sidey_sz + INTEGER :: stride_lat, stride_lon + INTEGER :: src_grid_lat_beg,src_grid_lat_end,src_grid_lon_beg,src_grid_lon_end + INTEGER :: src_grid_lon_beg1,src_grid_lon_end1,src_grid_lon_beg2,src_grid_lon_end2 + INTEGER :: grid_ct, lake_ct, co_gc, tmp + + REAL*8 :: lake_dpth_sum + LOGICAL :: two_section, enclosure_cnvx + + IF (tile_req /= 7) THEN + sidex_res = cs_res; sidey_res = cs_res + ELSE + sidex_res = res_x; sidey_res = res_y + ENDIF + + sidex_sz = 2*sidex_res+1; sidey_sz = 2*sidey_res+1 + + stride_lat = 1 + + lat_sz_max = 0.0 + lon_sz_max = 0.0 + + cs_lakestat = 0.0 + + DO tile_num = tile_beg, tile_end + row = 2 + sidex_sz*(tile_num-1); col = 2 + DO gp = 1, sidex_res*sidey_res + two_section = .false. + cs_grid_idx = (row-1)*sidex_sz+col + cs_data_idx = (tile_num-1)*sidex_res*sidey_res+gp + IF (abs(cs_grid(cs_grid_idx,1)) > 80.0 ) THEN !ignore lakes in very high latitude + cs_lakestat(cs_data_idx) = 0.0 + cs_lakedpth(cs_data_idx) = 0.0 + ! move to next cs cell + col = col + 2 + IF (col > sidex_sz) THEN + col = 2 + row = row + 2 + ENDIF + CYCLE + ENDIF + ! get the four corners of the cs cell + lolf(1) = cs_grid(cs_grid_idx-sidex_sz-1, 1) + lolf(2) = cs_grid(cs_grid_idx-sidex_sz-1, 2) + IF (lolf(2) > 180.0) lolf(2) = lolf(2) - 360.0 + lort(1) = cs_grid(cs_grid_idx-sidex_sz+1, 1) + lort(2) = cs_grid(cs_grid_idx-sidex_sz+1, 2) + IF (lort(2) > 180.0) lort(2) = lort(2) - 360.0 + uplf(1) = cs_grid(cs_grid_idx+sidex_sz-1,1) + uplf(2) = cs_grid(cs_grid_idx+sidex_sz-1,2) + IF (uplf(2) > 180.0) uplf(2) = uplf(2) - 360.0 + uprt(1) = cs_grid(cs_grid_idx+sidex_sz+1,1) + uprt(2) = cs_grid(cs_grid_idx+sidex_sz+1,2) + + v(1,1) = lolf(1); v(2,1) = lolf(2) + v(1,2) = lort(1); v(2,2) = lort(2) + v(1,3) = uprt(1); v(2,3) = uprt(2) + v(1,4) = uplf(1); v(2,4) = uplf(2) + v(:,:) = v(:,:) * d2r + + IF (uprt(2) > 180.0) uprt(2) = uprt(2) - 360.0 + ! gather the candidate indices in lakestat +#ifdef LIMIT_CAL + CALL find_limit (lolf, lort, sd_ltmn(1), sd_ltmx(1)) + CALL find_limit (lort, uprt, sd_ltmn(2), sd_ltmx(2)) + CALL find_limit (uprt, uplf, sd_ltmn(3), sd_ltmx(3)) + CALL find_limit (uplf, lolf, sd_ltmn(4), sd_ltmx(4)) + latmin = min(sd_ltmn(1),min(sd_ltmn(2),min(sd_ltmn(3),sd_ltmn(4)))) + latmax = max(sd_ltmx(1),max(sd_ltmx(2),max(sd_ltmx(3),sd_ltmx(4)))) +#endif + latmin = min(lolf(1),min(lort(1),min(uplf(1),uprt(1)))) + latmax = max(lolf(1),max(lort(1),max(uplf(1),uprt(1)))) + lonmin = min(lolf(2),min(lort(2),min(uplf(2),uprt(2)))) + lonmax = max(lolf(2),max(lort(2),max(uplf(2),uprt(2)))) +! lat_sz_max = max(lat_sz_max, (latmax-latmin)) +! lon_sz_max = max(lon_sz_max, (lonmax-lonmin)) + + src_grid_lat_beg = nint((90.0-latmax)*gppdeg+0.5) + src_grid_lat_end = nint((90.0-latmin)*gppdeg+0.5) + src_grid_lon_beg = nint((180.0+lonmin)*gppdeg+0.5) + src_grid_lon_end = nint((180.0+lonmax)*gppdeg+0.5) + + IF (src_grid_lat_beg > src_grid_lat_end) THEN + print*,'switch!!!' + tmp = src_grid_lat_beg + src_grid_lat_beg = src_grid_lat_end + src_grid_lat_end = tmp + ENDIF + IF (src_grid_lon_beg > src_grid_lon_end) THEN + print*,'switch!!!' + tmp = src_grid_lon_beg + src_grid_lon_beg = src_grid_lon_end + src_grid_lon_end = tmp + ENDIF + IF ((src_grid_lon_end - src_grid_lon_beg) > nlon*0.75) THEN + two_section = .true. + src_grid_lon_beg1 = src_grid_lon_end + src_grid_lon_end1 = nlon + src_grid_lon_beg2 = 1 + src_grid_lon_end2 = src_grid_lon_beg + ENDIF + +#ifdef DIAG_N_VERBOSE + PRINT*, 'cell centre lat/lon =', & + gp, cs_res*cs_res, cs_grid(cs_grid_idx,1),cs_grid(cs_grid_idx,2) + PRINT*, 'lat index range and stride', & + src_grid_lat_beg,src_grid_lat_end,stride_lat + PRINT*, 'lat range ', & + src_grid_lat(src_grid_lat_beg),src_grid_lat(src_grid_lat_end) +#endif + lake_ct = 0; grid_ct = 0 + lake_dpth_sum = 0.0 + DO j = src_grid_lat_beg, src_grid_lat_end, stride_lat + stride_lon = int(1.0/cos(src_grid_lat(j)*d2r)*REAL(stride_lat)) +#ifdef DIAG_N_VERBOSE + IF (j == src_grid_lat_beg) THEN + PRINT*, 'lon index range and stride', & + src_grid_lon_beg,src_grid_lon_end,stride_lon + PRINT*, 'lon range ', & + src_grid_lon(src_grid_lon_beg),src_grid_lon(src_grid_lon_end) + IF (two_section == .true.) THEN + PRINT*, 'section1 index lon range and stride', & + src_grid_lon_beg1,src_grid_lon_end1,stride_lon + PRINT*, 'section1 lon range ', & + src_grid_lon(src_grid_lon_beg1),src_grid_lon(src_grid_lon_end1) + PRINT*, 'section2 index lon range and stride', & + src_grid_lon_beg2,src_grid_lon_end2,stride_lon + PRINT*, 'section2 lon range ', & + src_grid_lon(src_grid_lon_beg2),src_grid_lon(src_grid_lon_end2) + ENDIF + ENDIF +#endif + IF (two_section == .false.) THEN + DO i = src_grid_lon_beg, src_grid_lon_end, stride_lon + p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) + p(:) = p(:)*d2r + IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN + grid_ct = grid_ct+1 + IF (lakestat((j-1)*nlon+i) /= 0) THEN + lake_ct = lake_ct+1 + IF (lakedpth((j-1)*nlon+i) < 0) THEN + IF (lakestat((j-1)*nlon+i) == 4) THEN + lake_dpth_sum = lake_dpth_sum+30.0 + ELSE + lake_dpth_sum = lake_dpth_sum+100.0 + ENDIF + ELSE + lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i)) + ENDIF + ENDIF + ENDIF + ENDDO + ELSE + DO i = src_grid_lon_beg1, src_grid_lon_end1, stride_lon + p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) + p(:) = p(:)*d2r + IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN + grid_ct = grid_ct+1 + IF (lakestat((j-1)*nlon+i) /= 0) THEN + lake_ct = lake_ct+1 + IF (lakedpth((j-1)*nlon+i) < 0) THEN + IF (lakestat((j-1)*nlon+i) == 4) THEN + lake_dpth_sum = lake_dpth_sum+30.0 + ELSE + lake_dpth_sum = lake_dpth_sum+100.0 + ENDIF + ELSE + lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i)) + ENDIF + ENDIF + ENDIF + ENDDO + DO i = src_grid_lon_beg2, src_grid_lon_end2, stride_lon + p(1) = src_grid_lat(j); p(2) = src_grid_lon(i) + p(:) = p(:)*d2r + IF(enclosure_cnvx(v, 4, p, co_gc) == .true.) THEN + grid_ct = grid_ct+1 + IF (lakestat((j-1)*nlon+i) /= 0) THEN + lake_ct = lake_ct+1 + IF (lakedpth((j-1)*nlon+i) < 0) THEN + IF (lakestat((j-1)*nlon+i) == 4) THEN + lake_dpth_sum = lake_dpth_sum+30.0 + ELSE + lake_dpth_sum = lake_dpth_sum+100.0 + ENDIF + ELSE + lake_dpth_sum = lake_dpth_sum+REAL(lakedpth((j-1)*nlon+i)) + ENDIF + ENDIF + ENDIF + ENDDO + ENDIF + ENDDO + cs_lakestat(cs_data_idx)=REAL(lake_ct)/REAL(grid_ct) + IF (lake_ct /= 0) THEN + cs_lakedpth(cs_data_idx)=lake_dpth_sum/REAL(lake_ct)/10.0 !convert to meter + ELSE + cs_lakedpth(cs_data_idx)=0.0 + ENDIF +#ifdef DIAG_N_VERBOSE + PRINT*, 'tile_num, row, col:', tile_num, row, col + PRINT*, 'grid_ct, lake_ct = ', grid_ct, lake_ct + PRINT*, 'lake_frac= ', cs_lakestat(cs_data_idx) + PRINT*, 'lake_depth (avg) = ', cs_lakedpth(cs_data_idx) +#endif + + ! move to the next control volume + col = col + 2 + IF (col > sidex_sz) THEN + col = 2 + row = row + 2 + ENDIF + ENDDO + PRINT "('*'$)" ! progress '*' + ENDDO + PRINT*, '' + +END SUBROUTINE cal_lake_frac_depth + + +SUBROUTINE read_cubed_sphere_grid(res, grid) + INTEGER, INTENT(IN) :: res + REAL, INTENT(OUT) :: grid(:,:) + + REAL*8, ALLOCATABLE :: lat(:), lon(:) + INTEGER :: side_sz, tile_sz, ncid, varid + INTEGER :: i, tile_beg, tile_end, stat + CHARACTER(len=256) :: gridfile_path,gridfile + CHARACTER(len=1) ich + CHARACTER(len=4) res_ch + + side_sz = 2*res+1 + tile_sz = side_sz*side_sz + ALLOCATE(lat(tile_sz), lon(tile_sz)) + + IF (tile_req == 0) THEN + tile_beg = 1; tile_end = 6 + ELSE + tile_beg = tile_req; tile_end = tile_req + ENDIF + WRITE(res_ch,'(I4)') res +! gridfile_path = "/scratch1/NCEPDEV/global/glopara/fix/fix_fv3/C"//trim(adjustl(res_ch))//"/" + gridfile_path = "./" + DO i = tile_beg, tile_end + WRITE(ich, '(I1)') i + gridfile = trim(gridfile_path)//"C"//trim(adjustl(res_ch))//"_grid.tile"//ich//".nc" + + PRINT*, 'Open cubed sphere grid spec file ', trim(gridfile) + + stat = nf90_open(trim(gridfile), NF90_NOWRITE, ncid) + CALL nc_opchk(stat,'nf90_open') + + stat = nf90_inq_varid(ncid,'x',varid) + CALL nc_opchk(stat,'nf90_inq_lon') + stat = nf90_get_var(ncid,varid,lon,start=(/1,1/),count=(/side_sz,side_sz/)) + CALL nc_opchk(stat,'nf90_get_var_lon') + stat = nf90_inq_varid(ncid,'y',varid) + CALL nc_opchk(stat,'nf90_inq_lat') + stat = nf90_get_var(ncid,varid,lat,start=(/1,1/),count=(/side_sz,side_sz/)) + CALL nc_opchk(stat,'nf90_get_var_lat') + stat = nf90_close(ncid) + CALL nc_opchk(stat,'nf90_close') + + tile_beg = (i-1)*tile_sz+1 + tile_end = i*tile_sz + grid(tile_beg:tile_end,1) = lat(1:tile_sz) + grid(tile_beg:tile_end,2) = lon(1:tile_sz) + END DO + DEALLOCATE(lat,lon) + +END SUBROUTINE read_cubed_sphere_grid + +SUBROUTINE read_cubed_sphere_reg_grid(res, grid, halo_depth, res_x, res_y) + INTEGER, INTENT(IN) :: res, halo_depth + INTEGER, INTENT(OUT) :: res_x, res_y + REAL, ALLOCATABLE, INTENT(OUT) :: grid(:,:) + + REAL*8, ALLOCATABLE :: lat(:), lon(:) + INTEGER :: sidex_sz, sidey_sz, tile_sz, ncid, varid, dimid + INTEGER :: x_start, y_start + INTEGER :: nxp, nyp, stat + CHARACTER(len=256) :: gridfile_path,gridfile + CHARACTER(len=1) ich + CHARACTER(len=4) res_ch + CHARACTER(len=8) dimname + + + WRITE(res_ch,'(I4)') res +! gridfile_path = '/scratch4/BMC/fim/fv3data/fix/reg/' + gridfile_path = './' + gridfile = trim(gridfile_path)//"C"//trim(adjustl(res_ch))//"_grid.tile7.nc" + + PRINT*, 'Open cubed sphere grid spec file ', trim(gridfile) + + stat = nf90_open(trim(gridfile), NF90_NOWRITE, ncid) + CALL nc_opchk(stat,'nf90_open') + + stat = nf90_inq_dimid(ncid,'nxp',dimid) + CALL nc_opchk(stat,'nf90_inq_dimid: nxp') + stat = nf90_inquire_dimension(ncid,dimid,dimname,len=nxp) + CALL nc_opchk(stat,'nf90_inquire_dimension: nxp') + + stat = nf90_inq_dimid(ncid,'nyp',dimid) + CALL nc_opchk(stat,'nf90_inq_dimid: nyp') + stat = nf90_inquire_dimension(ncid,dimid,dimname,len=nyp) + CALL nc_opchk(stat,'nf90_inquire_dimension: nyp') + +! sidex_sz = nxp-2*halo_depth +! sidey_sz = nyp-2*halo_depth + sidex_sz = nxp + sidey_sz = nyp + tile_sz = sidex_sz*sidey_sz + ALLOCATE(lat(tile_sz), lon(tile_sz)) + +! x_start = halo_depth+1; y_start = halo_depth+1 + x_start = 1; y_start = 1 + stat = nf90_inq_varid(ncid,'x',varid) + CALL nc_opchk(stat,'nf90_inq_lon') + stat = nf90_get_var(ncid,varid,lon,start=(/x_start,y_start/),count=(/sidex_sz,sidey_sz/)) + CALL nc_opchk(stat,'nf90_get_var_lon') + stat = nf90_inq_varid(ncid,'y',varid) + CALL nc_opchk(stat,'nf90_inq_lat') + stat = nf90_get_var(ncid,varid,lat,start=(/x_start,y_start/),count=(/sidex_sz,sidey_sz/)) + CALL nc_opchk(stat,'nf90_get_var_lat') + stat = nf90_close(ncid) + CALL nc_opchk(stat,'nf90_close') + + ALLOCATE(grid(tile_sz,2)) + grid(1:tile_sz,1) = lat(1:tile_sz) + grid(1:tile_sz,2) = lon(1:tile_sz) + + res_x = sidex_sz/2; res_y = sidey_sz/2 + DEALLOCATE(lat,lon) + +END SUBROUTINE read_cubed_sphere_reg_grid + +SUBROUTINE read_lakedata(lakedata_path,lake_stat,lake_dpth,nlat,nlon) + CHARACTER(len=256), INTENT(IN) :: lakedata_path + INTEGER*1, INTENT(OUT) :: lake_stat(:) + INTEGER*2, INTENT(OUT) :: lake_dpth(:) + INTEGER, INTENT(IN) :: nlat, nlon + + CHARACTER(len=256) lakefile + INTEGER :: data_sz, i + + data_sz = nlon*nlat + + ! read in 30 arc seconds lake data + lakefile = trim(lakedata_path) // "GlobalLakeStatus.dat" + OPEN(10, file=lakefile, form='unformatted', access='direct', recl=data_sz*1) + READ(10,rec=1) lake_stat + CLOSE(10) + + lakefile = trim(lakedata_path) // "GlobalLakeDepth.dat" + OPEN(10, file=lakefile, form='unformatted', access='direct', recl=data_sz*2) + READ(10,rec=1) lake_dpth + CLOSE(10) + +END SUBROUTINE read_lakedata + + +SUBROUTINE write_lakedata_to_orodata(cs_res, cs_lakestat, cs_lakedpth) + USE netcdf + INTEGER, INTENT(IN) :: cs_res + REAL, INTENT(IN) :: cs_lakestat(:) + REAL, INTENT(IN) :: cs_lakedpth(:) + + INTEGER :: tile_sz, tile_num + INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2) + INTEGER :: lake_frac_id, lake_depth_id + INTEGER :: land_frac_id, slmsk_id, inland_id, geolon_id, geolat_id + CHARACTER(len=256) :: filename,string + CHARACTER(len=1) :: ich + CHARACTER(len=4) res_ch + REAL :: lake_frac(cs_res*cs_res),lake_depth(cs_res*cs_res) + REAL :: geolon(cs_res*cs_res),geolat(cs_res*cs_res) + REAL :: land_frac(cs_res*cs_res),slmsk(cs_res*cs_res),inland(cs_res*cs_res) + real, parameter :: epsil=1.e-6 ! numerical min for lake_frac/land_frac + real :: land_cutoff=1.e-4 ! land_frac=0 if it is < land_cutoff + + INTEGER :: i, j + +! include "netcdf.inc" + tile_sz = cs_res*cs_res + + WRITE(res_ch,'(I4)') cs_res + DO tile_num = tile_beg, tile_end + WRITE(ich, '(I1)') tile_num +! filename = "C" // trim(adjustl(res_ch)) // "_oro_data.tile" // ich // ".nc" +! filename = "oro_data.tile" // ich // ".nc" + filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc" + print *,'Read, update, and write ',trim(filename) + stat = nf90_open(filename, NF90_WRITE, ncid) + CALL nc_opchk(stat, "nf90_open oro_data.nc") + stat = nf90_inq_dimid(ncid, "lon", x_dimid) + CALL nc_opchk(stat, "nf90_inq_dim: x") + stat = nf90_inq_dimid(ncid, "lat", y_dimid) + CALL nc_opchk(stat, "nf90_inq_dim: y") +! dimids = (/ y_dimid, x_dimid /) +! original orodata netcdf file uses (y, x) order, so we made change to match it. + dimids = (/ x_dimid, y_dimid /) + stat = nf90_inq_varid(ncid, "land_frac", land_frac_id) + CALL nc_opchk(stat, "nf90_inq_varid: land_frac") + stat = nf90_inq_varid(ncid, "slmsk", slmsk_id) + CALL nc_opchk(stat, "nf90_inq_varid: slmsk") +! define 2 new variables + stat = nf90_redef(ncid) + CALL nc_opchk(stat, "nf90_redef") + stat = nf90_def_var(ncid,"lake_frac",NF90_FLOAT,dimids,lake_frac_id) + CALL nc_opchk(stat, "nf90_def_var: lake_frac") +#ifdef ADD_ATT_FOR_NEW_VAR + stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates") + stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name") + stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") + write(string,'(a,f5.2)') 'based on GLDBv2 (Choulga et al. 2014); missing lakes & + added based on land_frac in this dataset; lake_frac cutoff is',lake_cutoff + stat = nf90_put_att(ncid, lake_frac_id,'description',trim(string)) + CALL nc_opchk(stat, "nf90_put_att: lake_frac:description") +#endif + stat = nf90_def_var(ncid,"lake_depth",NF90_FLOAT,dimids,lake_depth_id) + CALL nc_opchk(stat, "nf90_def_var: lake_depth") +#ifdef ADD_ATT_FOR_NEW_VAR + stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates") + stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") + stat = nf90_put_att(ncid, lake_depth_id,'unit','meter') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") + stat = nf90_put_att(ncid, lake_depth_id,'description', & + 'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m & + (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") +#endif + + write(string,'(a,es8.1)') 'land_frac and lake_frac are adjusted such that their sum is 1 at points where inland=1; land_frac cutoff is',land_cutoff + stat = nf90_put_att(ncid, land_frac_id,'description',trim(string)) + CALL nc_opchk(stat, "nf90_put_att: land_frac:description") + + write(string,'(a)') 'slmsk = nint(land_frac)' + stat = nf90_put_att(ncid, slmsk_id,'description',trim(string)) + CALL nc_opchk(stat, "nf90_put_att: slmsk:description") + + stat = nf90_enddef(ncid) + CALL nc_opchk(stat, "nf90_enddef") + +! read in geolon and geolat and 2 variables from orog data file + stat = nf90_inq_varid(ncid, "geolon", geolon_id) + CALL nc_opchk(stat, "nf90_inq_varid: geolon") + stat = nf90_inq_varid(ncid, "geolat", geolat_id) + CALL nc_opchk(stat, "nf90_inq_varid: geolat") + stat = nf90_inq_varid(ncid, "land_frac", land_frac_id) + CALL nc_opchk(stat, "nf90_inq_varid: land_frac") + stat = nf90_inq_varid(ncid, "slmsk", slmsk_id) + CALL nc_opchk(stat, "nf90_inq_varid: slmsk") + stat = nf90_inq_varid(ncid, "inland", inland_id) + CALL nc_opchk(stat, "nf90_inq_varid: inland") + stat = nf90_get_var(ncid, geolon_id, geolon, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_get_var: geolon") + stat = nf90_get_var(ncid, geolat_id, geolat, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_get_var: geolat") + stat = nf90_get_var(ncid, land_frac_id, land_frac, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_get_var: land_frac") + stat = nf90_get_var(ncid, slmsk_id, slmsk, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_get_var: slmsk") + stat = nf90_get_var(ncid, inland_id, inland, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_get_var: inland") + + lake_frac (:) = cs_lakestat ((tile_num-1)*tile_sz+1:tile_num*tile_sz) + lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz) + +! add Caspian Sea and Aral Sea to lake_frac and lake_depth + IF (tile_num == 2 .or. tile_num == 3) THEN + DO i = 1, tile_sz + IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN + IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. & + geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN + lake_frac(i) = 1.-land_frac(i) + lake_depth(i) = 211.0 + ENDIF + IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. & + geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN + lake_frac(i) = 1.-land_frac(i) + lake_depth(i) = 10.0 + ENDIF + ENDIF + ENDDO + ENDIF + +! adjust land_frac and lake_frac, and make sure land_frac+lake_frac=1 at inland points + DO i = 1, tile_sz + if (lake_frac(i) >= lake_cutoff) then ! non-zero lake_frac dominates over land_frac + land_frac(i) = max(0., min(1., 1.-lake_frac(i))) + elseif (inland(i) == 1.) then ! land_frac dominates over lake_frac at inland points + lake_frac(i) = max(0., min(1., 1.-land_frac(i))) + end if + +! epsil is "numerical" nonzero min for lake_frac/land_frac +! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac + if (min(lake_cutoff,land_cutoff) < epsil) then + print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...' + lake_cutoff=max(epsil,lake_cutoff) + land_cutoff=max(epsil,land_cutoff) + end if + + if (lake_frac(i)< lake_cutoff) lake_frac(i)=0. + if (lake_frac(i)>1.-lake_cutoff) lake_frac(i)=1. + if (inland(i) == 1.) land_frac(i) = 1.-lake_frac(i) + if (land_frac(i)< land_cutoff) land_frac(i)=0. + if (land_frac(i)>1.-land_cutoff) land_frac(i)=1. + + if (lake_frac(i) < lake_cutoff) then + lake_depth(i)=0. + elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then + lake_depth(i)=10. + end if + slmsk(i) = nint(land_frac(i)) + ENDDO +! write 2 new variables + stat = nf90_put_var(ncid, lake_frac_id, lake_frac, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_put_var: lake_frac") + + stat = nf90_put_var(ncid, lake_depth_id, lake_depth, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_put_var: lake_depth") + +! write back 2 modified variables: land_frac and slmsk + stat = nf90_put_var(ncid, land_frac_id, land_frac, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_put_var: land_frac") + + stat = nf90_put_var(ncid, slmsk_id, slmsk, & + start = (/ 1, 1 /), count = (/ cs_res, cs_res /) ) + CALL nc_opchk(stat, "nf90_put_var: slmsk") + + stat = nf90_close(ncid) + CALL nc_opchk(stat, "nf90_close") + ENDDO + +END SUBROUTINE write_lakedata_to_orodata + +SUBROUTINE write_reg_lakedata_to_orodata(cs_res, tile_x_dim, tile_y_dim, cs_lakestat, cs_lakedpth) + USE netcdf + INTEGER, INTENT(IN) :: cs_res, tile_x_dim, tile_y_dim + REAL, INTENT(IN) :: cs_lakestat(:) + REAL, INTENT(IN) :: cs_lakedpth(:) + + INTEGER :: tile_sz, tile_num + INTEGER :: stat, ncid, x_dimid, y_dimid, varid, dimids(2) + INTEGER :: lake_frac_id, lake_depth_id + INTEGER :: land_frac_id, slmsk_id, geolon_id, geolat_id + CHARACTER(len=256) :: filename,string + CHARACTER(len=1) :: ich + CHARACTER(len=4) res_ch + + REAL, ALLOCATABLE :: lake_frac(:), lake_depth(:), geolon(:), geolat(:) + REAL, ALLOCATABLE :: land_frac(:), slmsk(:) + + real, parameter :: epsil=1.e-6 ! numerical min for lake_frac/land_frac + real :: land_cutoff=1.e-6 ! land_frac=0 if it is < land_cutoff + + INTEGER :: i, j, var_id + +! include "netcdf.inc" + tile_sz = tile_x_dim*tile_y_dim + + ALLOCATE(lake_frac(tile_sz), lake_depth(tile_sz)) + ALLOCATE(geolon(tile_sz), geolat(tile_sz)) + ALLOCATE(land_frac(tile_sz), slmsk(tile_sz)) + + WRITE(res_ch,'(I4)') cs_res + tile_num = 7 + WRITE(ich, '(I1)') tile_num +! filename = "C" // trim(adjustl(res_ch)) // "_oro_data.tile" // ich // ".halo4.nc" + filename = "oro.C" // trim(adjustl(res_ch)) // ".tile" // ich // ".nc" + PRINT*, 'Open and write regional orography data netCDF file ', trim(filename) + stat = nf90_open(filename, NF90_WRITE, ncid) + CALL nc_opchk(stat, "nf90_open oro_data.nc") + stat = nf90_inq_dimid(ncid, "lon", x_dimid) + CALL nc_opchk(stat, "nf90_inq_dim: x") + stat = nf90_inq_dimid(ncid, "lat", y_dimid) + CALL nc_opchk(stat, "nf90_inq_dim: y") + dimids = (/ x_dimid, y_dimid /) + + stat = nf90_inq_varid(ncid, "land_frac", land_frac_id) + CALL nc_opchk(stat, "nf90_inq_varid: land_frac") + stat = nf90_inq_varid(ncid, "slmsk", slmsk_id) + CALL nc_opchk(stat, "nf90_inq_varid: slmsk") + +! define 2 new variables, lake_frac and lake_depth + stat = nf90_redef(ncid) + CALL nc_opchk(stat, "nf90_redef") + + IF (nf90_inq_varid(ncid, "lake_frac",lake_frac_id) /= 0) THEN + stat = nf90_def_var(ncid,"lake_frac",NF90_FLOAT,dimids,lake_frac_id) + CALL nc_opchk(stat, "nf90_def_var: lake_frac") +#ifdef ADD_ATT_FOR_NEW_VAR + stat = nf90_put_att(ncid, lake_frac_id,'coordinates','geolon geolat') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:coordinates") + stat = nf90_put_att(ncid, lake_frac_id,'long_name','lake fraction') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:long_name") + stat = nf90_put_att(ncid, lake_frac_id,'unit','fraction') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:unit") + stat = nf90_put_att(ncid, lake_frac_id,'description', & + 'based on GLDBv2 (Choulga et al. 2014); missing lakes & + added based on land_frac in this dataset.') + CALL nc_opchk(stat, "nf90_put_att: lake_frac:description") +#endif + ENDIF + IF (nf90_inq_varid(ncid, "lake_depth",lake_depth_id) /= 0) THEN + stat = nf90_def_var(ncid,"lake_depth",NF90_FLOAT,dimids,lake_depth_id) + CALL nc_opchk(stat, "nf90_def_var: lake_depth") +#ifdef ADD_ATT_FOR_NEW_VAR + stat = nf90_put_att(ncid, lake_depth_id,'coordinates','geolon geolat') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:coordinates") + stat = nf90_put_att(ncid, lake_depth_id,'long_name','lake depth') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") + stat = nf90_put_att(ncid, lake_depth_id,'unit','meter') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:long_name") + stat = nf90_put_att(ncid, lake_depth_id,'description', & + 'based on GLDBv2 (Choulga et al. 2014); missing depth set to 10m & + (except to 211m in Caspian Sea); spurious large pos. depths are left unchanged.') + CALL nc_opchk(stat, "nf90_put_att: lake_depth:description") +#endif + ENDIF + write(string,'(a,es8.1)') 'land_frac is adjusted to 1-lake_frac where lake_frac>0 but left unchanged where lake_frac=0. This could lead to land_frac+lake_frac<1 at some inland points; land_frac cutoff is',land_cutoff + stat = nf90_put_att(ncid, land_frac_id,'description',trim(string)) + CALL nc_opchk(stat, "nf90_put_att: land_frac:description") + + write(string,'(a)') 'slmsk = nint(land_frac)' + stat = nf90_put_att(ncid, slmsk_id,'description',trim(string)) + CALL nc_opchk(stat, "nf90_put_att: slmsk:description") + + stat = nf90_enddef(ncid) + CALL nc_opchk(stat, "nf90_enddef") + +! read in geolon and geolat and 2 variables from orog data file + stat = nf90_inq_varid(ncid, "geolon", geolon_id) + CALL nc_opchk(stat, "nf90_inq_varid: geolon") + stat = nf90_inq_varid(ncid, "geolat", geolat_id) + CALL nc_opchk(stat, "nf90_inq_varid: geolat") + stat = nf90_inq_varid(ncid, "land_frac", land_frac_id) + CALL nc_opchk(stat, "nf90_inq_varid: land_frac") + stat = nf90_inq_varid(ncid, "slmsk", slmsk_id) + CALL nc_opchk(stat, "nf90_inq_varid: slmsk") + + stat = nf90_get_var(ncid, geolon_id, geolon, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_get_var: geolon") + stat = nf90_get_var(ncid, geolat_id, geolat, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_get_var: geolat") + stat = nf90_get_var(ncid, land_frac_id, land_frac, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_get_var: land_frac") + stat = nf90_get_var(ncid, slmsk_id, slmsk, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_get_var: slmsk") + + tile_num = 1 + lake_frac(:) = cs_lakestat((tile_num-1)*tile_sz+1:tile_num*tile_sz) + lake_depth(:) = cs_lakedepth((tile_num-1)*tile_sz+1:tile_num*tile_sz) + +! add Caspian Sea and Aral Sea to lake_frac and lake_depth + IF (tile_num == 2 .or. tile_num == 3) THEN + DO i = 1, tile_sz + IF (land_frac(i) < 0.9 .AND. lake_frac(i) < 0.1) THEN + IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. & + geolon(i) > 45.0 .AND. geolon(i) <= 55.0) THEN + lake_frac(i) = 1.-land_frac(i) + lake_depth(i) = 211.0 + ENDIF + IF (geolat(i) > 35.0 .AND. geolat(i) <= 50.0 .AND. & + geolon(i) > 57.0 .AND. geolon(i) <= 63.0) THEN + lake_frac(i) = 1.-land_frac(i) + lake_depth(i) = 10.0 + ENDIF + ENDIF + ENDDO + ENDIF + + DO i = 1, tile_sz +! epsil is "numerical" nonzero min for lake_frac/land_frac +! lake_cutoff/land_cutoff is practical min for lake_frac/land_frac + if (min(lake_cutoff,land_cutoff) < epsil) then + print *,'lake_cutoff/land_cutoff cannot be smaller than epsil, reset...' + lake_cutoff=max(epsil,lake_cutoff) + land_cutoff=max(epsil,land_cutoff) + end if + + if (lake_frac(i)< lake_cutoff) lake_frac(i)=0. + if (lake_frac(i)>1.-lake_cutoff) lake_frac(i)=1. + if (land_frac(i)< land_cutoff) land_frac(i)=0. + if (land_frac(i)>1.-land_cutoff) land_frac(i)=1. + + if (lake_frac(i) >= lake_cutoff) then ! non-zero lake_frac dominates over land_frac + land_frac(i) = max(0., min(1., 1.-lake_frac(i))) + end if + + if (lake_frac(i) < lake_cutoff) then + lake_depth(i)=0. + elseif (lake_frac(i) >= lake_cutoff .and. lake_depth(i)==0.) then + lake_depth(i)=10. + end if + slmsk(i) = nint(land_frac(i)) + ENDDO + +! write 2 new variables + stat = nf90_put_var(ncid, lake_frac_id, lake_frac, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_put_var: lake_frac") + + stat = nf90_put_var(ncid, lake_depth_id, lake_depth, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_put_var: lake_depth") + +! write back 2 modified variables: land_frac and slmsk + stat = nf90_put_var(ncid, land_frac_id, land_frac, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_put_var: land_frac") + + stat = nf90_put_var(ncid, slmsk_id, slmsk, & + start = (/ 1, 1 /), count = (/ tile_x_dim, tile_y_dim /) ) + CALL nc_opchk(stat, "nf90_put_var: slmsk") + + stat = nf90_close(ncid) + CALL nc_opchk(stat, "nf90_close") + +END SUBROUTINE write_reg_lakedata_to_orodata + +SUBROUTINE nc_opchk(stat,opname) + USE netcdf + IMPLICIT NONE + INTEGER stat + CHARACTER(len=*) opname + CHARACTER(64) msg + + IF (stat .NE.0) THEN + msg = trim(opname) // ' Error, status code and message:' + PRINT*,trim(msg), stat, nf90_strerror(stat) + STOP + END IF + +END SUBROUTINE nc_opchk + +END PROGRAM lake_frac diff --git a/sorc/orog.fd/.gitignore b/sorc/orog_mask_tools.fd/orog.fd/.gitignore similarity index 100% rename from sorc/orog.fd/.gitignore rename to sorc/orog_mask_tools.fd/orog.fd/.gitignore diff --git a/sorc/orog.fd/CMakeLists.txt b/sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt similarity index 100% rename from sorc/orog.fd/CMakeLists.txt rename to sorc/orog_mask_tools.fd/orog.fd/CMakeLists.txt diff --git a/sorc/orog.fd/README b/sorc/orog_mask_tools.fd/orog.fd/README similarity index 100% rename from sorc/orog.fd/README rename to sorc/orog_mask_tools.fd/orog.fd/README diff --git a/sorc/orog.fd/README_OC b/sorc/orog_mask_tools.fd/orog.fd/README_OC similarity index 100% rename from sorc/orog.fd/README_OC rename to sorc/orog_mask_tools.fd/orog.fd/README_OC diff --git a/sorc/orog.fd/machine.h b/sorc/orog_mask_tools.fd/orog.fd/machine.h similarity index 100% rename from sorc/orog.fd/machine.h rename to sorc/orog_mask_tools.fd/orog.fd/machine.h diff --git a/sorc/orog.fd/mtnlm7_oclsm.f b/sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.f similarity index 100% rename from sorc/orog.fd/mtnlm7_oclsm.f rename to sorc/orog_mask_tools.fd/orog.fd/mtnlm7_oclsm.f diff --git a/sorc/orog.fd/netcdf_io.F90 b/sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90 similarity index 100% rename from sorc/orog.fd/netcdf_io.F90 rename to sorc/orog_mask_tools.fd/orog.fd/netcdf_io.F90 diff --git a/sorc/orog.fd/resevod.h b/sorc/orog_mask_tools.fd/orog.fd/resevod.h similarity index 100% rename from sorc/orog.fd/resevod.h rename to sorc/orog_mask_tools.fd/orog.fd/resevod.h diff --git a/sorc/sfc_climo_gen.fd/model_grid.F90 b/sorc/sfc_climo_gen.fd/model_grid.F90 index b123b6bfa..bd0c0b5f0 100644 --- a/sorc/sfc_climo_gen.fd/model_grid.F90 +++ b/sorc/sfc_climo_gen.fd/model_grid.F90 @@ -334,7 +334,7 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim) real(esmf_kind_r4), intent(out) :: lat2d(idim,jdim) real(esmf_kind_r4), intent(out) :: lon2d(idim,jdim) - integer :: error, lat, lon + integer :: error, lat, lon, i, j integer :: ncid, id_dim, id_var real(kind=4), allocatable :: dummy(:,:) @@ -365,12 +365,37 @@ subroutine get_model_info(orog_file, mask, lat2d, lon2d, idim, jdim) allocate(dummy(idim,jdim)) - print*,"- READ LAND MASK" - error=nf90_inq_varid(ncid, 'slmsk', id_var) - call netcdf_err(error, "READING SLMSK ID") - error=nf90_get_var(ncid, id_var, dummy) - call netcdf_err(error, "READING SLMSK") - mask = nint(dummy) +!----------------------------------------------------------------------- +! If the lake maker was used, there will be a 'lake_frac' record. +! In that case, land/non-land is determined by 'land_frac'. +! +! If the lake maker was not used, use 'slmsk', which is defined +! as the nint(land_frac). +!----------------------------------------------------------------------- + + error=nf90_inq_varid(ncid, 'lake_frac', id_var) + if (error /= 0) then + print*,"- READ LAND MASK (SLMSK)" + error=nf90_inq_varid(ncid, 'slmsk', id_var) + call netcdf_err(error, "READING SLMSK ID") + error=nf90_get_var(ncid, id_var, dummy) + call netcdf_err(error, "READING SLMSK") + mask = nint(dummy) + else + print*,"- READ LAND FRACTION" + error=nf90_inq_varid(ncid, 'land_frac', id_var) + call netcdf_err(error, "READING LAND_FRAC ID") + error=nf90_get_var(ncid, id_var, dummy) + call netcdf_err(error, "READING LAND_FRAC") + mask = 0 + do j = 1, lat + do i = 1, lon + if (dummy(i,j) > 0.0) then + mask(i,j) = 1 + endif + enddo + enddo + endif print*,"- READ LATITUDE" error=nf90_inq_varid(ncid, 'geolat', id_var) diff --git a/ush/fv3gfs_driver_grid.sh b/ush/fv3gfs_driver_grid.sh index dc37d28ff..202a8e6a6 100755 --- a/ush/fv3gfs_driver_grid.sh +++ b/ush/fv3gfs_driver_grid.sh @@ -42,8 +42,12 @@ export machine=${machine:?} #---------------------------------------------------------------------------------- export res=${res:-96} # resolution of tile: 48, 96, 128, 192, 384, 768, 1152, 3072 -export gtype=${gtype:-uniform} # grid type: uniform, stretch, nest, regional_gfdl - # or regional_esg +export gtype=${gtype:-uniform} # grid type: uniform, stretch, nest, regional_gfdl, + # or regional_esg + +export add_lake=${add_lake:-false} # add lake fraction and depth. uniform only. +export lake_cutoff=${lake_cutoff:-0.20} # lake fractions less than lake_cutoff + # are ignored. if [ $gtype = uniform ]; then echo "Creating global uniform grid" @@ -88,6 +92,15 @@ else exit 9 fi +if [ $add_lake = true ]; then + if [ $gtype != uniform ]; then + set +x + echo "Adding lake data to orography data is only available" + echo "for uniform grids." + exit 1 + fi +fi + export TMPDIR=${TMPDIR:?} export out_dir=${out_dir:?} @@ -194,6 +207,14 @@ if [ $gtype = uniform ] || [ $gtype = stretch ] || [ $gtype = nest ]; then done fi + if [ $add_lake = true ]; then + $script_dir/fv3gfs_make_lake.sh + err=$? + if [ $err != 0 ]; then + exit $err + fi + fi + set +x echo "End uniform orography generation at `date`" set -x diff --git a/ush/fv3gfs_make_lake.sh b/ush/fv3gfs_make_lake.sh new file mode 100755 index 000000000..ad0d81f14 --- /dev/null +++ b/ush/fv3gfs_make_lake.sh @@ -0,0 +1,99 @@ +#!/bin/ksh + +echo +echo "CREATE AND ADD INLAND, LAKE_STATUS, AND LAKE_DEPTH TO THE OROGRAPHY FILES" +echo + +set -ax + +outdir=$orog_dir +indir=$topo + +if [ $gtype != uniform ]; then + echo "lake_frac has not been implemented for FV3 stand-alone Regional (SAR) model" + exit 1 +fi + +exe_add_lake=$exec_dir/lakefrac +if [ ! -s $exe_add_lake ]; then + echo "executable to add lake does not exist" + exit 1 +fi + +exe_inland=$exec_dir/inland +if [ ! -s $exe_inland ]; then + echo "executable to create inland mask does not exist" + exit 1 +fi + +workdir=$TMPDIR/C${res}/orog/tiles +if [ ! -s $workdir ]; then mkdir -p $workdir ;fi + +# Make lake data - +# Create and add inland, lake_status, and lake_depth to the orography files + +echo "workdir = $workdir" +echo "outdir = $outdir" +echo "indir = $indir" + +cd $workdir + +# link all required files to the current work directory + +tile_beg=1 +tile_end=6 +tile=$tile_beg +while [ $tile -le $tile_end ]; do + outfile=oro.C${res}.tile${tile}.nc + ln -s $outdir/$outfile . + outgrid="C${res}_grid.tile${tile}.nc" + ln -s $grid_dir/$outgrid . + tile=$(( $tile + 1 )) +done + +# create inland mask and save it to the orography files + +cutoff=0.99 +rd=7 +$APRUN $exe_inland $res $cutoff $rd +err=$? +if [ $err != 0 ]; then + set +x + echo ERROR CREATING INLAND MASK + exit $err +fi + +# create lake data for FV3 grid and save it to the orography files + +if [ $machine = WCOSS_C ]; then + touch ./lake.txt + tile=$tile_beg + while [ $tile -le $tile_end ]; do + echo "$exe_add_lake ${tile} ${res} ${indir} ${lake_cutoff}" >> ./lake.txt + tile=$(( $tile + 1 )) + done + aprun -j 1 -n 6 -N 6 -d 1 -cc depth cfp ./lake.txt + err=$? + if [ $err != 0 ]; then + set +x + echo ERROR CREATING LAKE FRACTION + exit $err + fi + rm ./lake.txt +else + tile=$tile_beg + while [ $tile -le $tile_end ]; do + outfile=oro.C${res}.tile${tile}.nc + $APRUN $exe_add_lake ${tile} ${res} ${indir} ${lake_cutoff} + err=$? + if [ $err != 0 ]; then + set +x + echo ERROR CREATING LAKE FRACTION FOR TILE $tile + exit $err + fi + echo "lake fraction is added to $outfile" + tile=$(( $tile + 1 )) + done +fi + +exit 0