diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index 66c1b551f1ff..70b33640f724 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -134,7 +134,6 @@ module subroutine MAPL_FieldAllocCommit(field, dims, location, typekind, & integer :: ub1, ub2, ub3 ! SSI - character(len=ESMF_MAXSTR) :: name type(ESMF_Pin_Flag) :: pinflag type(ESMF_VM) :: vm logical :: ssiSharedMemoryEnabled @@ -2598,7 +2597,6 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc) real(ESMF_KIND_R4) :: lonloc,latloc logical :: localSearch real(ESMF_KIND_R8), allocatable :: tmp_lons(:),tmp_lats(:) - real(ESMF_KIND_R8) :: lonRe(npts), latRe(npts) ! returned from reversed_shmidt transform type(ESMF_CoordSys_Flag) :: coordSys character(len=ESMF_MAXSTR) :: grid_type @@ -2747,7 +2745,7 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, real(ESMF_KIND_R8), allocatable, dimension (:,:) :: xyz real(ESMF_KIND_R8), allocatable, dimension (:) :: x,y,z real(ESMF_KIND_R8), allocatable :: max_abs(:) - real(ESMF_KIND_R8) :: dalpha + real(ESMF_KIND_R8) :: dalpha, shift0 real(ESMF_KIND_R8), allocatable :: lons(:), lats(:) ! sqrt(2.0d0), distance from center to the mid of an edge for a 2x2x2 cube @@ -2786,7 +2784,9 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, endif ! shift the grid away from Japan Fuji Mt. - lons = lons + shift + shift0 = shift + if (stretched) shift0 = 0 + lons = lons + shift0 ! get xyz from sphere surface allocate(xyz(3, npts), max_abs(npts)) @@ -2869,7 +2869,8 @@ function grid_is_ok(grid) result(OK) integer :: I1, I2, J1, J2, j real(ESMF_KIND_R8), pointer :: corner_lons(:,:), corner_lats(:,:) real(ESMF_KIND_R8), allocatable :: lonRe(:), latRe(:) - real(ESMF_KIND_R8) :: accurate_lat, accurate_lon, stretch_factor, target_lon, target_lat + real(ESMF_KIND_R8), allocatable :: accurate_lat(:), accurate_lon(:) + real(ESMF_KIND_R8) :: stretch_factor, target_lon, target_lat, shift0 real :: tolerance tolerance = epsilon(1.0) @@ -2881,36 +2882,38 @@ function grid_is_ok(grid) result(OK) call ESMF_GridGetCoord(grid,localDE=0,coordDim=2,staggerloc=ESMF_STAGGERLOC_CORNER, & farrayPtr=corner_lats, rc=status) - allocate(lonRe(j2-j1+2), latRe(j2-j1+2)) - call MAPL_Reverse_Schmidt(grid, stretched, J2-J1+1, lonR8=corner_lons(1,:), & - latR8=corner_lats(1,:), lonRe=lonRe, latRe=latRe, _RC) - if ( I1 ==1 .and. J2<=IM_WORLD ) then - if (J1 == 1) then - accurate_lon = 1.750d0*MAPL_PI_R8 - shift - if (abs(accurate_lon - lonRe(1)) > tolerance) then + if ( I1 == 1 .and. J1 == 1 ) then + allocate(lonRe(j2-j1+1), latRe(j2-j1+1)) + call MAPL_Reverse_Schmidt(grid, stretched, J2-J1+1, lonR8=corner_lons(1,:), & + latR8=corner_lats(1,:), lonRe=lonRe, latRe=latRe, _RC) + + allocate(accurate_lon(j2-j1+1), accurate_lat(j2-j1+1)) + + shift0 = shift + if (stretched) shift0 = 0 + + accurate_lon = 1.750d0*MAPL_PI_R8 - shift0 + accurate_lat = [(-alpha + (j-1)*dalpha, j = j1, j2)] + + if (any(abs(accurate_lon - lonRe) > 2.0* tolerance)) then print*, "accurate_lon: ", accurate_lon - print*, "corner_lon : ", lonRe(1) + print*, "corner_lon : ", lonRe print*, "Error: Grid should have pi/18 shift" OK = .false. return - endif endif - do j = J1+1, J2 - accurate_lat = -alpha + (j-1)*dalpha - if ( abs(accurate_lat - latRe(j-J1+1)) > 5.0*tolerance) then - print*, "accurate_lat: ", accurate_lat - print*, "edge_lat : ", latRe(j-J1+1) - print*, "edge point : ", j - print*, "Error: It could be " - print*, " 1)Grid is NOT gnomonic_ed;" - print*, " 2)lats lons from MAPL_GridGetCorners are NOT accurate (single precision from ESMF)" - OK = .false. - return - endif - enddo - endif + if ( any(abs(accurate_lat - latRe) > 2.0*tolerance)) then + print*, "accurate_lat: ", accurate_lat + print*, "edge_lat : ", latRe + print*, "Error: It could be " + print*, " 1)Grid is NOT gnomonic_ed;" + print*, " 2)lats lons from MAPL_GridGetCorners are NOT accurate (single precision from ESMF)" + OK = .false. + return + endif + endif end function end subroutine MAPL_GetGlobalHorzIJIndex @@ -3207,7 +3210,7 @@ module subroutine MAPL_FieldSplit(field, fields, aliasName, rc) ! local vars integer :: status integer :: k, n - integer :: k1,k2,kk + integer :: k1,k2 logical :: has_ungrd integer :: ungrd_cnt integer :: fieldRank @@ -3440,8 +3443,11 @@ module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts, lon, lat, lonR8, l c2p1 = 1 + stretch_factor*stretch_factor c2m1 = 1 - stretch_factor*stretch_factor - half_pi = MAPL_PI/2 - two_pi = MAPL_PI*2 + half_pi = MAPL_PI_R8/2 + two_pi = MAPL_PI_R8*2 + + target_lon = target_lon*MAPL_DEGREES_TO_RADIANS_R8 + target_lat = target_lat*MAPL_DEGREES_TO_RADIANS_R8 x = cos(latRe)*cos(lonRe - target_lon) y = cos(latRe)*sin(lonRe - target_lon)