diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt index c966e4df1..c52656386 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/CMakeLists.txt @@ -2,7 +2,6 @@ esma_set_this() set (srcs LDAS_DateTimeMod.F90 -EASE_conv.F90 mod_process_hres_data.F90 rasterize.F90 read_riveroutlet.F90 diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_conv.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_conv.F90 deleted file mode 100644 index 2060d77d5..000000000 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/EASE_conv.F90 +++ /dev/null @@ -1,786 +0,0 @@ - -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/makebcs/TileFile_ASCII_to_nc4.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/TileFile_ASCII_to_nc4.F90 index 8f489537e..2e04f5afd 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/TileFile_ASCII_to_nc4.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/TileFile_ASCII_to_nc4.F90 @@ -7,9 +7,8 @@ program TileFile_ASCII_to_nc4 use, intrinsic :: iso_fortran_env, only: REAL64 - use MAPL - use LogRectRasterizeMod, only: MAPL_UNDEF_R8 - use EASE_conv, only: ease_extent + use MAPL, only: MAPL_WriteTilingNC4, MAPL_ease_extent + use LogRectRasterizeMod, only: MAPL_UNDEF_R8 implicit none @@ -114,7 +113,7 @@ program TileFile_ASCII_to_nc4 ! rTable(:,4) is tile area fraction within grid cell (fr), convert to area; ! get fr back in WriteTilingNC4 - call ease_extent(gName1, tmp_in1, tmp_in2, cell_area=cell_area) ! get EASE grid cell area + call MAPL_ease_extent(gName1, tmp_in1, tmp_in2, cell_area=cell_area) ! get EASE grid cell area rTable(:,3) = rTable(:,4)*cell_area rTable(:,4) = cell_area diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py index 725cc34a6..386d834c4 100755 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs_questionary.py @@ -75,7 +75,7 @@ def get_configs_from_answers(answers): if lbcsv in ['F25', 'GM4', 'ICA']: maskfile = 'global.cat_id.catch.GreatLakesCaspian_Updated.DL' if (maskfile == ''): - print(" \!\!\!\! Invalid Ocean Resolution, stopping ") + print(" !!! Invalid Ocean Resolution, stopping !!!") exit() if 'EASEv1' == grid_type or 'EASEv2' == grid_type: diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 index 942d85ff1..fcafbf333 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkCatchParam.F90 @@ -1,3 +1,6 @@ +#define I_AM_MAIN +#include "MAPL_ErrLog.h" + PROGRAM mkCatchParam ! !INTERFACE: @@ -18,10 +21,10 @@ PROGRAM mkCatchParam ! ! Sarith Mahanama - March 23, 2012 ! Email: sarith.p.mahanama@nasa.gov - use MAPL - use EASE_conv + use MAPL, only: MAPL_ease_extent, MAPL_ReadTilingNC4 use rmTinyCatchParaMod use process_hres_data + use MAPL_ExceptionHandling ! use module_irrig_params, ONLY : create_irrig_params @@ -190,7 +193,7 @@ PROGRAM mkCatchParam if (index(Gridname,'EASEv') /=0) then ! here Gridname has alias EASELabel - call ease_extent(Gridname, nc_ease, nr_ease ) + call MAPL_ease_extent(Gridname, nc_ease, nr_ease, _RC) write(nc_string, '(i0)') nc_ease write(nr_string, '(i0)') nr_ease Gridname = trim(Gridname)//'_'//trim(nc_string)//'x'//trim(nr_string) diff --git a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 index 7cb59e4bb..5856811b8 100644 --- a/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 +++ b/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mkEASETilesParam.F90 @@ -29,17 +29,16 @@ PROGRAM mkEASETilesParam ! - removed repetition of identical operations ! - added comments ! - white-space changes for improved readability - - use EASE_conv, only : EASE_extent, EASE_convert, EASE_inverse + use, intrinsic :: iso_fortran_env, only: REAL64 use rmTinyCatchParaMod, only : i_raster, j_raster use rmTinyCatchParaMod, only : RegridRasterReal use rmTinyCatchParaMod, only : Target_mean_land_elev use process_hres_data, only : histogram use LogRectRasterizeMod, only : SRTM_maxcat, MAPL_UNDEF_R8 ! rasterize.F90 use MAPL_SortMod - use MAPL_ConstantsMod + use MAPL_ConstantsMod, only : MAPL_PI_r8, MAPL_RADIUS use MAPL_ExceptionHandling - use MAPL + use MAPL, only : MAPL_ease_extent, MAPL_ease_convert, MAPL_ease_inverse, MAPL_WriteTilingNC4 use netcdf implicit none @@ -158,7 +157,7 @@ PROGRAM mkEASETilesParam EASElabel = trim(EASELabel_) - call ease_extent( EASELabel, nc_ease, nr_ease ) + call MAPL_ease_extent( EASELabel, nc_ease, nr_ease, _RC) write(nc_string, '(i0)') nc_ease write(nr_string, '(i0)') nr_ease @@ -373,7 +372,7 @@ PROGRAM mkEASETilesParam ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) - call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) + call MAPL_EASE_convert(EASELabel, clat, clon, r_ease, s_ease, _RC) ind_col = nint(r_ease) + 1 ind_row = nint(s_ease) + 1 ! can be negative or greater than nr_ease (lat near N/S pole) @@ -543,7 +542,7 @@ PROGRAM mkEASETilesParam ! do j =nr ,1 ,-1 ! ! clat = -90. + float(j-1)*dy + dy/2. -! call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) +! call MAPL_EASE_convert(EASELabel, clat, clon, r_ease, s_ease) ! ! ind_col = nint(r_ease) + 1 ! ind_row = nint(s_ease) + 1 @@ -657,7 +656,7 @@ PROGRAM mkEASETilesParam ! get 1-based ind_col and ind_row indices of EASE grid cell that contains raster grid cell (i,j) - call EASE_convert(EASELabel, clat, clon, r_ease, s_ease) + call MAPL_EASE_convert(EASELabel, clat, clon, r_ease, s_ease, _RC) ind_col = nint(r_ease) + 1 ind_row = nint(s_ease) + 1 ! NOTE: can be zero or negative or greater than nr_ease (lat near N/S pole) @@ -884,20 +883,20 @@ PROGRAM mkEASETilesParam ! get min/max lat/lon of EASE grid cell ! BUG: This is *not* the desired min/max lat/lon of the land tile!!! - call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) + call MAPL_ease_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon, _RC) mnx = clon - dx_ease mxx = clon + dx_ease jgv = real(jg-1) + 0.5 - call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) + call MAPL_ease_inverse( EASELabel, real(ig-1), jgv, clat, clon, _RC) mny = clat jgv = real(jg-1) - 0.5 - call EASE_inverse( EASELabel, real(ig-1), jgv, clat, clon ) + call MAPL_ease_inverse( EASELabel, real(ig-1), jgv, clat, clon, _RC) mxy = clat @@ -921,7 +920,7 @@ PROGRAM mkEASETilesParam ! contributing raster grid cells, which is *not* the same for all EASE grid cells; ! that is, cannot use exact (globally constant) area of EASE grid cell. - call EASE_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon ) + call MAPL_ease_inverse( EASELabel, real(ig-1), real(jg-1), clat, clon, _RC) fr_gcm = tile_area(nn) / ease_grid_area((jg-1)*nc_ease+ig) @@ -1032,7 +1031,7 @@ PROGRAM mkEASETilesParam !!! 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) +!!! call MAPL_ease_inverse(MGRID, x, y, yout, xout) !!! ys (i,j) = dble(yout) !!! xs (i,j) = dble(xout) !!! end do