diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt index 5fa2e6df1..e20ff706b 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/CMakeLists.txt @@ -2,9 +2,8 @@ esma_set_this(OVERRIDE raster) set (srcs date_time_util.F90 leap_year.F90 -easeV1_conv.F90 +EASE_conv.F90 mod_process_hres_data.F90 -easeV2_conv.F90 rasterize.F90 read_riveroutlet.F90 CubedSphere_GridMod.F90 @@ -47,7 +46,6 @@ ecbuild_add_executable (TARGET mkMITAquaRaster.x SOURCES mkMITAquaRaster.F90 LIB ecbuild_add_executable (TARGET mkMOMAquaRaster.x SOURCES mkMOMAquaRaster.F90 LIBS MAPL ${this}) ecbuild_add_executable (TARGET FillMomGrid.x SOURCES FillMomGrid.F90 LIBS MAPL ${this}) ecbuild_add_executable (TARGET mk_runofftbl.x SOURCES mk_runofftbl.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkSMAPTilesPara.x SOURCES mkSMAPTilesPara.F90 LIBS MAPL ${this}) -ecbuild_add_executable (TARGET mkSMAPTilesPara_v2.x SOURCES mkSMAPTilesPara_v2.F90 LIBS MAPL ${this}) +ecbuild_add_executable (TARGET mkEASETilesParam.x SOURCES mkEASETilesParam.F90 LIBS MAPL ${this}) install(PROGRAMS make_bcs clsm_plots.pro plot_curves.pro create_README.csh plot_curves.csh DESTINATION bin) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/EASE_conv.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/EASE_conv.F90 new file mode 100644 index 000000000..2060d77d5 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/EASE_conv.F90 @@ -0,0 +1,786 @@ + +module EASE_conv + + ! ===================================================================================== + ! + ! EASE_conv.F90 - FORTRAN routines for conversion of Equal-Area Scalable Earth (EASE) + ! grid coordinates (lat/lon <--> row/col indices) + ! + ! Implemented for global cylindrical ('Mxx') EASE grids only. + ! + ! Works for EASE[v1] and EASEv2 grids. + ! + ! ------------------------------------------------------------------------------------- + ! + ! CHANGELOG (easeV1_conv.F90): + ! ============================ + ! + ! easeV1_conv.F90 - FORTRAN routines for conversion of azimuthal + ! equal area and equal area cylindrical grid coordinates + ! + ! 30-Jan-1992 H.Maybee + ! 20-Mar-1992 Ken Knowles 303-492-0644 knowles@kryos.colorado.edu + ! 16-Dec-1993 MJ Brodzik 303-492-8263 brodzik@jokull.colorado.edu + ! Copied from nsmconv.f, changed resolutions from + ! 40-20-10 km to 25-12.5 km + ! 21-Dec-1993 MJ Brodzik 303-492-8263 brodzik@jokull.colorado.edu + ! Fixed sign of Southern latitudes in ease_inverse. + ! 12-Sep-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu + ! Changed grid cell size. Changed "c","f" to "l","h" + ! 25-Oct-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu + ! Changed row size from 587 to 586 for Mercator projection + ! 11-May-2011 reichle: Changed "smap" to "easeV1". + ! Added SSM/I and AMSR-E "M25" grid. + ! So far ONLY for cylindrical grids. + ! Converted from *.f to *.F90 module + ! + ! $Log$ + ! Revision 1.1 2011-05-11 21:58:46 rreichle + ! Adding utilities to map between EASE grids and lat/lon coordinates. + ! + ! Revision 1.3 1994/11/01 23:40:43 brodzik + ! Replaced all references to 'ease' with 'smap' + ! Replaced all references to 'smap' with 'easeV1' -- reichle + ! + ! + ! CHANGELOG (easeV2_conv.F90): + ! ============================ + ! + ! easeV2_conv.F90 - FORTRAN routines for converting grid coordinates + ! (latitude/longitude <--> row/column indices) + ! of the Equal Area Scalable Earth, version 2 (EASEv2) grid + ! + ! ***** ONLY cylindrical ('M') projection implemented ***** + ! + ! Ported from Steven Chan's matlab code (smapease2inverse.m, + ! smapease2forward.m), which has been ported from NSIDC's IDL code + ! (wgs84_convert.pro, wgs84_inverse.pro) available from + ! ftp://sidads.colorado.edu/pub/tools/easegrid/geolocation_tools/ + ! + ! Official references: + ! doi:10.3390/ijgi1010032 + ! doi:10.3390/ijgi3031154 -- correction of M25 "map_scale_m" parameters! + ! + ! 04-Apr-2013 - reichle + ! 11-Sep-2018 - reichle, mgirotto -- added 'M25' grid parameters + ! + ! + ! CHANGELOG (EASE_conv.F90): + ! ========================== + ! + ! 2022-09-13, wjiang+reichle: + ! merged easeV1_conv.F90 and easeV2_conv.F90 into EASE_conv.F90 + ! - using different values for PI in easeV1 and easeV2 calcs as in old easeV[x]_conv.F90 modules; + ! in contrast, LDAS_EASE_conv.F90 in GEOSldas used only a single value for PI. + ! - bug fix in easeV2_get_params() for EASEv2/M25 (to compute s0, divide by 2.0 not by integer 2) + ! + ! + ! ========================================================================== + + implicit none + + private + + public :: ease_convert + public :: ease_inverse + public :: ease_extent + + ! ======================================================================= + ! + ! EASEv1 global constants + + ! ***NEVER*** change these constants to GEOS MAPL constants!!!! + + ! radius of the earth (km), authalic sphere based on International datum + + real*8, parameter :: easeV1_RE_km = 6371.228 + + ! scale factor for standard paralles at +/-30.00 degrees + + real*8, parameter :: easeV1_COS_PHI1 = .866025403 + + real*8, parameter :: easeV1_PI = 3.141592653589793 + + ! ======================================================================= + ! + ! EASEv2 global constants + + ! ***NEVER*** change these constants to GEOS MAPL constants!!!! + + ! radius of the earth (m) and map eccentricity + + real*8, parameter :: map_equatorial_radius_m = 6378137.0 + + real*8, parameter :: map_eccentricity = 0.081819190843 + + real*8, parameter :: easeV2_PI = 3.14159265358979323846 + + real*8, parameter :: e2 = map_eccentricity * map_eccentricity + real*8, parameter :: e4 = e2 * e2 + real*8, parameter :: e6 = e2 * e4 + + real*8, parameter :: epsilon = 1.e-6 + + real*8, parameter :: map_reference_longitude = 0.0 ! 'M', 'N', 'S' + + ! constants for 'N' and 'S' (azimuthal) projections + + real*8, parameter :: N_map_reference_latitude = 90.0 + real*8, parameter :: S_map_reference_latitude = -90.0 + + ! constants for 'M' (cylindrical) projection + + real*8, parameter :: M_map_reference_latitude = 0.0 + real*8, parameter :: M_map_second_reference_latitude = 30.0 + + real*8, parameter :: M_sin_phi1 = sin(M_map_second_reference_latitude*easeV2_PI/180.) + real*8, parameter :: M_cos_phi1 = cos(M_map_second_reference_latitude*easeV2_PI/180.) + + real*8, parameter :: M_kz = M_cos_phi1/sqrt(1.0-e2*M_sin_phi1*M_sin_phi1) + + +contains + + ! ******************************************************************* + ! + ! GENERIC routines (public interface) + ! + ! ******************************************************************* + ! + ! EASELabel = *EASEv[x]_[p][yy]* (e.g., EASEv2_M09) + ! + ! version: x = { 1, 2 } + ! projection: p = { M } ! only cylindrical ("M") implemented + ! resolution: yy = { 01, 03, 09, 25, 36 } ! 12.5 km not yet implemented + ! + ! Coordinate arguments for ease_convert() and ease_inverse(): + ! + ! | map coords | 0-based index | # grid cells | + ! | | (real numbers!) | | + ! ---------------------------------------------------------- + ! | latitude | s | rows | + ! | longitude | r | cols | + ! + ! Indices are 0-based and run west to east (r) and north to south (s). + ! + ! -------------------------------------------------------------------- + + subroutine ease_convert (EASELabel, lat, lon, r, s) ! note odd/reversed order of (lat,lon) and (r,s) + + character*(*), intent(in) :: EASELabel + real, intent(in) :: lat, lon + real, intent(out) :: r, s ! r = lon index, s = lat index + + character(3) :: grid + + if ( index(EASELabel,'M36') /=0 ) then + grid='M36' + else if (index(EASELabel,'M25') /=0 ) then + grid='M25' + else if (index(EASELabel,'M09') /=0 ) then + grid='M09' + else if (index(EASELabel,'M03') /=0 ) then + grid='M03' + else if (index(EASELabel,'M01') /=0 ) then + grid='M01' + else + print*,"ease_convert(): unknown grid projection and resolution: "//trim(EASELabel)//" STOPPING." + stop + endif + + if( index(EASELabel,'EASEv2') /=0) then + call easeV2_convert(grid,lat,lon,r,s) + else if(index(EASELabel,'EASEv1') /=0) then + call easeV1_convert(grid,lat,lon,r,s) + else + print*,"ease_convert(): unknown grid version: "//trim(EASELabel)//" STOPPING." + stop + endif + + end subroutine ease_convert + + ! ******************************************************************* + + subroutine ease_inverse (EASELabel, r, s, lat, lon) ! note odd/reversed order of (r,s) and (lat,lon) + + ! Note: Get lat/lon of grid cell borders by using fractional indices. + ! E.g., s=-0.5 yields northern grid cell boundary of northernmost grid cells. + + character*(*), intent(in) :: EASELabel + real, intent(in) :: r, s ! r = lon index, s = lat index + real, intent(out) :: lat, lon + + character(3) :: grid + + if ( index(EASELabel,'M36') /=0 ) then + grid='M36' + else if (index(EASELabel,'M25') /=0 ) then + grid='M25' + else if (index(EASELabel,'M09') /=0 ) then + grid='M09' + else if (index(EASELabel,'M03') /=0 ) then + grid='M03' + else if (index(EASELabel,'M01') /=0 ) then + grid='M01' + else + print*,"ease_inverse(): unknown grid projection and resolution: "//trim(EASELabel)//" STOPPING." + stop + endif + + if( index(EASELabel,'EASEv2') /=0) then + call easeV2_inverse(grid,r,s,lat,lon) + else if(index(EASELabel,'EASEv1') /=0) then + call easeV1_inverse(grid,r,s,lat,lon) + else + print*,"ease_inverse(): unknown grid version: "//trim(EASELabel)//" STOPPING." + stop + endif + + end subroutine ease_inverse + + ! ******************************************************************* + + subroutine ease_extent (EASELabel, cols, rows, cell_area, ll_lon, ll_lat, ur_lon, ur_lat) + + ! get commonly used EASE grid parameters + + character*(*), intent(in) :: EASELabel + integer, intent(out) :: cols, rows ! number of grid cells in lon and lat direction, resp. + real, optional, intent(out) :: cell_area ! [m^2] + real, optional, intent(out) :: ll_lon ! lon of grid cell boundary in lower left corner + real, optional, intent(out) :: ll_lat ! lat of grid cell boundary in lower left corner + real, optional, intent(out) :: ur_lon ! lon of grid cell boundary in upper right corner + real, optional, intent(out) :: ur_lat ! lat of grid cell boundary in upper right corner + + ! --------------------------------------------------------------------- + + real*8 :: map_scale_m, CELL_km, r0, s0, Rg + real :: tmplon + character(3) :: grid + + if ( index(EASELabel,'M36') /=0 ) then + grid='M36' + else if (index(EASELabel,'M25') /=0 ) then + grid='M25' + else if (index(EASELabel,'M09') /=0 ) then + grid='M09' + else if (index(EASELabel,'M03') /=0 ) then + grid='M03' + else if (index(EASELabel,'M01') /=0 ) then + grid='M01' + else + print*,"ease_extent(): unknown grid projection and resolution: "//trim(EASELabel)//" STOPPING." + stop + endif + + if( index(EASELabel,'EASEv2') /=0) then + + call easeV2_get_params(grid, map_scale_m, cols, rows, r0, s0) + + if(present(cell_area)) cell_area = map_scale_m**2 + + else if(index(EASELabel,'EASEv1') /=0) then + + call easeV1_get_params(grid, CELL_km, cols, rows, r0, s0, Rg) + + if(present(cell_area)) cell_area = CELL_km**2 * 1000. * 1000. + + else + + print*,"ease_extent(): unknown grid version: "//trim(EASELabel)//" STOPPING." + stop + + endif + + ! get lat/lon of corner grid cells + ! + ! recall that EASE grid indexing is zero-based + + if (present(ll_lat)) call ease_inverse(EASElabel, 0., rows-0.5, ll_lat, tmplon) + if (present(ur_lat)) call ease_inverse(EASElabel, 0., -0.5, ur_lat, tmplon) + + if (present(ll_lon)) ll_lon = -180. + if (present(ur_lon)) ur_lon = 180. + + end subroutine ease_extent + + ! ******************************************************************* + ! + ! EASEv1 routines (private) + ! + ! ******************************************************************* + + subroutine easeV1_convert (grid, lat, lon, r, s) + + ! convert geographic coordinates (spherical earth) to + ! azimuthal equal area or equal area cylindrical grid coordinates + ! + ! status = easeV1_convert (grid, lat, lon, r, s) + ! + ! input : grid - projection name '[M][xx]' + ! where xx = approximate resolution [km] + ! ie xx = "01", "03", "09", "36" (SMAP) + ! or xx = "12", "25" (SSM/I, AMSR-E) + ! lat, lon = geo. coords. (decimal degrees) + ! + ! output: r, s - column, row coordinates + ! + ! result: status = 0 indicates normal successful completion + ! -1 indicates error status (point not on grid) + ! + ! -------------------------------------------------------------------------- + + character*(*), intent(in) :: grid + real, intent(in) :: lat, lon + real, intent(out) :: r, s + + ! local variables + + integer :: cols, rows + real*8 :: Rg, phi, lam, rho, CELL_km, r0, s0 + + real*8, parameter :: PI = easeV1_PI + + ! --------------------------------------------------------------------- + + call easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) + + phi = lat*PI/180. ! convert from degree to radians + lam = lon*PI/180. ! convert from degree to radians + + if (grid(1:1).eq.'N') then + rho = 2 * Rg * sin(PI/4. - phi/2.) + r = r0 + rho * sin(lam) + s = s0 + rho * cos(lam) + + else if (grid(1:1).eq.'S') then + rho = 2 * Rg * cos(PI/4. - phi/2.) + r = r0 + rho * sin(lam) + s = s0 - rho * cos(lam) + + else if (grid(1:1).eq.'M') then + r = r0 + Rg * lam * easeV1_COS_PHI1 + s = s0 - Rg * sin(phi) / easeV1_COS_PHI1 + + endif + + end subroutine easeV1_convert + + ! ******************************************************************* + + subroutine easeV1_inverse (grid, r, s, lat, lon) + + ! convert azimuthal equal area or equal area cylindrical + ! grid coordinates to geographic coordinates (spherical earth) + ! + ! status = easeV1_inverse (grid, r, s, lat, lon) + ! + ! input : grid - projection name '[M][xx]' + ! where xx = approximate resolution [km] + ! ie xx = "01", "03", "09", "36" (SMAP) + ! or xx = "12", "25" (SSM/I, AMSR-E) + ! r, s - column, row coordinates + ! + ! output: lat, lon = geo. coords. (decimal degrees) + ! + ! result: status = 0 indicates normal successful completion + ! -1 indicates error status (point not on grid) + ! + ! -------------------------------------------------------------------------- + + character*(*), intent(in) :: grid + real, intent(in) :: r, s + real, intent(out) :: lat, lon + + ! local variables + + integer :: cols, rows + real*8 :: Rg, phi, lam, rho, CELL_km, r0, s0 + real*8 :: gamma, beta, epsilon, x, y, c + real*8 :: sinphi1, cosphi1 + + real*8, parameter :: PI = easeV1_PI + + ! --------------------------------------------------------------------- + + call easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) + + x = r - r0 + y = -(s - s0) + + if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then + rho = sqrt(x*x + y*y) + if (rho.eq.0.0) then + if (grid(1:1).eq.'N') lat = 90.0 + if (grid(1:1).eq.'S') lat = -90.0 + lon = 0.0 + else + if (grid(1:1).eq.'N') then + sinphi1 = sin(PI/2.) + cosphi1 = cos(PI/2.) + if (y.eq.0.) then + if (r.le.r0) lam = -PI/2. + if (r.gt.r0) lam = PI/2. + else + lam = atan2(x,-y) + endif + else if (grid(1:1).eq.'S') then + sinphi1 = sin(-PI/2.) + cosphi1 = cos(-PI/2.) + if (y.eq.0.) then + if (r.le.r0) lam = -PI/2. + if (r.gt.r0) lam = PI/2. + else + lam = atan2(x,y) + endif + endif + gamma = rho/(2 * Rg) + if (abs(gamma) .gt. 1.) return + c = 2 * asin(gamma) + beta = cos(c) * sinphi1 + y * sin(c) * (cosphi1/rho) + if (abs(beta).gt.1.) return + phi = asin(beta) + lat = phi*180./PI ! convert from radians to degree + lon = lam*180./PI ! convert from radians to degree + endif + + else if (grid(1:1).eq.'M') then + + ! allow .5 cell tolerance in arcsin function + ! so that grid coordinates which are less than .5 cells + ! above 90.00N or below 90.00S are given a lat of 90.00 + + epsilon = 1 + 0.5/Rg + beta = y*easeV1_COS_PHI1/Rg + if (abs(beta).gt.epsilon) return + if (beta.le.-1.) then + phi = -PI/2. + else if (beta.ge.1.) then + phi = PI/2. + else + phi = asin(beta) + endif + lam = x/easeV1_COS_PHI1/Rg + lat = phi*180./PI ! convert from radians to degree + lon = lam*180./PI ! convert from radians to degree + endif + + end subroutine easeV1_inverse + + ! ******************************************************************* + + subroutine easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) + + implicit none + + character*(*), intent(in) :: grid + real*8, intent(out) :: CELL_km, r0, s0, Rg + integer, intent(out) :: cols, rows + + ! -------------------------------------------------------- + ! + ! r0,s0 are defined such that cells at all scales have + ! coincident center points + ! + !c r0 = (cols-1)/2. * scale + !c s0 = (rows-1)/2. * scale + ! + ! -------------------------------------------------------- + + if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then + + print *,'easeV1_get_params(): polar projections not implemented yet' + stop + + else if (grid(1:1).eq.'M') then + + if (grid .eq. 'M36') then ! SMAP 36 km grid + CELL_km = 36.00040279063 ! nominal cell size in kilometers + cols = 963 + rows = 408 + r0 = 481.0 + s0 = 203.5 + + else if (grid .eq. 'M25') then ! SSM/I, AMSR-E 25 km grid + CELL_km = 25.067525 ! nominal cell size in kilometers + cols = 1383 + rows = 586 + r0 = 691.0 + s0 = 292.5 + + else if (grid .eq. 'M09') then ! SMAP 9 km grid + CELL_km = 9.00010069766 ! nominal cell size in kilometers + cols = 3852 + rows = 1632 + r0 = 1925.5 + s0 = 815.5 + + else if (grid .eq. 'M03') then ! SMAP 3 km grid + CELL_km = 3.00003356589 ! nominal cell size in kilometers + cols = 11556 + rows = 4896 + r0 = 5777.5 + s0 = 2447.5 + + else if (grid .eq. 'M01') then ! SMAP 1 km grid + CELL_km = 1.00001118863 ! nominal cell size in kilometers + cols = 34668 + rows = 14688 + r0 = 17333.5 + s0 = 7343.5 + + else + + print *,'easeV1_get_params(): unknown resolution: ',grid + stop + + endif + + else + + print *, 'easeV1_get_params(): unknown projection: ', grid + stop + + endif + + Rg = easeV1_RE_km/CELL_km + + end subroutine easeV1_get_params + + + ! ******************************************************************* + ! + ! EASEv2 routines (private) + ! + ! ******************************************************************* + + subroutine easeV2_convert (grid, lat, lon, col_ind, row_ind) + + ! convert geographic coordinates (spherical earth) to + ! azimuthal equal area or equal area cylindrical grid coordinates + ! + ! *** NOTE order of calling arguments: "lat-lon-lon-lat" *** + ! + ! useage: call easeV2_convert (grid, lat, lon, r, s) + ! + ! input : grid - projection name '[M][xx]' + ! where xx = approximate resolution [km] + ! ie xx = "01", "03", "09", "36" (SMAP) + ! lat, lon = geo. coords. (decimal degrees) + ! + ! output: col_ind, row_ind - column, row coordinates + ! + ! -------------------------------------------------------------------------- + + character*(*), intent(in) :: grid + real, intent(in) :: lat, lon + real, intent(out) :: col_ind, row_ind + + ! local variables + + integer :: cols, rows + real*8 :: dlon, phi, lam, map_scale_m, r0, s0, ms, x, y, sin_phi, q + + real*8, parameter :: PI = easeV2_PI + + ! --------------------------------------------------------------------- + + call easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) + + dlon = lon + + if (abs(map_reference_longitude)>epsilon) then + + dlon = lon - map_reference_longitude + + end if + + if (dlon .lt. -180.0) dlon = dlon + 360.0 + if (dlon .gt. 180.0) dlon = dlon - 360.0 + + phi = lat*PI/180. ! convert from degree to radians + lam = dlon*PI/180. ! convert from degree to radians + + sin_phi = sin(phi) + + ms = map_eccentricity*sin_phi + + q = (1. - e2)* & + ( & + (sin_phi /(1. - e2*sin_phi*sin_phi)) & + - & + .5/map_eccentricity*log((1.-ms)/(1.+ms)) & + ) + + ! note: "qp" only needed for 'N' and 'S' projections + + if (grid(1:1).eq.'M') then + + x = map_equatorial_radius_m*M_kz*lam + + y = (map_equatorial_radius_m*q)/(2.*M_kz) + + else + + print *,'EASEv2_convert(): Polar projections not implemented yet' + stop + + endif + + row_ind = s0 - (y/map_scale_m) + col_ind = r0 + (x/map_scale_m) + + end subroutine easeV2_convert + + ! ******************************************************************* + + subroutine easeV2_inverse (grid, r, s, lat, lon) + + ! convert azimuthal equal area or equal area cylindrical + ! grid coordinates to geographic coordinates (spherical earth) + ! + ! *** NOTE order of calling arguments: "lon-lat-lat-lon" *** + ! + ! useage: call easeV1_inverse (grid, r, s, lat, lon) + ! + ! input : grid - projection name '[M][xx]' + ! where xx = approximate resolution [km] + ! ie xx = "01", "03", "09", "36" (SMAP) + ! r, s - column, row coordinates + ! + ! output: lat, lon = geo. coords. (decimal degrees) + ! + ! -------------------------------------------------------------------------- + + character*(*), intent(in) :: grid + real, intent(in) :: r, s + real, intent(out) :: lat, lon + + ! local variables + + integer :: cols, rows + real*8 :: phi, lam, map_scale_m, r0, s0, beta, x, y, qp + + real*8, parameter :: PI = easeV2_PI + + ! --------------------------------------------------------------------- + + call easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) + + x = (r - r0)*map_scale_m + y = -(s - s0)*map_scale_m + + qp = (1. - e2)* & + ( & + (1./(1.-e2)) & + - & + .5/map_eccentricity*log((1.-map_eccentricity)/(1.+map_eccentricity)) & + ) + + if (grid(1:1).eq.'M') then + + beta = asin(2.*y*M_kz/(map_equatorial_radius_m*qp)) + + lam = x/(map_equatorial_radius_m*M_kz) + + else + + print *,'EASEv2_inverse(): Polar projections not implemented yet' + stop + + endif + + phi = beta & + + ( ( e2/3. + 31./180.*e4 + 517./ 5040.*e6 )*sin(2.*beta) ) & + + ( ( 23./360.*e4 + 251./ 3780.*e6 )*sin(4.*beta) ) & + + ( ( 761./45360.*e6 )*sin(6.*beta) ) + + lat = phi*180./PI ! convert from radians to degree + lon = lam*180./PI + map_reference_longitude ! convert from radians to degree + + if (lon .lt. -180.0) lon = lon + 360.0 + if (lon .gt. 180.0) lon = lon - 360.0 + + end subroutine easeV2_inverse + + ! ******************************************************************* + + subroutine easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) + + implicit none + + character*(*), intent(in) :: grid + real*8, intent(out) :: map_scale_m, r0, s0 + integer, intent(out) :: cols, rows + + + if (grid(1:1).eq.'M') then + + if (grid .eq. 'M36') then ! SMAP 36 km grid + + map_scale_m = 36032.220840584 ! nominal cell size in meters + cols = 964 + rows = 406 + r0 = (cols-1)/2.0 + s0 = (rows-1)/2.0 + + + else if (grid .eq. 'M25') then ! 25 km grid + + map_scale_m = 25025.2600000 ! nominal cell size in meters (see doi:10.3390/ijgi3031154) + cols = 1388 + rows = 584 + r0 = (cols-1)/2.0 + s0 = (rows-1)/2.0 + + else if (grid .eq. 'M09') then ! SMAP 9 km grid + + map_scale_m = 9008.055210146 ! nominal cell size in meters + cols = 3856 + rows = 1624 + r0 = (cols-1)/2.0 + s0 = (rows-1)/2.0 + + else if (grid .eq. 'M03') then ! SMAP 3 km grid + + map_scale_m = 3002.6850700487 ! nominal cell size in meters + cols = 11568 + rows = 4872 + r0 = (cols-1)/2.0 + s0 = (rows-1)/2.0 + + else if (grid .eq. 'M01') then ! SMAP 1 km grid + + map_scale_m = 1000.89502334956 ! nominal cell size in meters + cols = 34704 + rows = 14616 + r0 = (cols-1)/2.0 + s0 = (rows-1)/2.0 + + else + + print *,'easeV2_get_params(): unknown resolution: ',grid + stop + + endif + + else if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then + + print *,'easeV2_get_params(): Polar projections not implemented yet' + stop + + else + + print *, 'easeV2_get_params(): unknown projection: ', grid + stop + + endif + + end subroutine easeV2_get_params + + ! ******************************************************************* + +end module EASE_conv + +! =============================== EOF ================================= + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/bcs_utils.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/bcs_utils.py new file mode 100755 index 000000000..59e8fce6a --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/bcs_utils.py @@ -0,0 +1,301 @@ +#!/usr/bin/env python3 +# +# source install/bin/g5_modules +# +# Newer GEOS code should load a module with GEOSpyD Python3 if not run: +# module load python/GEOSpyD/Min4.10.3_py3.9 +# + +import os +import socket +import subprocess +import shlex +import ruamel.yaml +import shutil +import questionary +import glob +from datetime import datetime + +def get_account(): + cmd = 'id -gn' + p = subprocess.Popen(shlex.split(cmd), stdout=subprocess.PIPE) + (accounts, err) = p.communicate() + p_status = p.wait() + accounts = accounts.decode().split() + return accounts[0] + +def get_config_from_answers(answers): + + grid_type = answers['grid_type'] + lbcsv = answers['bcs_version'] + resolution = answers['resolution'] + orslvs = answers['orslvs'] + + im = {} + imo = {} + jm = {} + jmo = {} + NX = 8640 + NY = 4320 + NT = 232000000 + + hostname = socket.gethostname() + input_dir = '' + if 'discover' in hostname: + input_dir = '/discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/' + else: + input_dir = '/nobackup/gmao_SIteam/ModelData/l_data/LandBCs_files_for_mkCatchParam/V001/' + + maskfile = '' + + if orslvs in['O1','T2','T3','T4','T1MOM6','T2MOM6','T4MOM6']: + maskfile = 'GEOS5_10arcsec_mask_freshwater-lakes.nc' + if lbcsv in ['F25', 'GM4', 'ICA']: + maskfile = 'global.cat_id.catch.DL' + + if orslvs in['O2','O3','CS']: + maskfile = 'GEOS5_10arcsec_mask.nc' + if lbcsv in ['F25', 'GM4', 'ICA']: + maskfile = 'global.cat_id.catch.GreatLakesCaspian_Updated.DL' + + if grid_type in ['EASEv1', 'EASEv2']: + maskfile = 'GEOS5_10arcsec_mask.nc' + + imo['O1'] = 360 + jmo['O1'] = 180 + + imo['O2'] = 1440 + jmo['O2'] = 720 + + imo['O3'] = 2880 + jmo['O3'] = 1440 + + imo['T2'] = 360 + jmo['T2'] = 200 + + imo['T3'] = 720 + jmo['T3'] = 410 + + imo['T4'] = 1440 + jmo['T4'] = 1080 + + imo['T1MOM6'] = 72 + jmo['T1MOM6'] = 36 + + imo['T2MOM6'] = 360 + jmo['T2MOM6'] = 210 + + imo['T4MOM6'] = 1440 + jmo['T4MOM6'] = 1080 + + im['b'] = 144 + jm['b'] = 91 + + im['c'] = 288 + jm['c'] = 181 + + im['d'] = 576 + jm['d'] = 361 + + im['e'] = 1152 + jm['e'] = 721 + + cubed = [12,24,48, 90, 180,360,720, 768,1000,1152, 1440,1536,2880,3072,5760] + for i in cubed: + key = 'c'+str(i) + im[key] = i + jm[key] = i*6 + + if grid_type == 'EASEv2' : + im['M01'] = 34704 + jm['M01'] = 14616 + + im['M03'] = 11568 + jm['M03'] = 4872 + + im['M09'] = 3856 + jm['M09'] = 1624 + + im['M25'] = 1388 + jm['M25'] = 584 + + im['M36'] = 964 + jm['M36'] = 406 + + if grid_type == 'EASEv1' : + im['M01'] = 34668 + jm['M01'] = 14688 + + im['M03'] = 11556 + jm['M03'] = 4896 + + im['M09'] = 3852 + jm['M09'] = 1632 + + im['M25'] = 1383 + jm['M25'] = 586 + + im['M36'] = 963 + jm['M36'] = 408 + + if resolution in ['c768','c1440'] : + NX = 17280 + NY = 8640 + if resolution == 'c2800': + NX = 21600 + NY = 10800 + if resolution in ['c1536', 'c3072','c5760'] : + NX = 43200 + NY = 21600 + + if 'GEOS5_10arcsec_mask' in maskfile: + NX = 43200 + NY = 21600 + + config = {} + config['grid_type'] = grid_type + config['lbcsv'] = lbcsv + config['resolution']= resolution + config['orslvs'] = orslvs + config ['im'] = im[resolution] + config ['jm'] = jm[resolution] + config ['imo'] = imo[orslvs] + config ['jmo'] = jmo[orslvs] + config ['NX'] = NX + config ['NY'] = NY + config ['NT'] = NT + config ['MASKFILE'] = maskfile + user = os.getlogin() + config ['expdir'] = '/discover/nobackup/'+user+'/BCS_PACKAGE/'+lbcsv+'/' + now = datetime.now() + config ['outdir'] =now.strftime("%Y%m%d%H%M%S") + config ['inputdir'] = input_dir + config ['NCPUS'] = 20 + + return config + +def ask_questions(): + + questions = [ + + { + "type": "select", + "name": "bcs_version", + "message": "Select land BCS version: \n \ + *BCs produced by this code will differ from BCs in archived directories!!! \n \ + These differences are caused by compiler changes, code improvements and bug \n \ + fixes that were implemented since the archived BCs in the above-mentioned \n \ + directories were originally created. The impact of these differences on \n \ + science is insignificant, and the parameter files produced by current \n \ + code are scientifically equivalent to the corresponding archived BCs. \n", + "choices": [ \ + "F25 : Fortuna-2_5 (archived*: n/a)", \ + "GM4 : Ganymed-4_0 (archived*: /discover/nobackup/ltakacs/bcs/Ganymed-4_0/)", \ + "ICA : Icarus (archived*: /discover/nobackup/ltakacs/bcs/Icarus/)", \ + "NL3 : Icarus-NLv3 (archived*: /discover/nobackup/ltakacs/bcs/Icarus-NLv3/)", \ + "NL4 : NLv4 [SMAPL4] (archived*: /discover/nobackup/projects/gmao/smap/bcs_NLv4/NLv4/)", \ + "NL5 : NLv5 [SMAPL4] (archived*: /discover/nobackup/projects/gmao/smap/SMAP_L4/L4_SM/bcs/CLSM_params/Icarus-NLv5_EASE/)", \ + "DEV : Development version"], + "default": "NL3 : Icarus-NLv3 (archived*: /discover/nobackup/ltakacs/bcs/Icarus-NLv3/)", + }, + + { + "type": "select", + "name": "grid_type", + "message": "Select grid type: \n ", + "choices": ["Lat-Lon", "Cubed-Sphere", "EASEv2", "EASEv1"], + "default": "Cubed-Sphere", + }, + + { + "type": "select", + "name": "resolution", + "message": "Select lat-lon resolution: \n ", + "choices": ["b -- 2 deg", "c -- 1 deg", "d -- 1/2 deg","e -- 1/4 deg"], + "default": "d -- 1/2 deg", + "when": lambda x: x['grid_type'] == "Lat-Lon", + }, + + { + "type": "select", + "name": "resolution", + "message": "Select cubed-sphere resolution: \n ", + "choices": [ \ + "c12 -- 8 deg", \ + "c24 -- 4 deg", \ + "c48 -- 2 deg", \ + "c90 -- 1 deg", \ + "c180 -- 1/2 deg ( 56 km)", \ + "c360 -- 1/4 deg ( 28 km)", \ + "c720 -- 1/8 deg ( 14 km)", \ + "c768 -- 1/10 deg ( 12 km)", \ + "c1000 -- 1/10 deg ( 10 km)", \ + "c1152 -- 1/10 deg ( 8 km)", \ + "c1440 -- 1/16 deg ( 7 km)", \ + "c1536 -- 1/16 deg ( 7 km)", \ + "c2880 -- 1/32 deg ( 3 km)", \ + "c3072 -- 1/32 deg ( 3 km)", \ + "c5760 -- 1/64 deg ( 1.5 km)"], + "default": "c360 -- 1/4 deg ( 28 km)", + "when": lambda x: x['grid_type'] == "Cubed-Sphere", + }, + + { + "type": "select", + "name": "resolution", + "message": "Select EASE grid resolution: \n ", + "choices": [ \ + "M01 -- 1km", \ + "M03 -- 3km", \ + "M09 -- 9km", \ + "M25 -- 25km", \ + "M36 -- 36km"], + "default": "M09 -- 9km", + "when": lambda x: x['grid_type'] == "EASEv2" or x['grid_type'] == "EASEv1", + }, + + { + "type": "select", + "name": "orslvs", + "message": "Select ocean resolution: \n ", + "choices": [ \ + "O1 -- Low-Resolution Reynolds 1 deg (Lon/Lat Data-Ocean: 360x180 )", \ + "O2 -- Med-Resolution Reynolds 1/4 deg (Lon/Lat Data-Ocean: 1440x720 )", \ + "O3 -- High-Resolution OSTIA 1/8 deg (Lon/Lat Data-Ocean: 2880x1440)", \ + "T2 -- Med-Resolution Tripolar 1 deg (MOM-Tripolar-Ocean: 360x200 )", \ + "T3 -- High-Resolution Tripolar 1/2 deg (MOM-Tripolar-Ocean: 720x410 )", \ + "T4 -- High-Resolution Tripolar 1/4 deg (MOM-Tripolar-Ocean: 1440x1080)", \ + "T1MOM6 -- Low-Resolution Tripolar 5 deg (MOM6-Tripolar-Ocean: 72x36 )", \ + "T2MOM6 -- Med-Resolution Tripolar 1 deg (MOM6-Tripolar-Ocean: 360x210 )", \ + "T4MOM6 -- High-Resolution Tripolar 1/4 deg (MOM6-Tripolar-Ocean: 1440x1080)", \ + "CS -- Cubed-Sphere Ocean (Cubed-Sphere Data-Ocean )"], + "when": lambda x: x['grid_type'] == "Lat-Lon" or x['grid_type'] == "Cubed-Sphere", + }, + + + ] + answers_ = questionary.prompt(questions) + answers = {} + for key, value in answers_.items(): + answers[key] = value.split()[0] + if 'EASE' in answers['grid_type'] : answers['orslvs'] = 'O1' + + return answers + +def print_config( config, indent = 0 ): + for k, v in config.items(): + if isinstance(v, dict): + print(" " * indent, f"{k}:") + print_config(v, indent+1) + else: + print(" " * indent, f"{k}: {v}") + +if __name__ == "__main__": + + answers = ask_questions() + for key, value in answers.items(): + print(key, value) + + config = get_config_from_answers(answers) + print('\n print config file:\n') + print_config(config) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/clsm_plots.pro b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/clsm_plots.pro index a15739012..cfb19b470 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/clsm_plots.pro +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/clsm_plots.pro @@ -257,7 +257,7 @@ dy = NR/NR_plot catrow = lonarr(nc) cat = lonarr(nc,dy) -rst_file=path + '/rst/' + gfile+'.rst' +rst_file=path + '/rst/' + gfile+'*.rst' openr,1,rst_file,/F77_UNFORMATTED for j = 0l, NR_plot -1 do begin @@ -570,7 +570,7 @@ dy = NR/NR_movie cat = lonarr(nc,dy) catrow = lonarr(nc) -rst_file=path + '/rst/' + gfile+'.rst' +rst_file=path + '/rst/' + gfile+'*.rst' openr,1,rst_file,/F77_UNFORMATTED for j = 0l, NR_movie -1 do begin @@ -1340,7 +1340,7 @@ pfcr=0l cat =lonarr(nc) catp=lonarr(nc) -rst_file=path + '/rst/' + gfile+'.rst' +rst_file=path + '/rst/' + gfile+'*.rst' idum=0l openr,1,rst_file,/F77_UNFORMATTED @@ -1708,7 +1708,7 @@ dy = JM/NR catrow = lonarr (nc) tile_id = lonarr (NC, nr) -rst_file= path + '/rst/' + gfile+'.rst' +rst_file= path + '/rst/' + gfile+'*.rst' openr,1,rst_file,/F77_UNFORMATTED diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/easeV1_conv.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/easeV1_conv.F90 deleted file mode 100644 index 859c467a1..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/easeV1_conv.F90 +++ /dev/null @@ -1,305 +0,0 @@ - -module easeV1_conv - - ! ========================================================================== - ! - ! easeV1_conv.F90 - FORTRAN routines for conversion of azimuthal - ! equal area and equal area cylindrical grid coordinates - ! - ! 30-Jan-1992 H.Maybee - ! 20-Mar-1992 Ken Knowles 303-492-0644 knowles@kryos.colorado.edu - ! 16-Dec-1993 MJ Brodzik 303-492-8263 brodzik@jokull.colorado.edu - ! Copied from nsmconv.f, changed resolutions from - ! 40-20-10 km to 25-12.5 km - ! 21-Dec-1993 MJ Brodzik 303-492-8263 brodzik@jokull.colorado.edu - ! Fixed sign of Southern latitudes in ease_inverse. - ! 12-Sep-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu - ! Changed grid cell size. Changed "c","f" to "l","h" - ! 25-Oct-1994 David Hoogstrate 303-492-4116 hoogstra@jokull.colorado.edu - ! Changed row size from 587 to 586 for Mercator projection - ! 11-May-2011 reichle: Changed "smap" to "easeV1". - ! Added SSM/I and AMSR-E "M25" grid. - ! So far ONLY for cylindrical grids. - ! Converted from *.f to *.F90 module - ! - ! $Log$ - ! Revision 1.1.4.2 2012/04/18 18:04:57 smahanam - ! updated on 4-18-2012 to process global data on native grids and time steps - ! - ! Revision 1.1.2.1 2012-01-27 20:22:10 smahanam - ! Added Richard's equation solver, replaced SMAP f77 files with F90 module - ! - ! Revision 1.1 2011-05-11 21:58:46 rreichle - ! - ! Adding utilities to map between EASE grids and lat/lon coordinates. - ! - ! Revision 1.3 1994/11/01 23:40:43 brodzik - ! Replaced all references to 'ease' with 'smap' - ! Replaced all references to 'smap' with 'easeV1' -- reichle - ! - ! ========================================================================== - - implicit none - - private - - public :: easeV1_convert - public :: easeV1_inverse - - - ! ***NEVER*** change these constants to GEOS-5 MAPL constants!!!! - - ! radius of the earth (km), authalic sphere based on International datum - - real*8, parameter :: RE_km = 6371.228 - - ! scale factor for standard paralles at +/-30.00 degrees - - real*8, parameter :: COS_PHI1 = .866025403 - - real*8, parameter :: PI = 3.141592653589793 - - -contains - - ! ******************************************************************* - - subroutine easeV1_convert (grid, lat, lon, r, s) - - ! convert geographic coordinates (spherical earth) to - ! azimuthal equal area or equal area cylindrical grid coordinates - ! - ! status = easeV1_convert (grid, lat, lon, r, s) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! or xx = "12", "25" (SSM/I, AMSR-E) - ! lat, lon = geo. coords. (decimal degrees) - ! - ! output: r, s - column, row coordinates - ! - ! result: status = 0 indicates normal successful completion - ! -1 indicates error status (point not on grid) - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: lat, lon - real, intent(out) :: r, s - - ! local variables - - integer :: cols, rows, scale - real*8 :: Rg, phi, lam, rho, CELL_km, r0, s0 - - ! --------------------------------------------------------------------- - - call easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) - - phi = lat*PI/180. ! convert from degree to radians - lam = lon*PI/180. ! convert from degree to radians - - if (grid(1:1).eq.'N') then - rho = 2 * Rg * sin(PI/4. - phi/2.) - r = r0 + rho * sin(lam) - s = s0 + rho * cos(lam) - - else if (grid(1:1).eq.'S') then - rho = 2 * Rg * cos(PI/4. - phi/2.) - r = r0 + rho * sin(lam) - s = s0 - rho * cos(lam) - - else if (grid(1:1).eq.'M') then - r = r0 + Rg * lam * COS_PHI1 - s = s0 - Rg * sin(phi) / COS_PHI1 - - endif - - end subroutine easeV1_convert - - ! ******************************************************************* - - subroutine easeV1_inverse (grid, r, s, lat, lon) - - ! convert azimuthal equal area or equal area cylindrical - ! grid coordinates to geographic coordinates (spherical earth) - ! - ! status = easeV1_inverse (grid, r, s, lat, lon) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! or xx = "12", "25" (SSM/I, AMSR-E) - ! r, s - column, row coordinates - ! - ! output: lat, lon = geo. coords. (decimal degrees) - ! - ! result: status = 0 indicates normal successful completion - ! -1 indicates error status (point not on grid) - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: r, s - real, intent(out) :: lat, lon - - ! local variables - - integer :: cols, rows, scale - real*8 :: Rg, phi, lam, rho, CELL_km, r0, s0 - real*8 :: gamma, beta, epsilon, x, y, c - real*8 :: sinphi1, cosphi1 - - ! --------------------------------------------------------------------- - - call easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) - - x = r - r0 - y = -(s - s0) - - if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then - rho = sqrt(x*x + y*y) - if (rho.eq.0.0) then - if (grid(1:1).eq.'N') lat = 90.0 - if (grid(1:1).eq.'S') lat = -90.0 - lon = 0.0 - else - if (grid(1:1).eq.'N') then - sinphi1 = sin(PI/2.) - cosphi1 = cos(PI/2.) - if (y.eq.0.) then - if (r.le.r0) lam = -PI/2. - if (r.gt.r0) lam = PI/2. - else - lam = atan2(x,-y) - endif - else if (grid(1:1).eq.'S') then - sinphi1 = sin(-PI/2.) - cosphi1 = cos(-PI/2.) - if (y.eq.0.) then - if (r.le.r0) lam = -PI/2. - if (r.gt.r0) lam = PI/2. - else - lam = atan2(x,y) - endif - endif - gamma = rho/(2 * Rg) - if (abs(gamma) .gt. 1.) return - c = 2 * asin(gamma) - beta = cos(c) * sinphi1 + y * sin(c) * (cosphi1/rho) - if (abs(beta).gt.1.) return - phi = asin(beta) - lat = phi*180./PI ! convert from radians to degree - lon = lam*180./PI ! convert from radians to degree - endif - - else if (grid(1:1).eq.'M') then - - ! allow .5 cell tolerance in arcsin function - ! so that grid coordinates which are less than .5 cells - ! above 90.00N or below 90.00S are given a lat of 90.00 - - epsilon = 1 + 0.5/Rg - beta = y*COS_PHI1/Rg - if (abs(beta).gt.epsilon) return - if (beta.le.-1.) then - phi = -PI/2. - else if (beta.ge.1.) then - phi = PI/2. - else - phi = asin(beta) - endif - lam = x/COS_PHI1/Rg - lat = phi*180./PI ! convert from radians to degree - lon = lam*180./PI ! convert from radians to degree - endif - - end subroutine easeV1_inverse - - ! ******************************************************************* - - subroutine easeV1_get_params( grid, CELL_km, cols, rows, r0, s0, Rg ) - - implicit none - - character*(*), intent(in) :: grid - real*8, intent(out) :: CELL_km, r0, s0, Rg - integer, intent(out) :: cols, rows - - ! -------------------------------------------------------- - ! - ! r0,s0 are defined such that cells at all scales have - ! coincident center points - ! - !c r0 = (cols-1)/2. * scale - !c s0 = (rows-1)/2. * scale - ! - ! -------------------------------------------------------- - - if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then - - print *,'Polar projections not implemented yet' - stop - - else if (grid(1:1).eq.'M') then - - if (grid .eq. 'M36') then ! SMAP 36 km grid - CELL_km = 36.00040279063 ! nominal cell size in kilometers - cols = 963 - rows = 408 - r0 = 481.0 - s0 = 203.5 - - else if (grid .eq. 'M25') then ! SSM/I, AMSR-E 25 km grid - CELL_km = 25.067525 ! nominal cell size in kilometers - cols = 1383 - rows = 586 - r0 = 691.0 - s0 = 292.5 - - else if (grid .eq. 'M09') then ! SMAP 9 km grid - CELL_km = 9.00010069766 ! nominal cell size in kilometers - cols = 3852 - rows = 1632 - r0 = 1925.5 - s0 = 815.5 - - else if (grid .eq. 'M03') then ! SMAP 3 km grid - CELL_km = 3.00003356589 ! nominal cell size in kilometers - cols = 11556 - rows = 4896 - r0 = 5777.5 - s0 = 2447.5 - - else if (grid .eq. 'M01') then ! SMAP 1 km grid - CELL_km = 1.00001118863 ! nominal cell size in kilometers - cols = 34668 - rows = 14688 - r0 = 17333.5 - s0 = 7343.5 - - else - - print *,'easeV1_convert: unknown resolution: ',grid - stop - - endif - - else - - print *, 'easeV1_convert: unknown projection: ', grid - stop - - endif - - Rg = RE_km/CELL_km - - end subroutine easeV1_get_params - - ! ******************************************************************* - -end module easeV1_conv - -! =============================== EOF ================================= - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/easeV2_conv.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/easeV2_conv.F90 deleted file mode 100644 index 5b7370e6b..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/easeV2_conv.F90 +++ /dev/null @@ -1,317 +0,0 @@ - -module easeV2_conv - - ! ========================================================================== - ! - ! easeV2_conv.F90 - FORTRAN routines for converting grid coordinates - ! (latitude/longitude <--> row/column indices) - ! of the Equal Area Scalable Earth, version 2 (EASEv2) grid - ! - ! ***** ONLY cylindrical ('M') projection implemented ***** - ! - ! Ported from Steven Chan's matlab code (smapease2inverse.m, - ! smapease2forward.m), which has been ported from NSIDC's IDL code - ! (wgs84_convert.pro, wgs84_inverse.pro) available from - ! ftp://sidads.colorado.edu/pub/tools/easegrid/geolocation_tools/ - ! - ! Official references: - ! doi:10.3390/ijgi1010032 - ! doi:10.3390/ijgi3031154 -- correction of M25 "map_scale_m" parameters! - ! - ! 04-Apr-2013 - reichle - ! 11-Sep-2018 - reichle, mgirotto -- added 'M25' grid parameters - ! - ! ========================================================================== - - implicit none - - private - - public :: easeV2_convert - public :: easeV2_inverse - public :: easeV2_extent - - ! ***NEVER*** change these constants to GEOS-5 MAPL constants!!!! - - ! radius of the earth (m) and map eccentricity - - real*8, parameter :: map_equatorial_radius_m = 6378137.0 - - real*8, parameter :: map_eccentricity = 0.081819190843 - - real*8, parameter :: PI = 3.14159265358979323846 - - real*8, parameter :: e2 = map_eccentricity * map_eccentricity - real*8, parameter :: e4 = e2 * e2 - real*8, parameter :: e6 = e2 * e4 - - real*8, parameter :: epsilon = 1.e-6 - - real*8, parameter :: map_reference_longitude = 0.0 ! 'M', 'N', 'S' - - ! constants for 'N' and 'S' (azimuthal) projections - - real*8, parameter :: N_map_reference_latitude = 90.0 - real*8, parameter :: S_map_reference_latitude = -90.0 - - ! constants for 'M' (cylindrical) projection - - real*8, parameter :: M_map_reference_latitude = 0.0 - real*8, parameter :: M_map_second_reference_latitude = 30.0 - - real*8, parameter :: M_sin_phi1 = sin(M_map_second_reference_latitude*PI/180.) - real*8, parameter :: M_cos_phi1 = cos(M_map_second_reference_latitude*PI/180.) - - real*8, parameter :: M_kz = M_cos_phi1/sqrt(1.0-e2*M_sin_phi1*M_sin_phi1) - - -contains - - ! ******************************************************************* - - subroutine easeV2_convert (grid, lat, lon, col_ind, row_ind) - - ! convert geographic coordinates (spherical earth) to - ! azimuthal equal area or equal area cylindrical grid coordinates - ! - ! *** NOTE order of calling arguments: "lat-lon-lon-lat" *** - ! - ! useage: call easeV2_convert (grid, lat, lon, r, s) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! lat, lon = geo. coords. (decimal degrees) - ! - ! output: col_ind, row_ind - column, row coordinates - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: lat, lon - real, intent(out) :: col_ind, row_ind - - ! local variables - - integer :: cols, rows - real*8 :: dlon, phi, lam, map_scale_m, r0, s0, ms, x, y, sin_phi, q - - ! --------------------------------------------------------------------- - - call easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) - - dlon = lon - - if (abs(map_reference_longitude)>epsilon) then - - dlon = lon - map_reference_longitude - - end if - - if (dlon .lt. -180.0) dlon = dlon + 360.0 - if (dlon .gt. 180.0) dlon = dlon - 360.0 - - phi = lat*PI/180. ! convert from degree to radians - lam = dlon*PI/180. ! convert from degree to radians - - sin_phi = sin(phi) - - ms = map_eccentricity*sin_phi - - q = (1. - e2)* & - ( & - (sin_phi /(1. - e2*sin_phi*sin_phi)) & - - & - .5/map_eccentricity*log((1.-ms)/(1.+ms)) & - ) - - ! note: "qp" only needed for 'N' and 'S' projections - - if (grid(1:1).eq.'M') then - - x = map_equatorial_radius_m*M_kz*lam - - y = (map_equatorial_radius_m*q)/(2.*M_kz) - - else - - print *,'Polar projections not implemented yet' - stop - - endif - - row_ind = s0 - (y/map_scale_m) - col_ind = r0 + (x/map_scale_m) - - end subroutine easeV2_convert - - ! ******************************************************************* - - subroutine easeV2_inverse (grid, r, s, lat, lon) - - ! convert azimuthal equal area or equal area cylindrical - ! grid coordinates to geographic coordinates (spherical earth) - ! - ! *** NOTE order of calling arguments: "lon-lat-lat-lon" *** - ! - ! useage: call easeV1_inverse (grid, r, s, lat, lon) - ! - ! input : grid - projection name '[M][xx]' - ! where xx = approximate resolution [km] - ! ie xx = "01", "03", "09", "36" (SMAP) - ! r, s - column, row coordinates - ! - ! output: lat, lon = geo. coords. (decimal degrees) - ! - ! -------------------------------------------------------------------------- - - character*(*), intent(in) :: grid - real, intent(in) :: r, s - real, intent(out) :: lat, lon - - ! local variables - - integer :: cols, rows - real*8 :: phi, lam, map_scale_m, r0, s0, beta, x, y, qp - - ! --------------------------------------------------------------------- - - call easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) - - x = (r - r0)*map_scale_m - y = -(s - s0)*map_scale_m - - qp = (1. - e2)* & - ( & - (1./(1.-e2)) & - - & - .5/map_eccentricity*log((1.-map_eccentricity)/(1.+map_eccentricity)) & - ) - - if (grid(1:1).eq.'M') then - - beta = asin(2.*y*M_kz/(map_equatorial_radius_m*qp)) - - lam = x/(map_equatorial_radius_m*M_kz) - - else - - print *,'Polar projections not implemented yet' - stop - - endif - - phi = beta & - + ( ( e2/3. + 31./180.*e4 + 517./ 5040.*e6 )*sin(2.*beta) ) & - + ( ( 23./360.*e4 + 251./ 3780.*e6 )*sin(4.*beta) ) & - + ( ( 761./45360.*e6 )*sin(6.*beta) ) - - lat = phi*180./PI ! convert from radians to degree - lon = lam*180./PI + map_reference_longitude ! convert from radians to degree - - if (lon .lt. -180.0) lon = lon + 360.0 - if (lon .gt. 180.0) lon = lon - 360.0 - - end subroutine easeV2_inverse - - ! ******************************************************************* - - subroutine easeV2_get_params( grid, map_scale_m, cols, rows, r0, s0 ) - - implicit none - - character*(*), intent(in) :: grid - real*8, intent(out) :: map_scale_m, r0, s0 - integer, intent(out) :: cols, rows - - - if (grid(1:1).eq.'M') then - - if (grid .eq. 'M36') then ! SMAP 36 km grid - - map_scale_m = 36032.220840584 ! nominal cell size in meters - cols = 964 - rows = 406 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else if (grid .eq. 'M25') then ! 25 km grid - - map_scale_m = 25025.2600000 ! nominal cell size in meters (see doi:10.3390/ijgi3031154) - cols = 1388 - rows = 584 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else if (grid .eq. 'M09') then ! SMAP 9 km grid - - map_scale_m = 9008.055210146 ! nominal cell size in meters - cols = 3856 - rows = 1624 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else if (grid .eq. 'M03') then ! SMAP 3 km grid - - map_scale_m = 3002.6850700487 ! nominal cell size in meters - cols = 11568 - rows = 4872 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else if (grid .eq. 'M01') then ! SMAP 1 km grid - - map_scale_m = 1000.89502334956 ! nominal cell size in meters - cols = 34704 - rows = 14616 - r0 = (cols-1)/2.0 - s0 = (rows-1)/2.0 - - else - - print *,'easeV2_convert: unknown resolution: ',grid - stop - - endif - - else if ((grid(1:1).eq.'N').or.(grid(1:1).eq.'S')) then - - print *,'Polar projections not implemented yet' - stop - - else - - print *, 'easeV2_convert: unknown projection: ', grid - stop - - endif - - end subroutine easeV2_get_params - - ! ******************************************************************* - - subroutine easeV2_extent( grid, N_cols, N_rows ) - - ! simple wrapper to get N_cols (N_lon) and N_rows (N_lat) - - implicit none - - character*(*), intent(in) :: grid - integer, intent(out) :: N_cols, N_rows - - ! local variables - - real*8 :: map_scale_m, r0, s0 - - ! ------------------------------------------------ - - call easeV2_get_params( grid, map_scale_m, N_cols, N_rows, r0, s0 ) - - end subroutine easeV2_extent - - ! ******************************************************************* - -end module easeV2_conv - -! =============================== EOF ================================= - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs index af40fb04b..8d6d0b50f 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs @@ -9,7 +9,7 @@ if ( "$1" == "-h" | "$1" == "-help" | "$1" == "--help" ) set HELPMODE = YES if ( $HELPMODE == YES ) then - echo "Usage: `basename $0` [option] " + echo "Usage: `basename $0` [-h] " echo " " echo "Boundary Conditions (BCs) Package: " echo " Creates surface tile and other model parameter input files " @@ -24,8 +24,9 @@ if ( $HELPMODE == YES ) then echo " Answer the following interactive questions: " echo " a) Select Land BCs version. " echo " b) Select atmospheric resolution(s). " + echo " (If applicable, select EASE grid version.) " echo " c) Select ocean resolution(s). " - echo " (Not relevant for land-only EASE-grid BCs.) " + echo " (Not relevant for land-only EASE grid BCs.) " echo " d) Enter BCs output directory. " echo " e) Enter sponsor code for computing account. " echo " " @@ -96,14 +97,11 @@ if ( $HELPMODE != YES ) then echo "----------------------------------------------------------------------------------------------" echo " Boundary Conditions (BCs) Package:" echo "----------------------------------------------------------------------------------------------" - echo " " - echo "Enter 3-character alphanumeric code for land BCs version:" - echo "(Select only one.)" - echo " " -else - echo "Options for land BCs version:" endif - + +echo " " +echo "Options for land BCs version:" +echo " " echo " ${C2}F25 = Fortuna-2_5 (archived${CR}${C1}*${CR}${C2}: n/a)${CR}" echo " ${C2}GM4 = Ganymed-4_0 (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/ltakacs/bcs/Ganymed-4_0/)${CR}" echo " ${C2}ICA = Icarus (archived${CR}${C1}*${CR}${C2}: /discover/nobackup/ltakacs/bcs/Icarus/)${CR}" @@ -125,16 +123,13 @@ if ( $HELPMODE != YES ) then echo " science is insignificant, and the parameter files produced by current" echo " code are scientifically equivalent to the corresponding archived BCs." echo " " - echo " OR press ENTER to select $lbcsv (current default).${CR}" + echo " Enter one 3-character alphanumeric code for land BCs version from the above list" + echo " OR press ENTER to select $lbcsv (current default):${CR}" echo " " -endif - -LBCSV: -if ( $HELPMODE != YES ) then + LBCSV: set dummy = `echo $<` - if( $dummy == 'F25' | \ $dummy == 'GM4' | \ $dummy == 'ICA' | \ @@ -150,7 +145,7 @@ if ( $HELPMODE != YES ) then echo $lbcsv else echo " " - echo " ${C1} Invalid choice, try again:${CR}" + echo " ${C1} Invalid choice. Note that entry is case-sensitive. Try again:${CR}" goto LBCSV endif @@ -160,19 +155,14 @@ endif # # Atmospheric Horizontal Resolution (HRCODE) -if ( $HELPMODE != YES ) then - echo " " - echo "Enter 1-5 character alphanumeric code(s) for atmospheric horizontal resolution:" - echo "(Separate multiple entries by spaces.)" -else - echo "Options for atmospheric horizontal resolution:" -endif +echo "Options for atmospheric horizontal resolution:" +echo " " echo " Lat/Lon Cubed-Sphere EASE (land-only)" -echo " ${C2}b -- 2 deg c12 -- 8 deg m1 -- 1km EASEv2 Grid ${CR}" -echo " ${C2}c -- 1 deg c24 -- 4 deg m3 -- 3km EASEv2 Grid ${CR}" -echo " ${C2}d -- 1/2 deg c48 -- 2 deg m9 -- 9km EASEv2 Grid ${CR}" -echo " ${C2}e -- 1/4 deg c90 -- 1 deg m36 -- 36km EASEv2 Grid ${CR}" -echo " ${C2} c180 -- 1/2 deg ${C1}( 56 km) ${CR} ${C2}m25 -- 25km EASEv1 Grid ${CR}" +echo " ${C2}b -- 2 deg c12 -- 8 deg M01 -- 1-km EASE Grid ${CR}" +echo " ${C2}c -- 1 deg c24 -- 4 deg M03 -- 3-km EASE Grid ${CR}" +echo " ${C2}d -- 1/2 deg c48 -- 2 deg M09 -- 9-km EASE Grid ${CR}" +echo " ${C2}e -- 1/4 deg c90 -- 1 deg M25 -- 25-km EASE Grid ${CR}" +echo " ${C2} c180 -- 1/2 deg ${C1}( 56 km) ${CR} ${C2}M36 -- 36-km EASE Grid ${CR}" echo " ${C2} c360 -- 1/4 deg ${C1}( 28 km) ${CR}" echo " ${C2} c720 -- 1/8 deg ${C1}( 14 km) ${CR}" echo " ${C2} c768 -- 1/10 deg ${C1}( 12 km) ${CR}" @@ -184,15 +174,18 @@ echo " ${C2} c2880 -- 1/32 deg ${C1}( 3 km) ${CR}" echo " ${C2} c3072 -- 1/32 deg ${C1}( 3 km) ${CR}" echo " ${C2} c5760 -- 1/64 deg ${C1}( 1.5 km) ${CR}" #echo " ${C2}o -- other${CR} Lat/Lon or Cube" -echo " " +echo " " -HRCODE: if ( $HELPMODE != YES ) then + echo " " + echo "Enter alphanumeric code(s) for atmospheric horizontal resolution from the above list" + echo "(separate multiple entries by spaces):" + + HRCODE: set dummy = `echo $<` if (${%dummy} == 0) set dummy = "BLANK" - set dummy = `echo $dummy | tr "[:upper:]" "[:lower:]"` set HRCODES = $dummy[1] foreach HRCODE ($dummy) @@ -216,12 +209,12 @@ if ( $HELPMODE != YES ) then $HRCODE != 'c2880' & \ $HRCODE != 'c3072' & \ $HRCODE != 'c5760' & \ - $HRCODE != 'm1' & \ - $HRCODE != 'm3' & \ - $HRCODE != 'm9' & \ - $HRCODE != 'm36' & \ - $HRCODE != 'm25') then - echo " ${C1} Invalid choice, try again:${CR}" + $HRCODE != 'M01' & \ + $HRCODE != 'M03' & \ + $HRCODE != 'M09' & \ + $HRCODE != 'M25' & \ + $HRCODE != 'M36') then + echo " ${C1} Invalid choice. Note that entries are case-sensitive. Try again:${CR}" goto HRCODE endif if( $HRCODE != "$HRCODES" ) set HRCODES = `echo ${HRCODES} ${HRCODE}` @@ -233,27 +226,64 @@ if ( $HELPMODE != YES ) then set Resolution = `echo $<` endif - # figure out if one or more of the selected atmospheric resolutions is an EASE grid - - set isEASE = 0 - - foreach HRCODE ($HRCODES) - if( $HRCODE == 'm1' | \ - $HRCODE == 'm3' | \ - $HRCODE == 'm9' | \ - $HRCODE == 'm36' | \ - $HRCODE == 'm25') then - set isEASE = 1 - endif - end - else set HRCODE = null set HRCODES = null -endif +endif # if ( $HELPMODE != YES ) + +# -------------------------------- +# +# EASE Grid Version (EASEVERSION) + +# set default EASE version +set EASEVERSION = EASEv2 + +echo " Options for EASE grid version (land-only):" +echo " ${C2}EASEv1${CR}" +echo " ${C2}EASEv2 (default)${CR}" +echo " " +set isEASE = 0 + +if ( $HELPMODE != YES ) then + + # figure out if one or more of the selected atmospheric resolutions is an EASE grid + + foreach HRCODE ($HRCODES) + if( $HRCODE == 'M01' | \ + $HRCODE == 'M03' | \ + $HRCODE == 'M09' | \ + $HRCODE == 'M25' | \ + $HRCODE == 'M36') then + set isEASE = 1 + endif + end + + if ( $isEASE > 0 ) then + + echo "Enter one EASE grid version from the above list" + echo " OR press ENTER to select $EASEVERSION (current default):${CR}" + echo " " + + EASE_VERSION: + set dummy = `echo $<` + + if( $dummy == 'EASEv1' | \ + $dummy == 'EASEv2') then + set EASEVERSION = $dummy + else if ( $dummy == '' ) then + echo $EASEVERSION + else + echo " " + echo " ${C1} Invalid choice. Note that entry is case-sensitive. Try again:${CR}" + echo " " + goto EASE_VERSION + endif + endif + +endif ####################################################################### # @@ -262,21 +292,16 @@ endif if ($isEASE > 0) then set orslvs = O1 - echo " ${C1}---------------------------------------------------------------------- ${CR}" - echo " ${C1} Selected atmospheric resolution(s) include(s) EASE grid(s). ${CR}" - echo " ${C1} Setting ocean resolution to $orslvs \!\!\! ${CR}" - echo " ${C1} For a choice of ocean resolutions, exclude EASE grid selection(s) ${CR}" - echo " ${C1} from selected atmospheric resolution(s) and process separately. ${CR}" - echo " ${C1}-----------------------------------------------------------------------${CR}" + echo " ${C1}----------------------------------------------------------------------- ${CR}" + echo " ${C1} Selected atmospheric resolution(s) include(s) at least one EASE grid.${CR}" + echo " ${C1} Setting ${CR}ocean resolution ${C1}to $orslvs \!\!\! ${CR}" + echo " ${C1} For a choice of ocean resolutions, exclude EASE grid selection(s) ${CR}" + echo " ${C1} from selected atmospheric resolution(s) and process separately. ${CR}" + echo " ${C1}------------------------------------------------------------------------${CR}" echo " " else - if ( $HELPMODE != YES ) then - echo " " - echo "Enter 2-6 character alphanumeric code(s) for ocean horizontal resolution:" - echo "(Separate multiple entries by spaces.)" - else - echo "Options for ocean horizontal resolution:" - endif + echo "Options for ocean horizontal resolution:" + echo " " echo " ${C2}O1 -- Low-Resolution Reynolds 1 deg${CR} (Lon/Lat Data-Ocean: 360x180 )" echo " ${C2}O2 -- Med-Resolution Reynolds 1/4 deg${CR} (Lon/Lat Data-Ocean: 1440x720 )" echo " ${C2}O3 -- High-Resolution OSTIA 1/8 deg${CR} (Lon/Lat Data-Ocean: 2880x1440)" @@ -291,13 +316,16 @@ else if ( $HELPMODE == YES ) exit + echo " " + echo "Enter 2-6 character alphanumeric code(s) for ocean horizontal resolution from the above list." + echo "(Separate multiple entries by spaces.)" + ORSLV: set dummy = `echo $<` if (${%dummy} == 0) set dummy = "BLANK" - set dummy = `echo $dummy | tr "[:lower:]" "[:upper:]"` set orslvs = $dummy[1] foreach orslv ($dummy) @@ -311,7 +339,7 @@ ORSLV: $orslv != 'T2MOM6' & \ $orslv != 'T4MOM6' & \ $orslv != 'CS') then - echo " ${C1} Invalid choice, try again:${CR}" + echo " ${C1} Invalid choice. Note that entries are case-sensitive. Try again:${CR}" goto ORSLV endif if( $orslv != "$orslvs" ) set orslvs = `echo ${orslvs} ${orslv}` @@ -364,11 +392,14 @@ endif ####################################################################### -echo "" -echo "${C1} Land BCs version:${CR} ${C2}$lbcsv${CR}" -echo "${C1} Atmospheric resolution:${CR} ${C2}$HRCODES${CR}" -echo "${C1} Ocean resolution:${CR} ${C2}$orslvs${CR}" -echo "${C1} Output directory:${CR} ${C2}$EXPDIR${CR}" +echo "" +echo "${C1} Land BCs version:${CR} ${C2}$lbcsv${CR}" +echo "${C1} Atmospheric resolution:${CR} ${C2}$HRCODES${CR}" +if ( $isEASE > 0) then + echo "${C1} EASE grid version:${CR} ${C2}$EASEVERSION${CR}" +endif +echo "${C1} Ocean resolution:${CR} ${C2}$orslvs${CR}" +echo "${C1} Output directory:${CR} ${C2}$EXPDIR${CR}" ####################################################################### ####################################################################### @@ -426,7 +457,7 @@ else endif -if($HRCODE == m1 | $HRCODE == m3 | $HRCODE == m9 | $HRCODE == m36 | $HRCODE == m25) then +if($HRCODE == M01 | $HRCODE == M03 | $HRCODE == M09 | $HRCODE == M25 | $HRCODE == M36) then set GLOBAL_CATCH_DATA = ${input_dir}/GEOS5_10arcsec_mask.nc endif @@ -628,36 +659,68 @@ if( $HRCODE == c5760 ) then @ NX = 43200 @ NY = 21600 endif -if( $HRCODE == m1 ) then - set im = 34704 - set jm = 14616 - set grid = ease - set MGRID = M01 -endif -if( $HRCODE == m3 ) then - set im = 11568 - set jm = 4872 - set grid = ease - set MGRID = M03 -endif -if( $HRCODE == m9 ) then - set im = 3856 - set jm = 1624 - set grid = ease - set MGRID = M09 -endif -if( $HRCODE == m36 ) then - set im = 964 - set jm = 406 - set grid = ease - set MGRID = M36 -endif -if( $HRCODE == m25 ) then - set im = 1383 - set jm = 586 - set grid = ease - set MGRID = M25 + +if ( $isEASE > 0 ) then + + set grid = ease + + # im, jm must be consistent with values Fortran module EASE_conv + + if ( $EASEVERSION == EASEv1 ) then + + if( $HRCODE == M01 ) then + set im = 34668 + set jm = 14688 + endif + if( $HRCODE == M03 ) then + set im = 11556 + set jm = 4896 + endif + if( $HRCODE == M09 ) then + set im = 3852 + set jm = 1632 + endif + if( $HRCODE == M25 ) then + set im = 1383 + set jm = 586 + endif + if( $HRCODE == M36 ) then + set im = 963 + set jm = 408 + endif + + else if ( $EASEVERSION == EASEv2 ) then + + if( $HRCODE == M01 ) then + set im = 34704 + set jm = 14616 + endif + if( $HRCODE == M03 ) then + set im = 11568 + set jm = 4872 + endif + if( $HRCODE == M09 ) then + set im = 3856 + set jm = 1624 + endif + if( $HRCODE == M25 ) then + set im = 1388 + set jm = 584 + endif + if( $HRCODE == M36 ) then + set im = 964 + set jm = 406 + endif + + else + + echo " \!\!\!\! Invalid EASE grid version, stopping \!\!\!\! " + exit + + endif endif + + if( $MASKFILE == GEOS5_10arcsec_mask_freshwater-lakes.nc | $MASKFILE == GEOS5_10arcsec_mask.nc ) then @ NX = 43200 @ NY = 21600 @@ -936,7 +999,6 @@ ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/1440x1080 data/M cd data ln -s $input_dir CATCH - cd ../ if( -e CF${NC}x6C_${DATENAME}${IMO}x${POLENAME}${JMO}.stdout ) /bin/rm -f CF${NC}x6C_${DATENAME}${IMO}x${POLENAME}${JMO}.stdout @@ -1262,41 +1324,35 @@ endif # End Cube Test ####################################################################### if( $grid == ease ) then -echo $HRCODE - -if( $HRCODE == m25 ) then -set EVERSION = EASE -else -set EVERSION = EASEv2 -endif - -set RS = ${im}x${jm} -set IM = `echo ${im} | awk '{printf "%4.4i", $1}'` -set JM = `echo ${jm} | awk '{printf "%4.4i", $1}'` - -set BCNAME = SMAP_${EVERSION}_${MGRID} -set BCDIR = $EXPDIR/$OUTDIR/$BCNAME.scratch -set BCJOB = $BCDIR/$BCNAME.j -set THISGRID = SMAP-${EVERSION}-${MGRID} - -set nfiles = `find $EXPDIR -maxdepth 5 -name ${BCNAME}".j" | wc -l` -if( $nfiles >= 1 ) then - echo "" - echo "${C1} ----------------------------------------------------${CR}" - echo "${C1} Abort: ${CR}" - echo "${C2} This BCS run $BCDIR ${CR}" - echo "${C2} will create resolution dir already present: $BCNAME ${CR}" - echo "${C2} Please delete run dir and same resolution BCS files and resubmit" - echo "${C1} ----------------------------------------------------${CR}" - echo "" - exit -endif + echo $HRCODE + + set RS = ${im}x${jm} + set IM = `echo ${im} | awk '{printf "%4.4i", $1}'` + set JM = `echo ${jm} | awk '{printf "%4.4i", $1}'` + + set BCNAME = SMAP_${EASEVERSION}_${HRCODE} + set BCDIR = $EXPDIR/$OUTDIR/$BCNAME.scratch + set BCJOB = $BCDIR/$BCNAME.j + set THISGRID = SMAP-${EASEVERSION}-${HRCODE} + + set nfiles = `find $EXPDIR -maxdepth 5 -name ${BCNAME}".j" | wc -l` + if( $nfiles >= 1 ) then + echo "" + echo "${C1} ----------------------------------------------------${CR}" + echo "${C1} Abort: ${CR}" + echo "${C2} This BCS run $BCDIR ${CR}" + echo "${C2} will create resolution dir already present: $BCNAME ${CR}" + echo "${C2} Please delete run dir and same resolution BCS files and resubmit" + echo "${C1} ----------------------------------------------------${CR}" + echo "" + exit + endif -mkdir -p $BCDIR -mkdir -p $EXPDIR/$OUTDIR/logs + mkdir -p $BCDIR + mkdir -p $EXPDIR/$OUTDIR/logs -echo "${C1} Creating:${CR} ${C2}$BCJOB${CR}" -/bin/rm -f $BCJOB + echo "${C1} Creating:${CR} ${C2}$BCJOB${CR}" + /bin/rm -f $BCJOB cat << _EOF_ > $BCJOB #!/bin/csh -x @@ -1318,38 +1374,31 @@ cd data ln -s $input_dir CATCH cd ../ limit stacksize unlimited -if ( $EVERSION == EASEv2 ) then + +#if ( $EASEVERSION == EASEv2 ) then + ## This section was used to make Irrigated Tiles + ##if(${HRCODE} == M09 | ${HRCODE} == M36) then + ## bin/mkLandRaster.x -x ${NX} -y ${NY} -v -t ${NT} + ## bin/mkEASETilesParam.x -ease_label ${HRCODE} -pfaf_til T + ## bin/CombineRasters.x -f 0 -t 232000000 ${THISGRID} Pfafstetter > /dev/null + ## bin/CombineRasters.x -t 232000000 ${THISGRID} ${THISGRID}-Pfafstetter + ## /bin/mv til/${THISGRID}_${THISGRID}-Pfafstetter.til til/${THISGRID}_${THISGRID}-Pfafstetter.ind + ##endif +#endif + setenv MASKFILE ${MASKFILE} -## This section was used to make Irrigated Tiles -##if(${MGRID} == M09 | ${MGRID} == M36) then -## bin/mkLandRaster.x -x ${NX} -y ${NY} -v -t ${NT} -## bin/mkSMAPTilesPara_v2.x -smap_grid ${MGRID} -pfaf_til T -## bin/CombineRasters.x -f 0 -t 232000000 ${THISGRID} Pfafstetter > /dev/null -## bin/CombineRasters.x -t 232000000 ${THISGRID} ${THISGRID}-Pfafstetter -## /bin/mv til/${THISGRID}_${THISGRID}-Pfafstetter.til til/${THISGRID}_${THISGRID}-Pfafstetter.ind -##endif setenv OMP_NUM_THREADS 1 -bin/mkSMAPTilesPara_v2.x -smap_grid ${MGRID} -v $lbcsv -setenv OMP_NUM_THREADS ${NCPUS} -head -1 clsm/mkCatchParam.log > smap_cmd -chmod 755 smap_cmd -./smap_cmd -chmod 755 bin/create_README.csh -bin/create_README.csh -else -setenv MASKFILE ${MASKFILE} +bin/mkEASETilesParam.x -ease_label ${BCNAME} + setenv OMP_NUM_THREADS 1 -bin/mkSMAPTilesPara.x -smap_grid ${MGRID} +bin/mkCatchParam.x -g ${BCNAME} -v $lbcsv -x ${NX} -y ${NY} setenv OMP_NUM_THREADS ${NCPUS} -head -1 clsm/mkCatchParam.log > smap_cmd -chmod 755 smap_cmd -./smap_cmd +bin/mkCatchParam.x -g ${BCNAME} -v $lbcsv -x ${NX} -y ${NY} chmod 755 bin/create_README.csh bin/create_README.csh -endif /bin/mv clsm clsm.${IM}x${JM} -/bin/cp til/SMAP_${EVERSION}_${MGRID}_${RS}.til clsm.${IM}x${JM} +/bin/cp til/SMAP_${EASEVERSION}_${HRCODE}_${RS}.til clsm.${IM}x${JM} ##/bin/cp til//${THISGRID}_${THISGRID}-Pfafstetter.TIL clsm.${IM}x${JM} cd clsm.${IM}x${JM} @@ -1365,9 +1414,9 @@ cd clsm.${IM}x${JM} cd ../ -/bin/rm -rf SMAP_${EVERSION}_${MGRID} -/bin/mv clsm.${IM}x${JM} SMAP_${EVERSION}_${MGRID} - cd SMAP_${EVERSION}_${MGRID} +/bin/rm -rf SMAP_${EASEVERSION}_${HRCODE} +/bin/mv clsm.${IM}x${JM} SMAP_${EASEVERSION}_${HRCODE} + cd SMAP_${EASEVERSION}_${HRCODE} mkdir clsm /bin/mv ar.new \ bf.dat \ @@ -1397,8 +1446,8 @@ cd ../ clsm cd ../ -/bin/mv rst SMAP_${EVERSION}_${MGRID} -/bin/mv til SMAP_${EVERSION}_${MGRID} +/bin/mv rst SMAP_${EASEVERSION}_${HRCODE} +/bin/mv til SMAP_${EASEVERSION}_${HRCODE} cd ../../ /bin/mv $BCDIR/$BCNAME . @@ -1408,7 +1457,7 @@ cd ../../ /bin/rm -r $OUTDIR #mkdir -p IRRIG/$BCNAME/clsm IRRIG/$BCNAME/rst -#bin/mkIrrigTiles.x -x 43200 -y 21600 -b $BCNAME -t SMAP_${EVERSION}_${MGRID}_${RS}.til -r _${RS}_DE -p $IRRIGTHRES +#bin/mkIrrigTiles.x -x 43200 -y 21600 -b $BCNAME -t SMAP_${EASEVERSION}_${HRCODE}_${RS}.til -r _${RS}_DE -p $IRRIGTHRES #/bin/cp $BCNAME/${THISGRID}_${THISGRID}-Pfafstetter.TIL IRRIG/$BCNAME/. _EOF_ diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs.py new file mode 100755 index 000000000..b0466dd35 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_bcs.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +# +# source install/bin/g5_modules +# + +import sys +import argparse +import textwrap +import ruamel.yaml +from bcs_utils import * +from make_ease_bcs import * + +# Define the argument parser +def parse_args(): + + program_description = textwrap.dedent(f''' + Usage: + + Boundary Conditions (BCs) Package: + Creates surface tile and other model parameter input files + (.til [tiles], .rst [raster], and land parameters) for + combinations of atmospheric resolution, ocean resolution, + and land BCs version. + + STEP 1: Build the model. (GCM or GEOSldas) + STEP 2: cd [install-path]/bin + STEP 3: source g5_modules (for bash or zsh use g5_modules.[z]sh) + STEP 4: ./make_bcs.py + Answer the following interactive questions: + a) Select Land BCs version. + b) Select atmospheric resolution(s). + c) Select ocean resolution(s). + (Not relevant for land-only EASE-grid BCs.) + d) Enter BCs output directory. + e) Enter sponsor code for computing account. + + To skip the generation of land parameter files (ie, mkCatchParam.x), + use: ./make_bcsi.pc -noland + This option saves time when additional bcs are created that have + the exact same land parameters as an existing set of bcs because + the only difference between the two sets of bcs is the [non-tripolar] + ocean resolution. ''') + + parser = argparse.ArgumentParser(description='make bcs',epilog=program_description,formatter_class=argparse.RawDescriptionHelpFormatter) + # define a mutually exclusive group of arguments + group = parser.add_mutually_exclusive_group() + group.add_argument('-c', '--config_file', help='YAML config file') + + # Parse using parse_known_args so we can pass the rest to the remap scripts + # If config_file is used, then extra_args will be empty + # If flattened_yaml is used, then extra_args will be populated + args, extra_args = parser.parse_known_args() + return args, extra_args + +def main(): + + question_flag = False + config = '' + + # Parse the command line arguments from parse_args() capturing the arguments and the rest + command_line_args, extra_args = parse_args() + print(f'command_line_args: {command_line_args}') + config_yaml = command_line_args.config_file + + if config_yaml: + config = yaml_to_config(config_yaml) + else: + answers = ask_questions() + config = get_config_from_answers(answers) + + make_ease_bcs(config) + + +if __name__ == '__main__' : + exit("The python version of make_bcs is not yet ready for general use. Until further notice, please use csh script make_bcs") + main() + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_cube_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_cube_bcs.py new file mode 100755 index 000000000..52cf736d5 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_cube_bcs.py @@ -0,0 +1,422 @@ +#!/usr/bin/env python3 +# +# source install/bin/g5_modules + +import os +from bcs_utils import * + +cube_template = """ +#!/bin/csh -x + +#SBATCH --output={EXPDIR}/{OUTDIR}/logs/{BCNAME}.log +#SBATCH --error={EXPDIR}/{OUTDIR}/logs/{BCNAME}.err +#SBATCH --account={account} +#SBATCH --time=12:00:00 +#SBATCH --ntasks=28 +#SBATCH --job-name={BCNAME}.j +#SBATCH --constraint=sky + +cd {BCDIR} + +/bin/ln -s {bin_dir} +source bin/g5_modules +mkdir -p til rst data/MOM5 data/MOM6 clsm/plots +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/360x200 data/MOM5/360x200 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/720x410 data/MOM5/720x410 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/1440x1080 data/MOM5/1440x1080 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/72x36 data/MOM6/72x36 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/1440x1080 data/MOM6/1440x1080 + +cd data +ln -s {input_dir} CATCH + +cd ../ +if( -e CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout ) /bin/rm -f CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout +setenv MASKFILE {MASKFILE} +limit stacksize unlimited +bin/mkCubeFVRaster.x -x {NX} -y {NY} {NC} >/dev/null +bin/mkLandRaster.x -x {NX} -y {NY} -v -t {NT} + +if( {LATLON_OCEAN} == TRUE ) then + bin/mkLatLonRaster.x -x {NX} -y {NY} -b DE -p PE -t 0 {IMO} {JMO} >/dev/null + bin/CombineRasters.x -f 0 -t {NT} DE{IMO}xPE{JMO} Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} CF{NC}x6C DE{IMO}xPE{JMO}-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif + +if( {TRIPOL_OCEAN} == TRUE ) then + bin/mkMOMAquaRaster.x -x {NX} -y {NY} data/{MOM_VERSION}/{imo}x{jmo}/MAPL_Tripolar.nc > /dev/null + /bin/cp til/Pfafstetter.til til/Pfafstetter-ORIG.til + /bin/cp rst/Pfafstetter.rst rst/Pfafstetter-ORIG.rst + bin/FillMomGrid.x -f 0 -g Pfafstetter-M {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter data/{MOM_VERSION}/{imo}x{jmo}/MAPL_Tripolar.nc + /bin/mv til/Pfafstetter-M.til til/Pfafstetter.til + /bin/mv rst/Pfafstetter-M.rst rst/Pfafstetter.rst + bin/CombineRasters.x -f 0 -t {NT} {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} CF{NC}x6C {DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + bin/mk_runofftbl.x CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif +if( {CUBED_SPHERE_OCEAN} == TRUE ) then + bin/CombineRasters.x -f 0 -t {NT} CF{NC}x6C Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} CF{NC}x6C CF{NC}x6C-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_CF{NC}x6C-Pfafstetter -v {lbcsv} + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_CF{NC}x6C-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif + +/bin/mv clsm clsm.C{NC} +/bin/cp til/CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til clsm.C{NC} +if( {TRIPOL_OCEAN} == TRUE ) /bin/cp til/CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.TRN clsm.C{NC} +/bin/rm clsm.C{NC}/CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.file + +cd clsm.C{NC} + /bin/mv irrig.dat irrigation_{RC}.dat + /bin/mv vegdyn.data vegdyn_{RC}.dat + /bin/mv nirdf.dat nirdf_{RC}.dat + /bin/mv visdf.dat visdf_{RC}.dat + /bin/mv lai.dat lai_clim_{RC}.data + /bin/mv green.dat green_clim_{RC}.data + /bin/mv lnfm.dat lnfm_clim_{RC}.data + /bin/mv ndvi.dat ndvi_clim_{RC}.data + +/bin/rm -f sedfile +if( {CUBED_SPHERE_OCEAN} == TRUE ) then +cat > sedfile << EOF +s/{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter/OC{nc}x{nc6}-CF/g +s/CF{NC}x6C/PE{nc}x{nc6}-CF/g +EOF +sed -f sedfile CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til > tile.file +/bin/mv -f tile.file CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til +/bin/rm -f sedfile +else +cat > sedfile << EOF +s/CF{NC}x6C/PE{nc}x{nc6}-CF/g +s/{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter/PE{imo}x{jmo}-{DATENAME}/g +EOF +sed -f sedfile CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til > tile.file +/bin/mv -f tile.file CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til +/bin/rm -f sedfile +endif + +cd ../ + +/bin/rm -rf CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv clsm.C{NC} CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} + cd CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} + mkdir clsm + /bin/mv ar.new \ + bf.dat \ + ts.dat \ + catchment.def \ + cti_stats.dat \ + tau_param.dat \ + soil_param.dat \ + mosaic_veg_typs_fracs \ + soil_param.first \ + README \ + bad_sat_param.tiles \ + lai.* \ + AlbMap* \ + plots \ + CLM_veg_typs_fracs \ + mkCatchParam.log \ + CLM_NDep_SoilAlb_T2m \ + CLM4.5_abm_peatf_gdp_hdm_fc \ + catch_params.nc4 \ + catchcn_params.nc4 \ + country_and_state_code.data \ + clsm + cd ../ + +/bin/mv rst CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv til CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv {BCJOB} CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}/. +cd ../../ +/bin/mv {BCDIR}/{BCNAME} . +/bin/mv {BCJOB} {BCNAME} +/bin/mv {EXPDIR}/{OUTDIR}/logs {BCNAME}/. +/bin/mv {BCNAME}/clsm/mkCatchParam.log {BCNAME}/logs/mkCatchParam.log +/bin/rm -r {OUTDIR} +""" + +cube_template_1 = """ +#!/bin/csh -x + +#SBATCH --output={EXPDIR}/{OUTDIR}/logs/{BCNAME}.log +#SBATCH --error={EXPDIR}/{OUTDIR}/logs/{BCNAME}.err +#SBATCH --account={account} +#SBATCH --time=12:00:00 +#SBATCH --ntasks=1 +#SBATCH --job-name={BCNAME}.j +#SBATCH --constraint=sky + +cd {BCDIR} + +/bin/ln -s {bin_dir} +source bin/g5_modules +mkdir -p til rst data/MOM5 data/MOM6 clsm/plots +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/360x200 data/MOM5/360x200 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/720x410 data/MOM5/720x410 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/1440x1080 data/MOM5/1440x1080 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/72x36 data/MOM6/72x36 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/1440x1080 data/MOM6/1440x1080 + +cd data +ln -s {input_dir} CATCH +cd ../ + +if( -e CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout ) /bin/rm -f CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout +setenv MASKFILE {MASKFILE} +limit stacksize unlimited +bin/mkCubeFVRaster.x -x {NX} -y {NY} {NC} >/dev/null +bin/mkLandRaster.x -x {NX} -y {NY} -v -t {NT} +if( {LATLON_OCEAN} == TRUE ) then + bin/mkLatLonRaster.x -x {NX} -y {NY} -b DE -p PE -t 0 {IMO} {JMO} >/dev/null + bin/CombineRasters.x -f 0 -t {NT} DE{IMO}xPE{JMO} Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} CF{NC}x6C DE{IMO}xPE{JMO}-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} +endif + +if( {TRIPOL_OCEAN} == TRUE ) then + bin/mkMOMAquaRaster.x -x {NX} -y {NY} data/{MOM_VERSION}/{imo}x{jmo}/MAPL_Tripolar.nc > /dev/null + /bin/cp til/Pfafstetter.til til/Pfafstetter-ORIG.til + /bin/cp rst/Pfafstetter.rst rst/Pfafstetter-ORIG.rst + bin/FillMomGrid.x -f 0 -g Pfafstetter-M {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter data/{MOM_VERSION}/{imo}x{jmo}/MAPL_Tripolar.nc + /bin/mv til/Pfafstetter-M.til til/Pfafstetter.til + /bin/mv rst/Pfafstetter-M.rst rst/Pfafstetter.rst + bin/CombineRasters.x -f 0 -t {NT} {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} CF{NC}x6C {DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + bin/mk_runofftbl.x CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} +endif + +if( {CUBED_SPHERE_OCEAN} == TRUE ) then + bin/CombineRasters.x -f 0 -t {NT} CF{NC}x6C Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} CF{NC}x6C CF{NC}x6C-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_CF{NC}x6C-Pfafstetter -v {lbcsv} +endif +""" + +cube_template_2 = """ +#!/bin/csh -x + +#SBATCH --output={EXPDIR}/{OUTDIR}/logs/{BCNAME}-2.log +#SBATCH --error={EXPDIR}/{OUTDIR}/logs/{BCNAME}-2.err +#SBATCH --account={account} +#SBATCH --time=12:00:00 +#SBATCH --ntasks=28 +#SBATCH --job-name={BCNAME}-2.j +#SBATCH --constraint=sky + +cd {BCDIR} + +source bin/g5_modules + +setenv MASKFILE {MASKFILE} +limit stacksize unlimited + +if( {LATLON_OCEAN} == TRUE ) then + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif + +if( {TRIPOL_OCEAN} == TRUE ) then + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif + +if( {CUBED_SPHERE_OCEAN} == TRUE ) then + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g CF{NC}x6C_CF{NC}x6C-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif +/bin/mv clsm clsm.C{NC} +/bin/cp til/CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til clsm.C{NC} +if( {TRIPOL_OCEAN} == TRUE ) /bin/cp til/CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.TRN clsm.C{NC} +/bin/rm clsm.C{NC}/CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.file + +cd clsm.C{NC} + /bin/mv irrig.dat irrigation_{RC}.dat + /bin/mv vegdyn.data vegdyn_{RC}.dat + /bin/mv nirdf.dat nirdf_{RC}.dat + /bin/mv visdf.dat visdf_{RC}.dat + /bin/mv lai.dat lai_clim_{RC}.data + /bin/mv green.dat green_clim_{RC}.data + /bin/mv lnfm.dat lnfm_clim_{RC}.data + /bin/mv ndvi.dat ndvi_clim_{RC}.data + +/bin/rm -f sedfile +if( {CUBED_SPHERE_OCEAN} == TRUE ) then +cat > sedfile << EOF +s/{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter/OC{nc}x{nc6}-CF/g +s/CF{NC}x6C/PE{nc}x{nc6}-CF/g +EOF +sed -f sedfile CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til > tile.file +/bin/mv -f tile.file CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til +/bin/rm -f sedfile +else +cat > sedfile << EOF +s/CF{NC}x6C/PE{nc}x{nc6}-CF/g +s/{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter/PE{imo}x{jmo}-{DATENAME}/g +EOF +sed -f sedfile CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til > tile.file +/bin/mv -f tile.file CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til +/bin/rm -f sedfile +endif +cd ../ + +/bin/rm -rf CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv clsm.C{NC} CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} + cd CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} + mkdir clsm + /bin/mv ar.new \ + bf.dat \ + ts.dat \ + catchment.def \ + cti_stats.dat \ + tau_param.dat \ + soil_param.dat \ + mosaic_veg_typs_fracs \ + soil_param.first \ + README \ + bad_sat_param.tiles \ + lai.* \ + AlbMap* \ + plots \ + CLM_veg_typs_fracs \ + mkCatchParam.log \ + CLM_NDep_SoilAlb_T2m \ + CLM4.5_abm_peatf_gdp_hdm_fc \ + catch_params.nc4 \ + catchcn_params.nc4 \ + country_and_state_code.data \ + clsm + cd ../ +/bin/mv rst CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv til CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv {BCJOB}-2 CF{NC}x6C_{DATENAME}{IMO}x{POLENAME}{JMO}/. +cd ../../ +/bin/mv {BCDIR}/{BCNAME} . +/bin/mv {BCJOB}-2 {BCNAME} +/bin/mv {EXPDIR}/{OUTDIR}/logs {BCNAME}/. +/bin/mv {BCNAME}/clsm/mkCatchParam.log {BCNAME}/logs/mkCatchParam.log +/bin/rm -r {OUTDIR} +""" + +def make_cube_bcs(config): + bin_dir = os.getcwd() + if 'install/bin' not in bin_dir: + print(" please run this program in installed bin directory") + return + + grid_type = config['grid_type'] + if 'Cubed' not in grid_type : + print('This is not a Cubed-Sphere grid') + return + + resolution = config['resolution'] + + account = get_account() + IMO = '%04d'%config['imo'] + JMO = '%04d'%config['jmo'] + NC = '%04d'%config['im'] + imo = config['imo'] + jmo = config['imo'] + nc = config['im'] + + bcname = 'CF'+NC+'x6C_'+DATENAME+imo+'x'+POLENAME+jmo + if config['CUBED_SPHERE_OCEAN'] : + bcname = 'CF'+NC'x6_CF'+NC+'6C' + DATENAME = 'CF' + POLENAME = '' + imo = NC + jmo = '6C' + + + + now = datetime.now() + tmp_dir =now.strftime("%Y%m%d%H%M%S") + expdir = config['expdir'] + scratch_dir = expdir+ tmp_dir+'/'+bcname+'.scratch/' + log_dir = expdir+'/'+tmp_dir+'/logs/' + bcjob = scratch_dir+'/'+bcname+'.j' + + if os.path.exists(bcjob): + print('please remove the run temprory directory: ' + expdir+'/'+ tmp_dir) + return + + os.makedirs(scratch_dir) + os.makedirs(log_dir) + + job_script = ease_template.format(\ + account = account, \ + EXPDIR = config['expdir'], \ + OUTDIR = tmp_dir, \ + BCNAME = bcname, \ + bin_dir = bin_dir, \ + input_dir = config['inputdir'], \ + BCJOB = bcjob, \ + EASEVERSION = grid_type, \ + HRCODE = resolution, \ + IM = config['im'], \ + JM = config['jm'], \ + MASKFILE = config['MASKFILE'], \ + lbcsv = config['lbcsv'], \ + NX = config['NX'], \ + NY = config['NY'], \ + RS = RS,\ + BCDIR = scratch_dir, \ + NCPUS = config['NCPUS']) + + + ease_job = open(bcjob,'wt') + ease_job.write(job_script) + ease_job.close() + + interactive = os.getenv('SLURM_JOB_ID', default = None) + if ( interactive ) : + print('interactive mode\n') + ntasks = os.getenv('SLURM_NTASKS', default = None) + if ( not ntasks): + nnodes = int(os.getenv('SLURM_NNODES', default = '1')) + ncpus = int(os.getenv('SLURM_CPUS_ON_NODE', default = '28')) + subprocess.call(['chmod', '755', bcjob]) + print(bcjob+ ' 1>' + log_name + ' 2>&1') + os.system(bcjob + ' 1>' + log_name+ ' 2>&1') + else: + print("sbatch " + bcjob +"\n") + subprocess.call(['sbatch', bcjob]) + + print( "cd " + bin_dir) + os.chdir(bin_dir) + + print(job_script) + +if __name__ == "__main__": + + answers = ask_questions() + config = get_config_from_answers(answers) + print("make_ease_bcs") + make_ease_bcs(config) + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_ease_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_ease_bcs.py new file mode 100755 index 000000000..524bfe7b6 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_ease_bcs.py @@ -0,0 +1,178 @@ +#!/usr/bin/env python3 +# +# source install/bin/g5_modules + +import os +from bcs_utils import * + + +ease_template = """#!/bin/csh -x + +#SBATCH --output={EXPDIR}/{OUTDIR}/logs/{BCNAME}.log +#SBATCH --error={EXPDIR}/{OUTDIR}/logs/{BCNAME}.err +#SBATCH --account={account} +#SBATCH --time=12:00:00 +#SBATCH --ntasks=28 +#SBATCH --job-name={BCNAME}.j +#SBATCH --constraint=sky + +cd {BCDIR} + +/bin/ln -s {bin_dir} +source bin/g5_modules +mkdir -p til rst data/MOM5 data/MOM6 clsm/plots +cd data +ln -s {input_dir} CATCH +cd ../ +limit stacksize unlimited + +setenv MASKFILE {MASKFILE} +setenv OMP_NUM_THREADS 1 +bin/mkEASETilesParam.x -ease_label {BCNAME} +setenv OMP_NUM_THREADS 1 +bin/mkCatchParam.x -g {BCNAME} -v {lbcsv} -x {NX} -y {NY} +setenv OMP_NUM_THREADS {NCPUS} +bin/mkCatchParam.x -g {BCNAME} -v {lbcsv} -x {NX} -y {NY} +chmod 755 bin/create_README.csh +bin/create_README.csh + +/bin/mv clsm clsm.{IM}x{JM} +/bin/cp til/SMAP_{EASEVERSION}_{HRCODE}_{RS}.til clsm.{IM}x{JM} + +cd clsm.{IM}x{JM} + /bin/mv irrig.dat irrigation_{RS}_DE.dat + /bin/mv vegdyn.data vegdyn_{RS}_DE.dat + /bin/mv nirdf.dat nirdf_{RS}_DE.dat + /bin/mv visdf.dat visdf_{RS}_DE.dat + /bin/mv lai.dat lai_clim_{RS}_DE.data + /bin/mv green.dat green_clim_{RS}_DE.data + /bin/mv lnfm.dat lnfm_clim_{RS}_DE.data + /bin/mv ndvi.dat ndvi_clim_{RS}_DE.data + + +cd ../ +/bin/rm -rf SMAP_{EASEVERSION}_{HRCODE} +/bin/mv clsm.{IM}x{JM} SMAP_{EASEVERSION}_{HRCODE} + cd SMAP_{EASEVERSION}_{HRCODE} + mkdir clsm + /bin/mv ar.new \ + bf.dat \ + ts.dat \ + catchment.def \ + cti_stats.dat \ + tau_param.dat \ + soil_param.dat \ + mosaic_veg_typs_fracs \ + soil_param.first \ + bad_sat_param.tiles \ + README \ + lai.* \ + AlbMap* \ + g5fmt \ + vegetation.hst2 \ + pfaf_fractions.dat \ + plots \ + CLM_veg_typs_fracs \ + mkCatchParam.log \ + Grid2Catch_TransferData.nc \ + CLM_NDep_SoilAlb_T2m \ + CLM4.5_abm_peatf_gdp_hdm_fc \ + catch_params.nc4 \ + catchcn_params.nc4 \ + country_and_state_code.data \ + clsm + cd ../ +/bin/mv rst SMAP_{EASEVERSION}_{HRCODE} +/bin/mv til SMAP_{EASEVERSION}_{HRCODE} + +cd ../../ +/bin/mv {BCDIR}/{BCNAME} . +/bin/mv {BCJOB} {BCNAME} +/bin/mv {EXPDIR}/{OUTDIR}/logs {BCNAME}/. +/bin/mv {BCNAME}/clsm/mkCatchParam.log {BCNAME}/logs/mkCatchParam.log +/bin/rm -r {OUTDIR} + +""" + +def make_ease_bcs(config): + bin_dir = os.getcwd() + if 'install/bin' not in bin_dir: + print(" please run this program in installed bin directory") + return + + grid_type = config['grid_type'] + if 'EASEv' not in grid_type : + print('This is not a EASE grid') + return + + resolution = config['resolution'] + + EASElabel = 'SMAP_'+grid_type+'_'+ resolution + now = datetime.now() + tmp_dir =now.strftime("%Y%m%d%H%M%S") + expdir = config['expdir'] + scratch_dir = expdir+ tmp_dir+'/'+EASElabel+'.scratch/' + log_dir = expdir+'/'+tmp_dir+'/logs/' + bcjob = scratch_dir+'/'+EASElabel+'.j' + if os.path.exists(bcjob): + print('please remove the run temprory directory: ' + expdir+'/'+ tmp_dir) + return + + os.makedirs(scratch_dir) + os.makedirs(log_dir) + + account = get_account() + ims = '%04d'%config['im'] + jms = '%04d'%config['jm'] + RS = ims+'x'+jms + job_script = ease_template.format(\ + account = account, \ + EXPDIR = config['expdir'], \ + OUTDIR = tmp_dir, \ + BCNAME = EASElabel, \ + bin_dir = bin_dir, \ + input_dir = config['inputdir'], \ + BCJOB = bcjob, \ + EASEVERSION = grid_type, \ + HRCODE = resolution, \ + IM = config['im'], \ + JM = config['jm'], \ + MASKFILE = config['MASKFILE'], \ + lbcsv = config['lbcsv'], \ + NX = config['NX'], \ + NY = config['NY'], \ + RS = RS,\ + BCDIR = scratch_dir, \ + NCPUS = config['NCPUS']) + + + ease_job = open(bcjob,'wt') + ease_job.write(job_script) + ease_job.close() + + interactive = os.getenv('SLURM_JOB_ID', default = None) + if ( interactive ) : + print('interactive mode\n') + ntasks = os.getenv('SLURM_NTASKS', default = None) + if ( not ntasks): + nnodes = int(os.getenv('SLURM_NNODES', default = '1')) + ncpus = int(os.getenv('SLURM_CPUS_ON_NODE', default = '28')) + subprocess.call(['chmod', '755', bcjob]) + print(bcjob+ ' 1>' + log_name + ' 2>&1') + os.system(bcjob + ' 1>' + log_name+ ' 2>&1') + else: + print("sbatch " + bcjob +"\n") + subprocess.call(['sbatch', bcjob]) + + print( "cd " + bin_dir) + os.chdir(bin_dir) + + print(job_script) + +if __name__ == "__main__": + + answers = ask_questions() + config = get_config_from_answers(answers) + print("make_ease_bcs") + make_ease_bcs(config) + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_latlon_bcs.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_latlon_bcs.py new file mode 100755 index 000000000..0591c1074 --- /dev/null +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/make_latlon_bcs.py @@ -0,0 +1,226 @@ +#!/usr/bin/env python3 +# +# source install/bin/g5_modules + +import os +from bcs_utils import * + + +ease_template = """#!/bin/csh -x + +#SBATCH --output={EXPDIR}/{OUTDIR}/logs/{BCNAME}.log +#SBATCH --error={EXPDIR}/{OUTDIR}/logs/{BCNAME}.err +#SBATCH --account={account} +#SBATCH --time=12:00:00 +#SBATCH --ntasks=28 +#SBATCH --job-name={BCNAME}.j +#SBATCH --constraint=sky + +cd {BCDIR} + +source bin/g5_modules +/bin/ln -s {bin_dir} +mkdir -p til rst data/MOM5 data/MOM6 clsm/plots +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/360x200 data/MOM5/360x200 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/720x410 data/MOM5/720x410 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM5/1440x1080 data/MOM5/1440x1080 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/72x36 data/MOM6/72x36 +ln -s /discover/nobackup/projects/gmao/ssd/aogcm/ocean_bcs/MOM6/1440x1080 data/MOM6/1440x1080 + +cd data +ln -s {input_dir} CATCH +cd ../ + +if( -e DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}.stdout ) /bin/rm -f DC{IM}xPC{JM}_{DATENAME}{IMO}{POLENAME}{JMO}.stdout +setenv MASKFILE {MASKFILE} +limit stacksize unlimited +bin/mkLatLonRaster.x -x {NX} -y {NY} -t -1 {IM} {JM} >/dev/null +bin/mkLandRaster.x -x {NX} -y {NY} -v -t {NT} + +if( {LATLON_OCEAN} == TRUE ) then + bin/mkLatLonRaster.x -x {NX} -y {NY} -b DE -p PE -t 0 {IMO} {JMO} >/dev/null + bin/CombineRasters.x -f 0 -t {NT} DE{IMO}xPE{JMO} Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} DC{IM}xPC{JM} DE{IMO}xPE{JMO}-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g DC{IM}xPC{JM}_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g DC{IM}xPC{JM}_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif +if( {TRIPOL_OCEAN} == TRUE ) then + bin/mkMOMAquaRaster.x -x {NX} -y {NY} data/{MOM_VERSION}/{imo}x{jmo}/MAPL_Tripolar.nc > /dev/null + /bin/cp til/Pfafstetter.til til/Pfafstetter-ORIG.til + /bin/cp rst/Pfafstetter.rst rst/Pfafstetter-ORIG.rst + bin/FillMomGrid.x -f 0 -g Pfafstetter-M {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter data/{MOM_VERSION}/{imo}x{jmo}/MAPL_Tripolar.nc + /bin/mv til/Pfafstetter-M.til til/Pfafstetter.til + /bin/mv rst/Pfafstetter-M.rst rst/Pfafstetter.rst + bin/CombineRasters.x -f 0 -t {NT} {DATENAME}{IMO}x{POLENAME}{JMO} Pfafstetter >/dev/null + bin/CombineRasters.x -t {NT} DC{IM}xPC{JM} {DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + bin/mk_runofftbl.x DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter + setenv OMP_NUM_THREADS 1 + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g DE{IMO}xPE{JMO}_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + setenv OMP_NUM_THREADS {NCPUS} + if ({SKIPLAND} != YES) bin/mkCatchParam.x -x {NX} -y {NY} -g DE{IMO}xPE{JMO}_DE{IMO}xPE{JMO}-Pfafstetter -v {lbcsv} + chmod 755 bin/create_README.csh + bin/create_README.csh +endif + +/bin/mv clsm clsm.{IM}x{JM} +/bin/cp til/DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til clsm.{IM}x{JM} +if( {TRIPOL_OCEAN} == TRUE ) /bin/cp til/DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.TRN clsm.{IM}x{JM} +/bin/rm clsm.{IM}x{JM}/DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.file + +cd clsm.{IM}x{JM} + /bin/mv irrig.dat irrigation_{RS}_DC.dat + /bin/mv vegdyn.data vegdyn_{RS}_DC.dat + /bin/mv nirdf.dat nirdf_{RS}_DC.dat + /bin/mv visdf.dat visdf_{RS}_DC.dat + /bin/mv lai.dat lai_clim_{RS}_DC.data + /bin/mv green.dat green_clim_{RS}_DC.data + /bin/mv lnfm.dat lnfm_clim_{RS}_DC.data + /bin/mv ndvi.dat ndvi_clim_{RS}_DC.data +/bin/rm -f sedfile +cat > sedfile << EOF +s/DC{IM}xPC{JM}/PC{im}x{jm}-DC/g +s/{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter/PE{imo}x{jmo}-{DATENAME}/g +EOF +sed -f sedfile DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til > tile.file +/bin/mv -f tile.file DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO}-Pfafstetter.til +/bin/rm -f sedfile + +cd ../ + +/bin/rm -rf DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv clsm.{IM}x{JM} DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO} + cd DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO} + mkdir clsm + /bin/mv ar.new \ + bf.dat \ + ts.dat \ + catchment.def \ + cti_stats.dat \ + tau_param.dat \ + soil_param.dat \ + mosaic_veg_typs_fracs \ + soil_param.first \ + README \ + bad_sat_param.tiles \ + lai.* \ + AlbMap* \ + plots \ + CLM_veg_typs_fracs \ + mkCatchParam.log \ + CLM_NDep_SoilAlb_T2m \ + CLM4.5_abm_peatf_gdp_hdm_fc \ + catch_params.nc4 \ + catchcn_params.nc4 \ + country_and_state_code.data \ + clsm + cd ../ + +/bin/mv rst DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO} +/bin/mv til DC{IM}xPC{JM}_{DATENAME}{IMO}x{POLENAME}{JMO} + +cd ../../ +/bin/mv {BCDIR}/{BCNAME} . +/bin/mv {BCJOB} {BCNAME} +/bin/mv {EXPDIR}/{OUTDIR}/logs {BCNAME}/. +/bin/mv {BCNAME}/clsm/mkCatchParam.log {BCNAME}/logs/mkCatchParam.log +/bin/rm -r {OUTDIR} + +""" + +def make_cube_bcs(config): + bin_dir = os.getcwd() + if 'install/bin' not in bin_dir: + print(" please run this program in installed bin directory") + return + + grid_type = config['grid_type'] + if 'Cubed' not in grid_type : + print('This is not a Cubed-Sphere grid') + return + + resolution = config['resolution'] + + account = get_account() + imo = '%04d'%config['imo'] + jmo = '%04d'%config['jmo'] + NC = '%04d'%config['im'] + + bcname = 'CF'+NC+'x6C_'+DATENAME+imo+'x'+POLENAME+jmo + if config['CUBED_SPHERE_OCEAN'] : + bcname = 'CF'+NC'x6_CF'+NC+'6C' + DATENAME = 'CF' + POLENAME = '' + imo = NC + jmo = '6C' + + + + now = datetime.now() + tmp_dir =now.strftime("%Y%m%d%H%M%S") + expdir = config['expdir'] + scratch_dir = expdir+ tmp_dir+'/'+bcname+'.scratch/' + log_dir = expdir+'/'+tmp_dir+'/logs/' + bcjob = scratch_dir+'/'+bcname+'.j' + + if os.path.exists(bcjob): + print('please remove the run temprory directory: ' + expdir+'/'+ tmp_dir) + return + + os.makedirs(scratch_dir) + os.makedirs(log_dir) + + job_script = ease_template.format(\ + account = account, \ + EXPDIR = config['expdir'], \ + OUTDIR = tmp_dir, \ + BCNAME = bcname, \ + bin_dir = bin_dir, \ + input_dir = config['inputdir'], \ + BCJOB = bcjob, \ + EASEVERSION = grid_type, \ + HRCODE = resolution, \ + IM = config['im'], \ + JM = config['jm'], \ + MASKFILE = config['MASKFILE'], \ + lbcsv = config['lbcsv'], \ + NX = config['NX'], \ + NY = config['NY'], \ + RS = RS,\ + BCDIR = scratch_dir, \ + NCPUS = config['NCPUS']) + + + ease_job = open(bcjob,'wt') + ease_job.write(job_script) + ease_job.close() + + interactive = os.getenv('SLURM_JOB_ID', default = None) + if ( interactive ) : + print('interactive mode\n') + ntasks = os.getenv('SLURM_NTASKS', default = None) + if ( not ntasks): + nnodes = int(os.getenv('SLURM_NNODES', default = '1')) + ncpus = int(os.getenv('SLURM_CPUS_ON_NODE', default = '28')) + subprocess.call(['chmod', '755', bcjob]) + print(bcjob+ ' 1>' + log_name + ' 2>&1') + os.system(bcjob + ' 1>' + log_name+ ' 2>&1') + else: + print("sbatch " + bcjob +"\n") + subprocess.call(['sbatch', bcjob]) + + print( "cd " + bin_dir) + os.chdir(bin_dir) + + print(job_script) + +if __name__ == "__main__": + + answers = ask_questions() + config = get_config_from_answers(answers) + print("make_ease_bcs") + make_ease_bcs(config) + diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 index e163fb2d5..db921e3e1 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkCatchParam.F90 @@ -4,14 +4,12 @@ PROGRAM mkCatchParam ! ! !ARGUMENTS: ! -! Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE" +! Usage = "mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV " ! -x: Size of longitude dimension of input raster. DEFAULT: 8640 ! -y: Size of latitude dimension of input raster. DEFAULT: 4320 ! -b: position of the dateline in the first box. DEFAULT: DC ! -g: Gridname (name of the .til or .rst file without file extension) ! -v: LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09) -! -e: EASE : This is optional if catchment.def file is available already or -! the til file format is pre-Fortuna-2. ! ! ! This program is good to generate @@ -20,7 +18,7 @@ PROGRAM mkCatchParam ! ! Sarith Mahanama - March 23, 2012 ! Email: sarith.p.mahanama@nasa.gov - + use EASE_conv use rmTinyCatchParaMod use process_hres_data ! use module_irrig_params, ONLY : create_irrig_params @@ -40,13 +38,12 @@ PROGRAM mkCatchParam character*1 :: opt character*7 :: PEATSOURCE = 'GDLHWSD' character*3 :: VEGZSOURCE = 'D&S' - character*4 :: EASE =' ' character*2 :: DL ='DC' integer :: II, JJ, Type integer :: I, J, command_argument_count, nxt real*8 :: dx, dy, lon0 logical :: regrid - character(len=400), dimension (8) :: Usage + character(len=400), dimension (6) :: Usage character*128 :: Grid2 character*2 :: poles character*128 :: GridNameR = '' @@ -61,6 +58,8 @@ PROGRAM mkCatchParam character*200 :: fname_tmp, fname_tmp2, fname_tmp3, fname_tmp4 integer :: N_tile logical :: process_snow_albedo = .false. + character(len=10) :: nc_string, nr_string + integer :: nc_ease, nr_ease ! --------- VARIABLES FOR *OPENMP* PARALLEL ENVIRONMENT ------------ ! @@ -103,14 +102,12 @@ PROGRAM mkCatchParam ! call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH') ! call execute_command_line('cd ..') - USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV -e EASE " - USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " - USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " - USAGE(4) =" -g: Gridname (name of the .til or .rst file without file extension) " - USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " - USAGE(6) =" -e: EASE : This is optional if catchment.def file is available already or " - USAGE(7) =" the til file format is pre-Fortuna-2. " - USAGE(8) =" -v LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08, v09) " + USAGE(1) ="Usage: mkCatchParam -x nx -y ny -g Gridname -b DL -v LBCSV " + USAGE(2) =" -x: Size of longitude dimension of input raster. DEFAULT: 8640 " + USAGE(3) =" -y: Size of latitude dimension of input raster. DEFAULT: 4320 " + USAGE(4) =" -g: Gridname (name of the .til or .rst file without file extension) " + USAGE(5) =" -b: Position of the dateline in the first grid box (DC or DE). DEFAULT: DC " + USAGE(6) =" -v LBCSV : Land bcs version (F25, GM4, ICA, NL3, NL4, NL5, v06, v07, v08 v09 ) " ! Process Arguments !------------------ @@ -159,9 +156,6 @@ PROGRAM mkCatchParam call init_bcs_config (trim(LBCSV)) ! get bcs details from version string case ('b') DL = trim(arg) - case ('e') - EASE = trim(arg) - if(EASE=='EASE') ease_grid=.true. case default do j = 1,size(usage) print "(sp,a100)", Usage(j) @@ -174,17 +168,26 @@ PROGRAM mkCatchParam call get_environment_variable ("MASKFILE" ,MaskFile ) - if(trim(Gridname) == '') then + if(trim(Gridname) == '') then write (log_file,'(a)')'Unable to create parameters without til/rst files.... !' stop endif regrid = nc/=i_raster .or. nr/=j_raster + + if (index(Gridname,'EASEv') /=0) then + ! here gridname has alias EASELabel + call ease_extent(Gridname, nc_ease, nr_ease ) + write(nc_string, '(i0)') nc_ease + write(nr_string, '(i0)') nr_ease + gridname = trim(Gridname)//'_'//trim(nc_string)//'x'//trim(nr_string) + ease_grid = .true. + endif if(index(Gridname,'Pfaf.notiny')/=0) then GridnameR='clsm/'//trim(Gridname) GridnameT='clsm/'//trim(Gridname) - else + else GridnameR='rst/'//trim(Gridname) GridnameT='til/'//trim(Gridname) endif diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkEASETilesParam.F90 similarity index 68% rename from GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 rename to GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkEASETilesParam.F90 index 55c1abc37..db0ae03db 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara_v2.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkEASETilesParam.F90 @@ -1,22 +1,29 @@ #include "Raster.h" -PROGRAM mkSMAPTilesPara_v2 - -! This program constructs land and lake tiles for the SMAP-EASE-M09 and M36 grids (just set MGRID) -! for CLSM implementation. -! f90 -c create_smap_tiles.f90 -! f90 -c smapconv.f -! f90 -o create_smap_tiles create_smap_tiles.o smapconv.o -! - use easeV2_conv - use rmTinyCatchParaMod +PROGRAM mkEASETilesParam + + ! This program constructs land and lake tiles for EASE grid tile spaces such as + ! those used for the SMAP Level-4 products and other offline projects + + ! This program resulted from the merger and cleanup of mkSMAPTilesPara.F90 + ! and mkSMAPTilesPara_v2.F90 in September 2022. + ! Before the merger and cleanup, the EASE grid parameters were hard-coded here. + ! For EASEv2 M25, the outdated scale value was used here. + ! This is program is renamed to mkEASETileParam from mkSMAPTilesPara_v2 + ! + ! - wjiang + reichle, 21 Sep 2022 + + + use EASE_conv + use rmTinyCatchParaMod, only : i_raster, j_raster, SRTM_maxcat + use rmTinyCatchParaMod, only : RegridRaster, RegridRaster1, RegridRasterReal use process_hres_data use MAPL_SortMod use MAPL_ConstantsMod use LogRectRasterizeMod + use netcdf implicit none - character*5 :: LBCSV = 'UNDEF' integer i,j,ig,jg,i0,iop,n,d1,d2,j1,j2,i1,i2,ix, jx,icount,pcount integer :: NC = i_raster, NR = j_raster, NT = 16330000, ND = 10000, ND_raster = 10000 @@ -30,17 +37,18 @@ PROGRAM mkSMAPTilesPara_v2 INTEGER, ALLOCATABLE, DIMENSION (:) :: density, loc_int logical, dimension (:), allocatable :: unq_mask integer, pointer , dimension (:,:) :: subset - integer, pointer , dimension (:) :: subset1, subset_smap + integer, pointer , dimension (:) :: subset1 real, pointer , dimension (:) :: subset2 integer :: dx_esa, dy_esa, NBINS, NPLUS integer*8, allocatable, dimension (:) :: SRTM_catid + real(kind=8),allocatable, dimension(:):: SRTM_catid_r8 integer,allocatable, dimension (:,:), target :: tileid_index,catid_index integer,allocatable, dimension (:,:) :: catid, iaster integer,allocatable, dimension (:) :: land_id,water_id,ice_id integer,allocatable, dimension (:) :: my_land, all_id - real, allocatable, dimension (:) :: smap_grid_area,tile_area,SRTM_CatchArea + real, allocatable, dimension (:) :: ease_grid_area,tile_area,SRTM_CatchArea integer*1,allocatable, dimension (:,:) :: veg, i2aster real*4, dimension (:,:), allocatable :: q0,raster REAL, dimension (:), allocatable :: tile_ele, tile_area_land @@ -49,137 +57,107 @@ PROGRAM mkSMAPTilesPara_v2 integer l,imn,imx,jmn,jmx,mval,l_index,i_index,w_index,typ,pfaf,cindex integer :: LakeType, IceType, OceanType character(3) :: easegrid - real :: clat, clon, r_smap, s_smap, smap_convert, da + real :: clat, clon, r_ease, s_ease, da real :: fr_gcm - integer :: ind_col, ind_row, status, ncid, nciv,nland_cells, DOM_INDX + integer :: ind_col, ind_row, status, ncid, varid, nciv,nland_cells, DOM_INDX REAL (kind=8), PARAMETER :: RADIUS=6378137.0,pi=3.14159265358979323846 character*100 :: veg_class (12) - character*5 :: MGRID character*100 :: gfile,gtopo30 - integer :: nc_smap,nr_smap, N_args, command_argument_count - real :: EASE_grid_area, CELL_km + integer :: nc_ease,nr_ease, N_args, command_argument_count REAL :: dx,dy,d2r,lats,mnx,mxx,mny,mxy,sum1,sum2,jgv, VDUM,pix_area - character(40) :: arg, EASElabel + character(40) :: arg, EASElabel_ + character(len=:), allocatable :: EASElabel character*200 :: tmpstring, tmpstring1, tmpstring2 logical :: regrid = .false. character*128 :: MaskFile logical :: pfaf_til = .false. character*1 :: PF - - include 'netcdf.inc' + character(len=6) :: EASE_Version + character(len=10) :: nc_string, nr_string + character(128) :: usage1, usage2 + + ! -------------------------------------------------------------------------------------- + + usage1 = 'USAGE : bin/mkEASETilesParam.x -ease_label EASELabel ' + usage2 = ' where EASELabel = *EASEv[x]_M[yy]*, x={1,2}, yy={01,03,09,25,36}' N_args = command_argument_count() if(N_args < 1) then - print *,'USAGE : bin/mkSMAPTiles -smap_grid MXX' - print *,'Allowed SMAP grids are: M01 M03 M09 M25 M36' + print *,trim(usage1) + print *,trim(usage2) stop end if i=0 - do while ( i < N_args ) i = i+1 call get_command_argument(i,arg) - if ( trim(arg) == '-smap_grid' ) then + if ( trim(arg) == '-ease_label' ) then i = i+1 - call get_command_argument(i,MGRID) + call get_command_argument(i,EASELabel_) - elseif ( trim(arg) == '-pfaf_til' ) then - i = i+1 - call get_command_argument(i,PF) - if (PF == 'T') pfaf_til = .true. + ! WY noted: this may be used in the future for irrigation tiles + !elseif ( trim(arg) == '-pfaf_til' ) then + ! i = i+1 + ! call get_command_argument(i,PF) + ! if (PF == 'T') pfaf_til = .true. - elseif ( trim(arg) == '-v' ) then - i = i+1 - call get_command_argument(i,LBCSV) - else ! stop for any other arguments - - print *,'USAGE : bin/mkSMAPTiles -smap_grid MXX -pfaf_til T' - print *,'Allowed SMAP grids are: M09 M36 Ml' + print *,trim(usage1) + print *,trim(usage2) stop - endif end do - call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH') - call execute_command_line('cd ..') + ! WY note, remove this verification. There can be all combination + ! if (MGRID /= 'M25' .and. EASE_version == 'EASEv1') then + ! stop ("EASEv1 only supports M25") + ! endif + ! WY noted: should do it in the script that calls this program + !call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH') + !call execute_command_line('cd ..') - ! Setting SMAP Grid specifications + ! Setting EASE Grid specifications ! -------------------------------- - - if (trim(MGRID) == 'M09') then - - CELL_km = 9.008055210146 ! nominal cell size in kilometers - nc_smap = 3856 - nr_smap = 1624 - gfile = 'SMAP_EASEv2_'//trim(MGRID)//'_3856x1624' - EASE_grid_area = CELL_km*CELL_km - EASElabel = 'SMAP-EASEv2-M09' - - elseif(trim(MGRID) == 'M36') then - - CELL_km = 36.032220840584 ! nominal cell size in kilometers - nc_smap = 964 - nr_smap = 406 - gfile = 'SMAP_EASEv2_'//trim(MGRID)//'_964x406' - EASE_grid_area = CELL_km*CELL_km - EASElabel = 'SMAP-EASEv2-M36' - - elseif(trim(MGRID) == 'M25') then - - CELL_km = 25.0252600081 ! nominal cell size in kilometers - nc_smap = 1388 - nr_smap = 584 - gfile = 'SMAP_EASEv2_M25_1388x584' - EASE_grid_area = CELL_km*CELL_km - - else if (trim(MGRID) .eq. 'M03') then ! SMAP 3 km grid - CELL_km = 3.0026850700487 ! nominal cell size in kilometers - nc_smap = 11568 - nr_smap = 4872 - gfile = 'SMAP_EASEv2_M03_11568x4872' - EASE_grid_area = CELL_km*CELL_km + + EASElabel = trim(EASELabel_) + + call ease_extent(EASELabel, nc_ease, nr_ease ) + write(nc_string, '(i0)') nc_ease + write(nr_string, '(i0)') nr_ease + gfile = trim(EASElabel)//'_'//trim(nc_string)//'x'//trim(nr_string) + + if (index(EASELabel,'M03') /=0) then ! EASE 3 km grid regrid = .true. NC = 21600 NR = 10800 NT = 500000000 - - else if (trim(MGRID) .eq. 'M01') then ! SMAP 1 km grid - CELL_km = 1.00089502334956 ! nominal cell size in kilometers - nc_smap = 34704 - nr_smap = 14616 - gfile = 'SMAP_EASEv2_M01_34704x14616' - EASE_grid_area = CELL_km*CELL_km + endif + if (index(EASELabel,'M03') /=0) then ! EASE 1 km grid regrid = .true. NC = 43200 - NR = 21600 + NR = 21600 NT = 1500000000 - - else ! - - print *,'Unknown SMAP Grid stopping..' - stop - endif allocate(land_id (1:NT)) allocate(water_id (1:NT)) allocate(ice_id (1:NT)) - land_id = 0 - water_id = 0 - ice_id = 0 - OceanType = 0 - IceType =11 - LakeType =10 - - ND = 10*10**(nint(log10(1.*nr_smap))) + + land_id = 0 + water_id = 0 + ice_id = 0 + OceanType = 0 + IceType = 11 + LakeType = 10 + + ND = 10*10**(nint(log10(1.*nr_ease))) ! Check for the 10 arc-sec MaskFile ! ----------------------------------- @@ -188,20 +166,20 @@ PROGRAM mkSMAPTilesPara_v2 print *, 'Using MaskFile ', trim(MaskFile) - if(pfaf_til) then + ! This section was used to make Irrigated Tiles + !if(pfaf_til) then - nc = 43200 ! Number of rows in raster file - nr = 21600 - call mkEASEv2Raster + ! nc = 43200 ! Number of rows in raster file + ! nr = 21600 + ! call mkEASEv2Raster !else - ! This section was used to make Irrigated Tiles ! if((trim(MGRID) == 'M09').or.(trim(MGRID) == 'M36'))call write_tilfile - endif + !endif if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then ! New ESA (Veg) + SRTM (catchments) based mask file - ! is overlaid on SMAP + ! is overlaid on the EASE grid ! ------------------------------------------------- nc = 43200 ! Number of rows in raster file @@ -213,6 +191,7 @@ PROGRAM mkSMAPTilesPara_v2 allocate(tileid_index(1:nc,1:nr)) allocate(SRTM_catid (1:SRTM_maxcat+2)) + allocate(SRTM_catid_r8(1:SRTM_maxcat+2), source = 0.d0) allocate(catid_index (1:nc,1:nr)) allocate(veg (1:nc,1:nr)) allocate(geos_msk (1:nc_esa,1:dy_esa)) @@ -238,25 +217,26 @@ PROGRAM mkSMAPTilesPara_v2 catid_index = 0 veg = 0 - status = NF_OPEN ('data/CATCH/GEOS5_10arcsec_mask.nc', NF_NOWRITE, ncid) - status = NF_GET_VARA_INT64 (ncid,3,(/1/),(/SRTM_maxcat/),SRTM_catid(1:SRTM_maxcat)) ! Read pfafstetter IDs + status = NF90_OPEN ('data/CATCH/GEOS5_10arcsec_mask.nc', NF90_NOWRITE, ncid) + status = nf90_inq_varid(ncid, name='PfafID', varid=varid) + status = nf90_get_var(ncid, varid, SRTM_catid_r8, (/1/),(/SRTM_maxcat/)) if(status /=0) then - PRINT *, NF_STRERROR(STATUS) - print *, 'Problem with NF_OPEN',trim(MaskFile) + PRINT *, NF90_STRERROR(STATUS) + print *, 'Problem with NF90_OPEN',trim(MaskFile) endif - + SRTM_catid = int8(SRTM_catid_r8) SRTM_catid (SRTM_maxcat + 1) = 190000000 SRTM_catid (SRTM_maxcat + 2) = 200000000 i1 = 0 ! count # of 30-arcsec pixels - + status = nf90_inq_varid(ncid, name='CatchIndex', varid=varid) do j=1,nr clat = -90. + float(j-1)*dy + dy/2. - - status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' + status = NF90_GET_VAR (ncid, varid, geos_msk, (/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/)) ! Read 10-arcsec rows that lie within the raster row 'j' + !status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' if(status /=0) then - PRINT *, NF_STRERROR(STATUS) + PRINT *, NF90_STRERROR(STATUS) print *, 'Problem with NF_GET_VARA_INT',trim(MaskFile),status endif @@ -304,14 +284,14 @@ PROGRAM mkSMAPTilesPara_v2 if(catid_index (i,j) == SRTM_maxcat + 2) veg (i,j) = IceType if((catid_index(i,j) >= 1).and.(catid_index (i,j) <= SRTM_maxcat)) i1 = i1 + 1 - ! count in if this is i,j pixel is a land, lake or ice within ind_col,ind_row SMAP grid cell + ! count in if this is i,j pixel is a land, lake or ice within ind_col,ind_row EASE grid cell - call easeV2_convert(trim(MGRID), clat, clon, r_smap, s_smap) + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) - ind_col = nint(r_smap) + 1 - ind_row = nint(s_smap) + 1 + ind_col = nint(r_ease) + 1 + ind_row = nint(s_ease) + 1 - if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_smap)) then + if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_ease)) then l= ind_row*ND + ind_col if(veg(i,j)==LakeType) then @@ -326,7 +306,7 @@ PROGRAM mkSMAPTilesPara_v2 end do enddo - status = NF_CLOSE (ncid) + status = NF90_CLOSE (ncid) deallocate (geos_msk) print *,'Read ', trim (MaskFile) @@ -335,7 +315,7 @@ PROGRAM mkSMAPTilesPara_v2 else ! Old IGBP (Veg) + HYDRO1k (catchments) based mask will - ! Overlaid on SMAP mask + ! Overlaid on EASE mask ! ----------------------------------------------------- allocate(iaster (i_raster,j_raster)) @@ -446,11 +426,11 @@ PROGRAM mkSMAPTilesPara_v2 print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) ! While looping through the nc x nr grid (tile raster), this section counts # of - ! SMAP grid cells that contain land, ice or water, seperately. - ! Each SMAP grid cell is assigned with an ID = ind_row*ND + ind_col. + ! EASE grid cells that contain land, ice or water, seperately. + ! Each EASE grid cell is assigned with an ID = ind_row*ND + ind_col. ! This is just the prelimiminery assessment in the process of assigning separate - ! tiles for land, water and ice fractions within the SMAP Grid cell - ! The program checks each nc x nr pixels whether there is a SMAP grid cell underneath, and counts + ! tiles for land, water and ice fractions within the EASE Grid cell + ! The program checks each nc x nr pixels whether there is a EASE grid cell underneath, and counts ! number of water, land and ice pixels as seen on veg raster. ! ----------------------------------------------------------------------------------------------- @@ -462,12 +442,12 @@ PROGRAM mkSMAPTilesPara_v2 do j =nr ,1 ,-1 clat = -90. + float(j-1)*dy + dy/2. - call easeV2_convert(trim(MGRID), clat, clon, r_smap, s_smap) + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) - ind_col = nint(r_smap) + 1 - ind_row = nint(s_smap) + 1 + ind_col = nint(r_ease) + 1 + ind_row = nint(s_ease) + 1 - if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_smap)) then + if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_ease)) then l= ind_row*ND + ind_col if(veg(i,j)==LakeType) then @@ -503,9 +483,9 @@ PROGRAM mkSMAPTilesPara_v2 deallocate (raster) - print *,'# of Land pixels in SMAP: ',sum (land_id) - print *,'# of water pixels in SMAP: ',sum (water_id) - print *,'# of ice pixels in SMAP: ',sum (ice_id) + print *,'# of Land pixels in EASE: ',sum (land_id) + print *,'# of water pixels in EASE: ',sum (water_id) + print *,'# of ice pixels in EASE: ',sum (ice_id) l_index=0 w_index=sum (land_id) @@ -514,7 +494,7 @@ PROGRAM mkSMAPTilesPara_v2 allocate(tile_area (1:i_index + sum (ice_id))) - allocate(smap_grid_area (1:NT)) + allocate(ease_grid_area (1:NT)) allocate(tile_ele (1:w_index)) allocate(tile_area_land(1:w_index)) allocate(my_land (1:i_index + sum (ice_id))) @@ -526,18 +506,18 @@ PROGRAM mkSMAPTilesPara_v2 my_land = 0 all_id = 0 - smap_grid_area = 0. + ease_grid_area = 0. tile_area_land = 0. tile_ele = 0. tile_area = 0. ! While looping through the nc x nr grid, this section derives land, ice and water tiles. - ! Each SMAP grid cell is assigned with an ID = ind_row*ND + ind_col - ! ind_col, ind_row are overlying SMAP grid cell indices + ! Each EASE grid cell is assigned with an ID = ind_row*ND + ind_col + ! ind_col, ind_row are overlying EASE grid cell indices ! Based on the above sums: ! l_index Grid cells have land fractions (sum(land_id)) - ! w_index SMAP Grid cells have inland water fractions (sum(water_id)) - ! i_index SMAP Grid cells have ice fractions (sum(ice_id)) + ! w_index EASE Grid cells have inland water fractions (sum(water_id)) + ! i_index EASE Grid cells have ice fractions (sum(ice_id)) ! hence, tile_index 1 to l_index represent land tiles ! tile_index l_index +1 to l_index + w_index represent water (lakes) tiles ! tile_index l_index + w_index +1 to l_index + w_index + i_index represent ice tiles @@ -554,15 +534,15 @@ PROGRAM mkSMAPTilesPara_v2 do j =nr ,1 ,-1 lats = -90._8 + (j - 0.5_8)*dy clat = -90. + float(j-1)*dy + dy/2. - call easeV2_convert(trim(MGRID), clat, clon, r_smap, s_smap) + call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) - ind_col = nint(r_smap) + 1 - ind_row = nint(s_smap) + 1 + ind_col = nint(r_ease) + 1 + ind_row = nint(s_ease) + 1 l= ind_row*ND + ind_col pix_area =(sin(d2r*(lats+0.5*dy)) -sin(d2r*(lats-0.5*dy)))*(dx*d2r) - if((ind_row.ge.1).and.(veg(i,j).ge.1).and.(ind_row.le.nr_smap)) then + if((ind_row.ge.1).and.(veg(i,j).ge.1).and.(ind_row.le.nr_ease)) then if(veg(i,j)==LakeType) then if(water_id(l)==0) then @@ -600,8 +580,8 @@ PROGRAM mkSMAPTilesPara_v2 all_id (tileid_index(i,j)) = j*ND_raster + i endif - if((ind_row.ge.1).and.(ind_row.le.nr_smap)) then - smap_grid_area(l) = smap_grid_area(l) + & + if((ind_row.ge.1).and.(ind_row.le.nr_ease)) then + ease_grid_area(l) = ease_grid_area(l) + & pix_area endif @@ -645,7 +625,7 @@ PROGRAM mkSMAPTilesPara_v2 print *,'Total Land Area :', sum(tile_area(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000., & sum(tile_area_land(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. - print *,'Creating ... ', trim(gfile)//'rst' + print *,'Creating ... ', trim(gfile)//'.rst' !------------------------------------------- @@ -658,7 +638,7 @@ PROGRAM mkSMAPTilesPara_v2 close (10,status='keep') - print *,'Creating ... ', trim(gfile)//'til ,catchment.def' + print *,'Creating ... ', trim(gfile)//'.til ,catchment.def' !----------------------------------------------------------- @@ -669,9 +649,9 @@ PROGRAM mkSMAPTilesPara_v2 open (10, file ='til/'//trim(gfile)//'.til',form='formatted',status='unknown',action='write') write (10,*)i_index,SRTM_maxcat, nc, nr write (10,*)1 - write (10,*)'SMAP-EASEv2-'//trim(MGRID) - write (10,*)nc_smap - write (10,*)nr_smap + write (10,*)EASELabel + write (10,*)nc_ease + write (10,*)nr_ease ! write (10,*)'NO-OCEAN' ! write (10,*) -9999 ! write (10,*) -9999 @@ -694,20 +674,20 @@ PROGRAM mkSMAPTilesPara_v2 if (l <= l_index) then typ = 100 - call easeV2_inverse (trim(MGRID), real(ig-1),real(jg-1), clat, clon) + call EASE_inverse (EASELabel, real(ig-1),real(jg-1), clat, clon) - mnx = clon - 180./real(nc_smap) - mxx = clon + 180./real(nc_smap) + mnx = clon - 180./real(nc_ease) + mxx = clon + 180./real(nc_ease) jgv = real(jg-1) + 0.5 - call easeV2_inverse (trim(MGRID), real(ig-1),jgv, clat, clon) + call EASE_inverse (EASELabel, real(ig-1),jgv, clat, clon) mny = clat jgv = real(jg-1) - 0.5 - call easeV2_inverse (trim(MGRID), real(ig-1),jgv, clat, clon) + call EASE_inverse (EASELabel, real(ig-1),jgv, clat, clon) mxy = clat @@ -715,9 +695,9 @@ PROGRAM mkSMAPTilesPara_v2 endif - call easeV2_inverse (trim(MGRID), real(ig-1), real(jg-1), clat, clon) + call EASE_inverse (EASELabel, real(ig-1), real(jg-1), clat, clon) - fr_gcm= tile_area(l)/smap_grid_area(jg*ND + ig) + fr_gcm= tile_area(l)/ease_grid_area(jg*ND + ig) if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then write(10,'(i10,i9,2f10.4,2i6,f19.12,i10,i15,e13.4)') & @@ -733,7 +713,7 @@ PROGRAM mkSMAPTilesPara_v2 close (11,status='keep') deallocate (tileid_index,catid_index,veg) - deallocate (tile_area, smap_grid_area, tile_ele, tile_area_land, my_land, all_id) + deallocate (tile_area, ease_grid_area, tile_ele, tile_area_land, my_land, all_id) if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then @@ -741,7 +721,7 @@ PROGRAM mkSMAPTilesPara_v2 !--------------------------------------------------- - deallocate (SRTM_CatchArea, SRTM_catid) + deallocate (SRTM_CatchArea, SRTM_catid, SRTM_catid_r8) endif @@ -753,89 +733,97 @@ PROGRAM mkSMAPTilesPara_v2 ! now run mkCatchParam ! -------------------- - tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBCSV) - write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr - tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1) - print *,trim(tmpstring) - - call execute_command_line (tmpstring) - - contains - - ! ------------------------------------------------------------------------------- - - SUBROUTINE mkEASEv2Raster - - implicit none - - integer :: i, j, i_ease, j_ease - real*8, allocatable :: xs(:,:), ys(:,:) - real :: x,y, xout, yout - - allocate (xs ( nc_smap+1, nr_smap+1)) - allocate (ys ( nc_smap+1, nr_smap+1)) - - do j = 1, nr_smap+1 - do i = 1, nc_smap+1 - x = real(i-1) -0.5 - y = real(nr_smap - j)+0.5 - call easeV2_inverse(MGRID, x, y, yout, xout) - ys (i,j) = dble(yout) - xs (i,j) = dble(xout) - end do - end do - - call LRRasterize(EASElabel,xs,ys,nc=nc,nr=nr,xmn = xs(1,1), xmx= xs(nc_smap+1, nr_smap+1), ymn=ys(1,1), ymx = ys(nc_smap+1, nr_smap+1), Here=.false., Verb=.false.) - - stop - end SUBROUTINE mkEASEv2Raster + ! WY Note: now mkCatchParam is run in the make_bcs script, not here + ! and nthread will be reset to run mkCatchParam - ! ------------------------------------------------------------ + ! tmpstring1 = '-e EASE -g '//trim(gfile)//' -v '//trim(LBCSV) + ! write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr + ! tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1) + ! print *,trim(tmpstring) - SUBROUTINE write_tilfile - - implicit none - - character*200 :: infile - integer :: NT, NF, NC, NR, NPF, NG, IDUM, i, N, icol, rcol - character*20 :: cdum - integer, dimension (:,:), allocatable :: iRtable - real, dimension (:,:), allocatable :: rRtable - - infile = 'til/'//trim(EASElabel)//'_'//trim(EASElabel)//'-Pfafstetter.' - - open (10,file = trim(infile)//'ind', form = 'formatted', action = 'read', status = 'old') - open (11,file = trim(infile)//'TIL', form = 'formatted', action = 'write', status = 'unknown') - - read (10, *) NT, NF, NC, NR - write (11,'(4I10)')NT, NF, NC, NR - read (10, *) NG - write(11, *) NG - - do n = 1, NG - read (10, '(a)') cdum - write(11, '(a)') trim (cdum) - read (10, *) IDUM - write(11, '(I10)') IDUM - read (10, *) IDUM - write(11, '(I10)') IDUM - end do - - icol = 7 - rcol = 5 - allocate (iRtable (1, 1:icol)) - allocate (rRtable (1, 1:rcol)) - - do n = 1, nt - read(10,'(I10,3E20.12,9(2I10,E20.12,I10))') iRtable (1,1),rRtable(1,1),rRtable(1,2),rRtable(1,3),iRtable (1,2),iRtable (1,3),rRtable(1,4),iRtable (1,4),& - iRtable (1,5),iRtable (1,6),rRtable(1,5),iRtable (1,7) - write(11,'(I10,3E20.12,9(2I10,E20.12,I10))') iRtable (1,1),rRtable(1,1),rRtable(1,2),rRtable(1,3),iRtable (1,2)-1,nr_smap - iRtable (1,3),rRtable(1,4),iRtable (1,4),& - iRtable (1,5),iRtable (1,6),rRtable(1,5),iRtable (1,7) - end do - - close (10, status = 'keep') - close (11, status = 'keep') - - END SUBROUTINE write_tilfile + ! call execute_command_line (tmpstring) + + +!!! commented out. It may be used in the future for irrigation tiles +!!! contains +!!! +!!! ! ------------------------------------------------------------------------------- +!!! +!!! SUBROUTINE mkEASEv2Raster +!!! +!!! implicit none +!!! +!!! integer :: i, j, i_ease, j_ease +!!! real*8, allocatable :: xs(:,:), ys(:,:) +!!! real :: x,y, xout, yout +!!! +!!! allocate (xs ( nc_ease+1, nr_ease+1)) +!!! allocate (ys ( nc_ease+1, nr_ease+1)) +!!! +!!! do j = 1, nr_ease+1 +!!! do i = 1, nc_ease+1 +!!! x = real(i-1) -0.5 +!!! y = real(nr_ease - j)+0.5 +!!! call EASE_inverse(MGRID, x, y, yout, xout) +!!! ys (i,j) = dble(yout) +!!! xs (i,j) = dble(xout) +!!! end do +!!! end do +!!! +!!! call LRRasterize(EASElabel,xs,ys,nc=nc,nr=nr,xmn = xs(1,1), xmx= xs(nc_ease+1, nr_ease+1), & +!!! ymn=ys(1,1), ymx = ys(nc_ease+1, nr_ease+1), Here=.false., Verb=.false.) +!!! +!!! stop +!!! end SUBROUTINE mkEASEv2Raster +!!! +!!! ! ------------------------------------------------------------ +!!! +!!! SUBROUTINE write_tilfile +!!! +!!! implicit none +!!! +!!! character*200 :: infile +!!! integer :: NT, NF, NC, NR, NPF, NG, IDUM, i, N, icol, rcol +!!! character*20 :: cdum +!!! integer, dimension (:,:), allocatable :: iRtable +!!! real, dimension (:,:), allocatable :: rRtable +!!! +!!! infile = 'til/'//trim(EASElabel)//'_'//trim(EASElabel)//'-Pfafstetter.' +!!! +!!! open (10,file = trim(infile)//'ind', form = 'formatted', action = 'read', status = 'old') +!!! open (11,file = trim(infile)//'TIL', form = 'formatted', action = 'write', status = 'unknown') +!!! +!!! read (10, *) NT, NF, NC, NR +!!! write (11,'(4I10)')NT, NF, NC, NR +!!! read (10, *) NG +!!! write(11, *) NG +!!! +!!! do n = 1, NG +!!! read (10, '(a)') cdum +!!! write(11, '(a)') trim (cdum) +!!! read (10, *) IDUM +!!! write(11, '(I10)') IDUM +!!! read (10, *) IDUM +!!! write(11, '(I10)') IDUM +!!! end do +!!! +!!! icol = 7 +!!! rcol = 5 +!!! allocate (iRtable (1, 1:icol)) +!!! allocate (rRtable (1, 1:rcol)) +!!! +!!! do n = 1, nt +!!! read(10,'(I10,3E20.12,9(2I10,E20.12,I10))') iRtable (1,1),rRtable(1,1), & +!!! rRtable(1,2),rRtable(1,3),iRtable (1,2),iRtable (1,3),rRtable(1,4),iRtable (1,4),& +!!! iRtable (1,5),iRtable (1,6),rRtable(1,5),iRtable (1,7) +!!! write(11,'(I10,3E20.12,9(2I10,E20.12,I10))') iRtable (1,1),rRtable(1,1), & +!!! rRtable(1,2),rRtable(1,3),iRtable (1,2)-1,nr_ease - iRtable (1,3),rRtable(1,4),iRtable (1,4),& +!!! iRtable (1,5),iRtable (1,6),rRtable(1,5),iRtable (1,7) +!!! end do +!!! +!!! close (10, status = 'keep') +!!! close (11, status = 'keep') +!!! +!!! END SUBROUTINE write_tilfile END PROGRAM diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90 deleted file mode 100644 index fc62e7718..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mkSMAPTilesPara.F90 +++ /dev/null @@ -1,696 +0,0 @@ -PROGRAM mkSMAPTilesPara_v1 -! This program constructs land and lake tiles for the SMAP-EASE-M09 and M36 grids (just set MGRID) -! for CLSM implementation. -! f90 -c create_smap_tiles.f90 -! f90 -c smapconv.f -! f90 -o create_smap_tiles create_smap_tiles.o smapconv.o -! - use easeV1_conv - use rmTinyCatchParaMod - use process_hres_data - use MAPL_SortMod - use MAPL_ConstantsMod - - implicit none - - integer i,j,ig,jg,i0,iop,n,d1,d2,j1,j2,i1,i2,ix, jx,icount,pcount - integer :: NC = i_raster, NR = j_raster, NT = 16330000, ND = 10000, ND_raster = 10000 - - integer, parameter :: nc_esa = 129600, nr_esa = 64800 - - ! For regridding - - integer, allocatable, target, dimension (:,:) & - :: geos_msk - REAL, allocatable, DIMENSION (:) :: loc_val - INTEGER, ALLOCATABLE, DIMENSION (:) :: density, loc_int - logical, dimension (:), allocatable :: unq_mask - integer, pointer , dimension (:,:) :: subset - integer, pointer , dimension (:) :: subset1, subset_smap - real, pointer , dimension (:) :: subset2 - integer :: dx_esa, dy_esa, NBINS, NPLUS - - integer*8, allocatable, dimension (:) :: SRTM_catid - - integer,allocatable, dimension (:,:), target :: tileid_index,catid_index - integer,allocatable, dimension (:,:) :: catid, iaster - integer,allocatable, dimension (:) :: land_id,water_id,ice_id - integer,allocatable, dimension (:) :: my_land, all_id - real, allocatable, dimension (:) :: smap_grid_area,tile_area,SRTM_CatchArea - integer*1,allocatable, dimension (:,:) :: veg, i2aster - real*4, dimension (:,:), allocatable :: q0,raster - REAL, dimension (:), allocatable :: tile_ele, tile_area_land - - INTEGER*8 :: PFAF_CODE - integer l,imn,imx,jmn,jmx,mval,l_index,i_index,w_index,typ,pfaf,cindex - integer :: LakeType, IceType, OceanType - character(3) :: easegrid - real :: clat, clon, r_smap, s_smap, smap_convert, da - real :: fr_gcm - integer :: ind_col, ind_row, status, ncid, nciv,nland_cells, DOM_INDX - REAL (kind=8), PARAMETER :: RADIUS=6378137.0,pi=3.14159265358979323846 - character*100 :: veg_class (12) - character*5 :: MGRID - character*100 :: gfile,gtopo30 - integer :: nc_smap,nr_smap, N_args, command_argument_count - real :: EASE_grid_area, CELL_km - REAL :: dx,dy,d2r,lats,mnx,mxx,mny,mxy,sum1,sum2,jgv, VDUM,pix_area - character(40) :: arg - character*200 :: tmpstring, tmpstring1, tmpstring2 - logical :: regrid = .false. - character*128 :: MaskFile - include 'netcdf.inc' - - N_args = command_argument_count() - - if(N_args < 1) then - print *,'USAGE : bin/mkSMAPTiles_v1 -smap_grid MXX' - print *,'Allowed SMAP grids are: M25' - stop - end if - - i=0 - - do while ( i < N_args ) - - i = i+1 - - call get_command_argument(i,arg) - - if ( trim(arg) == '-smap_grid' ) then - i = i+1 - call get_command_argument(i,MGRID) - - else ! stop for any other arguments - - print *,'USAGE : bin/mkSMAPTiles -smap_grid MXX' - print *,'Allowed SMAP grids are: M09 M36 Ml' - stop - - endif - - end do - - call execute_command_line('cd data/ ; ln -s /discover/nobackup/projects/gmao/ssd/land/l_data/LandBCs_files_for_mkCatchParam/V001/ CATCH') - call execute_command_line('cd ..') - - - ! Setting SMAP Grid specifications - ! -------------------------------- - - if (trim(MGRID) == 'M25') then - CELL_km = 25.067525 ! nominal cell size in kilometers - nc_smap = 1383 - nr_smap = 586 - gfile = 'SMAP_EASE_M25_1383x586' - EASE_grid_area = CELL_km*CELL_km - - else ! - - print *,'Unknown SMAP Grid stopping..' - stop - - endif - - allocate(land_id (1:NT)) - allocate(water_id (1:NT)) - allocate(ice_id (1:NT)) - land_id = 0 - water_id = 0 - ice_id = 0 - OceanType = 0 - IceType =11 - LakeType =10 - - ND = 10*10**(nint(log10(1.*nr_smap))) - - ! Check for the 10 arc-sec MaskFile - ! ----------------------------------- - - call get_environment_variable ("MASKFILE" ,MaskFile ) - - print *, 'Using MaskFile ', trim(MaskFile) - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - ! New ESA (Veg) + SRTM (catchments) based mask file - ! is overlaid on SMAP - ! ------------------------------------------------- - - nc = 43200 ! Number of rows in raster file - nr = 21600 ! Number of columns in raster file - - regrid = .true. - dx_esa = nc_esa / nc ! x-dimension (or # of ESA columns within the raster grid cell) - dy_esa = nr_esa / nr ! y-dimension (or # of ESA rows within the raster grid cell) - - allocate(tileid_index(1:nc,1:nr)) - allocate(SRTM_catid (1:SRTM_maxcat+2)) - allocate(catid_index (1:nc,1:nr)) - allocate(veg (1:nc,1:nr)) - allocate(geos_msk (1:nc_esa,1:dy_esa)) - allocate(SRTM_CatchArea (1:SRTM_maxcat)) - - OPEN (10, FILE = 'data/CATCH/SRTM-TopoData/Pfafcatch-routing.dat', & - FORM = 'FORMATTED',STATUS='OLD',ACTION='READ') - - READ (10,*) I - DO N = 1, I - READ (10, '(i8,i15,4(1x,f9.4),1x,e10.3,4(1x,e9.3),I8,6(1x,f9.4))') & - DOM_INDX,PFAF_CODE,VDUM,VDUM,VDUM,VDUM,VDUM, & - SRTM_CatchArea (N) - END DO - CLOSE (10, STATUS='KEEP') - - dx = 360._8/nc - dy = 180._8/nr - d2r = PI/180._8 - da = MAPL_radius*MAPL_radius*pi*pi*dx*dy/180./180./1000000. - - tileid_index = 0 - catid_index = 0 - veg = 0 - - status = NF_OPEN ('data/CATCH/GEOS5_10arcsec_mask.nc', NF_NOWRITE, ncid) - status = NF_GET_VARA_INT64 (ncid,3,(/1/),(/SRTM_maxcat/),SRTM_catid(1:SRTM_maxcat)) ! Read pfafstetter IDs - if(status /=0) then - PRINT *, NF_STRERROR(STATUS) - print *, 'Problem with NF_OPEN',trim(MaskFile) - endif - - SRTM_catid (SRTM_maxcat + 1) = 190000000 - SRTM_catid (SRTM_maxcat + 2) = 200000000 - i1 = 0 ! count # of 30-arcsec pixels - - do j=1,nr - - clat = -90. + float(j-1)*dy + dy/2. - - status = NF_GET_VARA_INT (ncid,4,(/1,(j-1)*dy_esa +1/),(/nc_esa,dy_esa/),geos_msk) ! Read 10-arcsec rows that lie within the raster row 'j' - - if(status /=0) then - PRINT *, NF_STRERROR(STATUS) - print *, 'Problem with NF_GET_VARA_INT',trim(MaskFile),status - endif - - do i = 1,nc - - clon = -180. + float(i-1)*dx + dx/2. - - if (associated (subset)) NULLIFY (subset) - subset => geos_msk ((i-1)*dx_esa + 1 : i*dx_esa, 1:dy_esa) ! rectangular array contains ESA pixels that lie within the raster grid cell at i,j - if(maxval (subset) > SRTM_maxcat) then - where (subset == 190000000) subset = SRTM_maxcat + 1 - where (subset == 200000000) subset = SRTM_maxcat + 2 - endif - - if (maxval (subset) > 0) then ! check whether there are Non-ocean ESA pixels - ! catID of the middle pixel - - veg (i,j) = 1 ! veg is set to land - - NPLUS = count(subset >= 1 .and. subset <= SRTM_maxcat + 2) ! Count non-ocean ESA pixels within - allocate (loc_int (1:NPLUS)) - allocate (unq_mask(1:NPLUS)) - loc_int = pack(subset,mask = (subset >= 1 .and. subset <= SRTM_maxcat + 2)) ! loc_int contains catch_indices of non-ocean ESA pixels - call MAPL_Sort (loc_int) - unq_mask = .true. - do n = 2,NPLUS - unq_mask(n) = .not.(loc_int(n) == loc_int(n-1)) ! count number of unique numbers in loc_int for binning - end do - NBINS = count(unq_mask) - - if (NBINS > 1) then - allocate(loc_val (1:NBINS)) - allocate(density (1:NBINS)) - loc_val = 1.*pack(loc_int,mask =unq_mask) ! loc_val contains available non-ocean catch_indices within the i,j grid cell, - ! Those numbers will be used as bin values - call histogram (dx_esa*dy_esa, NBINS, density, loc_val, real(subset)) ! density is the pixel count for each bin value - catid_index (i,j) = loc_val (maxloc(density,1)) ! picks maximum density as the dominant catchment_index at i,j - deallocate (loc_val, density) - else - catid_index (i,j) = loc_int (1) - endif - deallocate (loc_int, unq_mask) - - if(catid_index (i,j) == SRTM_maxcat + 1) veg (i,j) = LakeType - if(catid_index (i,j) == SRTM_maxcat + 2) veg (i,j) = IceType - if((catid_index(i,j) >= 1).and.(catid_index (i,j) <= SRTM_maxcat)) i1 = i1 + 1 - - ! count in if this is i,j pixel is a land, lake or ice within ind_col,ind_row SMAP grid cell - - call easeV1_convert(trim(MGRID), clat, clon, r_smap, s_smap) - - ind_col = nint(r_smap) + 1 - ind_row = nint(s_smap) + 1 - - if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_smap)) then - l= ind_row*ND + ind_col - - if(veg(i,j)==LakeType) then - water_id(l) = 1 - else if(veg(i,j)==IceType) then - ice_id (l) = 1 - else - land_id (l) = 1 - endif - endif - endif - end do - enddo - - status = NF_CLOSE (ncid) - deallocate (geos_msk) - - print *,'Read ', trim (MaskFile) - print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) - - else - - ! Old IGBP (Veg) + HYDRO1k (catchments) based mask will - ! Overlaid on SMAP mask - ! ----------------------------------------------------- - - allocate(iaster (i_raster,j_raster)) - allocate(i2aster (i_raster,j_raster)) - allocate(veg (1:nc,1:nr)) - allocate(catid (1:nc,1:nr)) - allocate(catid_index (1:nc,1:nr)) - allocate(tileid_index(1:nc,1:nr)) - - dx = 360._8/nc - dy = 180._8/nr - d2r = PI/180._8 - da = MAPL_radius*MAPL_radius*pi*pi*dx*dy/180./180./1000000. - - tileid_index = 0 - - ! Simple Biosphere 2 Model Legend - ! Value Class Name - ! (ftp://edcftp.cr.usgs.gov/pub/data/glcc/globe/latlon/sib22_0.leg) - ! the types vary 0-11 (array index minus 1) - - veg_class(1) = 'Ocean' - veg_class(2) = 'Broadleaf Evergreen Trees' - veg_class(3) = 'Broadleaf Deciduous Trees' - veg_class(4) = 'Broadleaf and Needleleaf Trees' - veg_class(5) = 'Needleleaf Evergreen Trees' - veg_class(6) = 'Needleleaf Deciduous Trees' - veg_class(7) = 'Short Vegetation/C4 Grassland' - veg_class(8) = 'Shrubs with Bare Soil' - veg_class(9) = 'Dwarf Trees and Shrubs' - veg_class(10) = 'Agriculture or C3 Grassland' - veg_class(11) = 'Water, Wetlands' - veg_class(12) = 'Ice/Snow' - - ! reading SiB2 land cover classification data - the origin of the - ! 2.5'x2.5' vegetation raster file is global 1min IGBP data - ! (ftp://edcftp.cr.usgs.gov/pub/data/glcc/globe/latlon/sib22_0.leg) - - open (10,file='data/CATCH/sib22.5_v2.0.dat', & - form='unformatted', & - action='read', convert='big_endian',status='old') - - READ(10)i2aster - - close (10,status='keep') - - if(regrid) then - call RegridRaster1 (i2aster,veg) - else - veg = i2aster - endif - - deallocate (i2aster) - - ! reading 2.5'x2.5' global raster file of Pfafstetter Catchment IDs - ! In this version, the dateline has been overlaid over the catchments those straddle - ! across. The numbers contain for - ! 1 global ocean catchment : Pfafstetter ID 0 - ! 36716 global land catchments : Pfafstetter IDs 1000-5999900 - ! 1 global inland water (lakes) catchment : Pfafstetter ID 6190000 - ! 1 global ice catchment : Pfafstetter ID 6200000 - - open (10,file='data/CATCH/global.cat_id.catch.DL', form='formatted', & - action='read', status='old')! - - do j=1,j_raster - read(10,*)(iaster(i,j),i=1,i_raster) - end do - - close (10,status='keep') - - if(regrid) then - call RegridRaster(iaster,catid) - else - catid = iaster - endif - - print *,'Read global.cat_id.catch.DL' - print *,'Min and Max of Pfafstetter IDs:', minval(catid),maxval(catid) - - ! reading the 2.5'x2.5' global raster file of tile indices for the - ! above Pfafstetter Catchments - ! 1 global ocean catchment : tile_index 36719 - ! 36716 global land catchments : tile_index 1-36716 - ! 1 global inland water (lakes) catchment : tile_index 36717 - ! 1 global ice catchment : tile_index 36718 - ! ------------------------------------------------------------ - - open (10,file='data/CATCH/' & - //'PfafstatterDL.rst', form='unformatted', & - action='read',convert='little_endian', status='old') - - do j=1,j_raster - read(10)(iaster(i,j),i=1,i_raster) - end do - - close (10,status='keep') - - if(regrid) then - call RegridRaster(iaster,catid_index) - else - catid_index = iaster - endif - - deallocate (iaster) - - print *,'Read PfafstatterDL.rst' - print *,'Min and Max of tile indices:',minval(catid_index),maxval(catid_index) - - ! While looping through the nc x nr grid (tile raster), this section counts # of - ! SMAP grid cells that contain land, ice or water, seperately. - ! Each SMAP grid cell is assigned with an ID = ind_row*ND + ind_col. - ! This is just the prelimiminery assessment in the process of assigning separate - ! tiles for land, water and ice fractions within the SMAP Grid cell - ! The program checks each nc x nr pixels whether there is a SMAP grid cell underneath, and counts - ! number of water, land and ice pixels as seen on veg raster. - ! ----------------------------------------------------------------------------------------------- - - - do i = 1 ,nc - - clon = -180. + float(i-1)*dx + dx/2. - - do j =nr ,1 ,-1 - - clat = -90. + float(j-1)*dy + dy/2. - call easeV1_convert(trim(MGRID), clat, clon, r_smap, s_smap) - - ind_col = nint(r_smap) + 1 - ind_row = nint(s_smap) + 1 - - if((ind_row.ge.1).and.(veg(i,j).ne.OceanType).and.(ind_row.le.nr_smap)) then - l= ind_row*ND + ind_col - - if(veg(i,j)==LakeType) then - water_id(l) = 1 - else if(veg(i,j)==IceType) then - ice_id (l) = 1 - else - land_id (l) = 1 - endif - endif - end do - end do - - endif - - ! Reading SRTM elevation data - to be consistent with AGCM - ! -------------------------------------------------------- - - allocate(raster (i_raster,j_raster)) - allocate(q0(nc,nr)) - - gtopo30 = 'data/CATCH/srtm30_withKMS_2.5x2.5min.data' - - open (10,file=trim(gtopo30),form='unformatted',status='old',convert='little_endian') - read (10) raster - close (10,status='keep') - - if(regrid) then - call RegridRasterReal(raster,q0) - else - q0 = raster - endif - - deallocate (raster) - - print *,'# of Land pixels in SMAP: ',sum (land_id) - print *,'# of water pixels in SMAP: ',sum (water_id) - print *,'# of ice pixels in SMAP: ',sum (ice_id) - - l_index=0 - w_index=sum (land_id) - i_index=sum (land_id) + sum (water_id) - nland_cells = w_index - - - allocate(tile_area (1:i_index + sum (ice_id))) - allocate(smap_grid_area (1:NT)) - allocate(tile_ele (1:w_index)) - allocate(tile_area_land(1:w_index)) - allocate(my_land (1:i_index + sum (ice_id))) - allocate(all_id (1:i_index + sum (ice_id))) - - land_id = 0 - water_id= 0 - ice_id = 0 - - my_land = 0 - all_id = 0 - smap_grid_area = 0. - tile_area_land = 0. - tile_ele = 0. - tile_area = 0. - - ! While looping through the nc x nr grid, this section derives land, ice and water tiles. - ! Each SMAP grid cell is assigned with an ID = ind_row*ND + ind_col - ! ind_col, ind_row are overlying SMAP grid cell indices - ! Based on the above sums: - ! l_index Grid cells have land fractions (sum(land_id)) - ! w_index SMAP Grid cells have inland water fractions (sum(water_id)) - ! i_index SMAP Grid cells have ice fractions (sum(ice_id)) - ! hence, tile_index 1 to l_index represent land tiles - ! tile_index l_index +1 to l_index + w_index represent water (lakes) tiles - ! tile_index l_index + w_index +1 to l_index + w_index + i_index represent ice tiles - ! global nc x nr array of tileid_index(nc,nr) contains corresponding tile_index values which - ! is derived in the below loop - - ND_raster = 10*10**(nint(log10(1.*NR))) - i2 = 1 - - do i = 1 ,nc - - clon = -180. + float(i-1)*dx + dx/2. - - do j =nr ,1 ,-1 - lats = -90._8 + (j - 0.5_8)*dy - clat = -90. + float(j-1)*dy + dy/2. - call easeV1_convert(trim(MGRID), clat, clon, r_smap, s_smap) - - ind_col = nint(r_smap) + 1 - ind_row = nint(s_smap) + 1 - - l= ind_row*ND + ind_col - pix_area =(sin(d2r*(lats+0.5*dy)) -sin(d2r*(lats-0.5*dy)))*(dx*d2r) - - if((ind_row.ge.1).and.(veg(i,j).ge.1).and.(ind_row.le.nr_smap)) then - - if(veg(i,j)==LakeType) then - if(water_id(l)==0) then - w_index = w_index + 1 - water_id(l) = w_index - tileid_index(i,j)= water_id(l) - else - tileid_index(i,j)= water_id(l) - endif - endif - - if(veg(i,j)==IceType) then - if(ice_id(l)==0) then - i_index = i_index + 1 - ice_id (l) = i_index - tileid_index(i,j)= ice_id (l) !i_index - else - tileid_index(i,j)= ice_id (l) !i_index - endif - endif - - if(veg(i,j).lt.LakeType) then - if(land_id(l)==0) then - l_index = l_index + 1 - land_id (l) = l_index - tileid_index(i,j)= land_id (l) !1-l_index - else - tileid_index(i,j)= land_id (l) !1-l_index - endif - - endif - tile_area(tileid_index(i,j))= tile_area(tileid_index(i,j)) + & - pix_area - my_land(tileid_index(i,j)) = l - all_id (tileid_index(i,j)) = j*ND_raster + i - endif - - if((ind_row.ge.1).and.(ind_row.le.nr_smap)) then - smap_grid_area(l) = smap_grid_area(l) + & - pix_area - endif - - ! computing tile area/elevation - ! ----------------------------- - - if((tileid_index(i,j) > 0).and.(tileid_index(i,j) <= nland_cells))then - tile_ele(tileid_index(i,j)) = tile_ele(tileid_index(i,j)) + q0(i,j) * & - pix_area - tile_area_land(tileid_index(i,j)) = tile_area_land(tileid_index(i,j)) + & - pix_area - endif - end do - end do - - deallocate(land_id, q0) - deallocate(water_id) - deallocate(ice_id ) - - tile_ele = tile_ele/tile_area_land - - ! adjustment Global Mean Topography to 614.649 (615.662 GTOPO 30) m - ! ----------------------------------------------------------------- - sum1=0. - sum2=0. - - do j=1,l_index - sum1 = sum1 + tile_ele(j)*tile_area(j) - enddo - - if(sum1/sum(tile_area(1:l_index)).ne. 614.649D0 ) then - print *,'Global Mean Elevation (over land): ', sum1/sum(tile_area(1:l_index)) - tile_ele =tile_ele*(614.649D0 / (sum1/sum(tile_area(1:l_index)))) - sum1=0. - sum2=0. - do j=1,l_index - sum1 = sum1 + tile_ele(j)*tile_area(j) - enddo - print *,'Global Mean Elevation after scaling to SRTM : ',sum1/sum(tile_area(1:l_index)) - endif - print *,'Total Land Area :', sum(tile_area(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000., & - sum(tile_area_land(1:l_index))* MAPL_RADIUS * MAPL_RADIUS/1000./1000. - - print *,'Creating ... ', trim(gfile)//'rst' - - !------------------------------------------- - - open (10, file ='rst/'//trim(gfile)//'.rst',form='unformatted',status='unknown', & - action='write') - - do j=1,nr - write(10)(tileid_index(i,j),i=1,nc) - end do - - close (10,status='keep') - - print *,'Creating ... ', trim(gfile)//'til ,catchment.def' - - !----------------------------------------------------------- - - open (11,file='clsm/catchment.def', & - form='formatted',status='unknown') - write(11,*)l_index - - open (10, file ='til/'//trim(gfile)//'.til',form='formatted',status='unknown',action='write') - write (10,*)i_index, SRTM_maxcat, nc, nr - write (10,*)1 - write (10,*)'SMAP-EASEv2-'//trim(MGRID) - write (10,*)nc_smap - write (10,*)nr_smap -! write (10,*)'NO-OCEAN' -! write (10,*) -9999 -! write (10,*) -9999 - - do l=1,i_index - - ig = my_land(l)-ND*(my_land(l)/ND) - jg = my_land(l)/ND - - cindex= catid_index(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - pfaf = cindex - else - pfaf = catid(all_id(l)-ND_raster*(all_id(l)/ND_raster),all_id(l)/ND_raster) - endif - - if ((l > l_index).and.(l <= w_index)) typ =19 - if (l > w_index) typ = 20 - - if (l <= l_index) then - typ = 100 - call easeV1_inverse (trim(MGRID), real(ig-1),real(jg-1), clat, clon) - - mnx = clon - 180./real(nc_smap) - mxx = clon + 180./real(nc_smap) - - jgv = real(jg-1) + 0.5 - - call easeV1_inverse (trim(MGRID), real(ig-1),jgv, clat, clon) - - mny = clat - - jgv = real(jg-1) - 0.5 - - call easeV1_inverse (trim(MGRID), real(ig-1),jgv, clat, clon) - - mxy = clat - - write (11,'(i8,i8,5(2x,f9.4), i4)')l,pfaf,mnx,mxx,mny,mxy,tile_ele(l) - - endif - - call easeV1_inverse (trim(MGRID), real(ig-1), real(jg-1), clat, clon) - - fr_gcm= tile_area(l)/smap_grid_area(jg*ND + ig) - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - write(10,'(i10,i9,2f10.4,2i5,f19.12,i10,i15,e13.4)') & - typ,pfaf,clon,clat,ig-1,jg-1,fr_gcm ,pfaf,SRTM_catid(cindex) - else - write(10,'(i10,i9,2f10.4,2i5,f19.12,i10,e13.4,i8)') & - typ,pfaf,clon,clat,ig-1,jg-1,fr_gcm ,cindex - endif - - end do - - close (10,status='keep') - close (11,status='keep') - - deallocate (tileid_index,catid_index,veg) - deallocate (tile_area, smap_grid_area, tile_ele, tile_area_land, my_land, all_id) - - if (index(MaskFile,'GEOS5_10arcsec_mask') /= 0) then - - print *,'Creating SMAP-Catch_TransferData.nc files.' - - !--------------------------------------------------- - - deallocate (SRTM_CatchArea, SRTM_catid) - - endif - - ! create Grid2Catch transfer file - ! ------------------------------- - - ! CALL CREATE_ROUT_PARA_FILE (NC, NR, trim(gfile), MGRID=MGRID) - - ! now run mkCatchParam - ! -------------------- - - tmpstring1 = '-e EASE -g '//trim(gfile) - write(tmpstring2,'(2(a2,x,i5,x))')'-x',nc,'-y',nr - tmpstring = 'bin/mkCatchParam.x '//trim(tmpstring2)//' '//trim(tmpstring1) - print *,trim(tmpstring) - - call execute_command_line (tmpstring) - - END PROGRAM mkSMAPTilesPara_v1 - diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 index f7d47a432..8e4a16005 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/rmTinyCatchParaMod.F90 @@ -177,6 +177,11 @@ SUBROUTINE init_bcs_config (LBCSV) use_PEATMAP = .true. jpl_height = .false. + case default + + print *,'init_bcs_config(): unknown land boundary conditions version (LBCSV)' + stop + end select END SUBROUTINE init_bcs_config diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 index 1d71c3f12..4e6b2d94f 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/mk_restarts/mk_CatchRestarts.F90 @@ -571,7 +571,6 @@ SUBROUTINE read_and_write_rst (NTILES, SURFLAY, OutIsOld, NTILES_IN, idi, rc) endif HAVE if (filetype == 0) then - call MAPL_VarWrite(OutFmt,names(1),BF1(Idx)) call MAPL_VarWrite(OutFmt,names(2),BF2(Idx)) call MAPL_VarWrite(OutFmt,names(3),BF3(Idx))