From 7a979448706f03273cd417bc2bd284b05934a8df Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Tue, 28 May 2024 15:38:29 -0400 Subject: [PATCH 1/7] buble up error from MAPL_getHorzIJindex function --- CHANGELOG.md | 1 + base/Base/Base_Base_implementation.F90 | 7 +++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 33c9f9cf8461..c34f0a892028 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Changed +- Propgated the error message from MAPL_Horzijindex function - fixed a bug in generate_newnxy in MAPL_SwathGridFactory.F90 (NX*NY=Ncore) - pFIO Clients don't send "Done" message when there is no request - Update `components.yaml` diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index 20b27e36deea..05eaadc22063 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -2629,7 +2629,7 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc) call ESMF_AttributeGet(grid, name='GridType', value=grid_type, _RC) if(trim(grid_type) == "Cubed-Sphere") then - call MAPL_GetGlobalHorzIJIndex(npts, II, JJ, lon=lon, lat=lat, lonR8=lonR8, latR8=latR8, Grid=Grid, rc=rc) + call MAPL_GetGlobalHorzIJIndex(npts, II, JJ, lon=lon, lat=lat, lonR8=lonR8, latR8=latR8, Grid=Grid, _RC) call MAPL_Grid_Interior(Grid,i1,i2,j1,j2) ! convert index to local, if it is not in domain, set it to -1 just as the legacy code @@ -2770,7 +2770,10 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, ! make sure the grid can be used in this subroutine good_grid = grid_is_ok(grid) - _ASSERT( good_grid, "MAPL_GetGlobalHorzIJIndex cannot handle this grid") + + if ( .not. good_grid ) then + _FAIL( "MAPL_GetGlobalHorzIJIndex cannot handle this grid") + endif allocate(lons(npts),lats(npts)) if (present(lon) .and. present(lat)) then From 2d84cea0663dcca03d2f752e327915b4a607a946 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 30 Sep 2024 09:49:17 -0400 Subject: [PATCH 2/7] rename loacl target_lons and target_lats to avoid confusion of the stretched grid --- base/Base/Base_Base_implementation.F90 | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index 2eb23d4303c0..c5780649a185 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -2594,8 +2594,8 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc) real(ESMF_KIND_R8), allocatable :: elats(:) integer :: i,iiloc,jjloc, i1, i2, j1, j2 real(ESMF_KIND_R4) :: lonloc,latloc - logical :: localSearch - real(ESMF_KIND_R8), allocatable :: target_lons(:),target_lats(:) + logical :: localSearch + real(ESMF_KIND_R8), allocatable :: tmp_lons(:),tmp_lats(:) type(ESMF_CoordSys_Flag) :: coordSys character(len=ESMF_MAXSTR) :: grid_type @@ -2616,13 +2616,13 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc) localSearch = .false. end if - allocate(target_lons(npts),target_lats(npts)) + allocate(tmp_lons(npts),tmp_lats(npts)) if (present(lon) .and. present(lat)) then - target_lons = lon - target_lats = lat + tmp_lons = lon + tmp_lats = lat else if (present(lonR8) .and. present(latR8)) then - target_lons = lonR8 - target_lats = latR8 + tmp_lons = lonR8 + tmp_lats = latR8 end if !AOO change tusing GridType atribute if (im_world*6==jm_world) then @@ -2661,8 +2661,8 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc) ! lat-lon grid goes from -180 to 180 shift if we must ! BMA this -180 to 180 might change at some point do i=1,npts - lonloc = target_lons(i) - latloc = target_lats(i) + lonloc = tmp_lons(i) + latloc = tmp_lats(i) if (lonloc > MAPL_PI) lonloc = lonloc - 2.0*MAPL_PI IIloc = ijsearch(elons,im+1,lonloc,.false.) JJloc = ijsearch(elats,jm+1,latloc,.false.) @@ -2672,7 +2672,7 @@ module subroutine MAPL_GetHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, rc) deallocate(elons,elats) end if - deallocate(target_lons, target_lats) + deallocate(tmp_lons, tmp_lats) _RETURN(ESMF_SUCCESS) contains From 922f616539a9762756a97b6b33bdb07209c8f481 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 30 Sep 2024 14:38:30 -0400 Subject: [PATCH 3/7] Added MAPL_Reverse_Schmidt to reverse the stretched grid --- CHANGELOG.md | 2 + base/Base/Base_Base.F90 | 16 ++++ base/Base/Base_Base_implementation.F90 | 121 ++++++++++++++++++++++--- 3 files changed, 124 insertions(+), 15 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 793e069197e2..2eadb4a79e67 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ### Added +- Added MAPL_Reverse_Schmidt to reverse the stretched grid for indices computation + ### Changed - Propagated the error message from MAPL_HorzIJIndex subroutine diff --git a/base/Base/Base_Base.F90 b/base/Base/Base_Base.F90 index dcef19267c47..e5baa73393c5 100644 --- a/base/Base/Base_Base.F90 +++ b/base/Base/Base_Base.F90 @@ -61,6 +61,7 @@ module MAPL_Base public MAPL_FieldBundleDestroy public MAPL_GetHorzIJIndex public MAPL_GetGlobalHorzIJIndex + public MAPL_Reverse_Schmidt public MAPL_GenGridName public MAPL_GenXYOffset public MAPL_GeosNameNew @@ -712,6 +713,21 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, integer, optional, intent(out ) :: rc ! return code end subroutine MAPL_GetGlobalHorzIJIndex + module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts,lon,lat,lonR8,latR8, lonRe, latRe, rc) + use ESMF, only: ESMF_KIND_R8, ESMF_GRid + implicit none + !ARGUMENTS: + type(ESMF_Grid), intent(inout) :: Grid ! ESMF grid + logical, intent(out ) :: stretched + integer, intent(in ) :: npts ! number of points in lat and lon arrays + real, optional, intent(in ) :: lon(npts) ! array of longitudes in radians + real, optional, intent(in ) :: lat(npts) ! array of latitudes in radians + real(ESMF_KIND_R8), optional, intent(in ) :: lonR8(npts) ! array of longitudes in radians + real(ESMF_KIND_R8), optional, intent(in ) :: latR8(npts) ! array of latitudes in radians + real(ESMF_KIND_R8), optional, intent(out ) :: lonRe(npts) ! array of longitudes in radians + real(ESMF_KIND_R8), optional, intent(out ) :: latRe(npts) ! array of latitudes in radians + integer, optional, intent(out ) :: rc ! return code + end subroutine MAPL_Reverse_Schmidt module subroutine MAPL_GenGridName(im, jm, lon, lat, xyoffset, gridname, geos_style) integer :: im, jm diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index c5780649a185..4fecd89782e8 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -24,6 +24,8 @@ ! !USES: ! use ESMF + use ESMFL_Mod + use MAPL_FieldUtils use MAPL_Constants use MAPL_RangeMod @@ -2596,6 +2598,7 @@ 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 @@ -2756,7 +2759,7 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, ! MAPL_PI_R8/18, Japan Fuji mountain shift real(ESMF_KIND_R8), parameter :: shift= 0.174532925199433d0 - logical :: good_grid + logical :: good_grid, stretched ! Return if no local points _RETURN_IF(npts == 0) @@ -2769,6 +2772,10 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, JM_World = dims(2) _ASSERT( IM_WORLD*6 == JM_WORLD, "It only works for cubed-sphere grid") + allocate(lons(npts),lats(npts)) + + call MAPL_Reverse_Schmidt(Grid, stretched, npts, lon=lon, lat=lat, lonR8=lonR8, latR8=latR8, lonRe=lons, latRe=lats, _RC) + dalpha = 2.0d0*alpha/IM_WORLD ! make sure the grid can be used in this subroutine @@ -2778,15 +2785,6 @@ module subroutine MAPL_GetGlobalHorzIJIndex(npts,II,JJ,lon,lat,lonR8,latR8,Grid, _FAIL( "MAPL_GetGlobalHorzIJIndex cannot handle this grid") endif - allocate(lons(npts),lats(npts)) - if (present(lon) .and. present(lat)) then - lons = lon - lats = lat - else if (present(lonR8) .and. present(latR8)) then - lons = lonR8 - lats = latR8 - end if - ! shift the grid away from Japan Fuji Mt. lons = lons + shift @@ -2870,7 +2868,8 @@ function grid_is_ok(grid) result(OK) logical :: OK integer :: I1, I2, J1, J2, j real(ESMF_KIND_R8), pointer :: corner_lons(:,:), corner_lats(:,:) - real(ESMF_KIND_R8) :: accurate_lat, accurate_lon + real(ESMF_KIND_R8), allocatable :: lonRe(:), latRe(:) + real(ESMF_KIND_R8) :: accurate_lat, accurate_lon, stretch_factor, target_lon, target_lat real :: tolerance tolerance = epsilon(1.0) @@ -2882,12 +2881,15 @@ 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 - corner_lons(1,1)) > tolerance) then + if (abs(accurate_lon - lonRe(1)) > tolerance) then print*, "accurate_lon: ", accurate_lon - print*, "corner_lon : ", corner_lons(1,1) + print*, "corner_lon : ", lonRe(1) print*, "Error: Grid should have pi/18 shift" OK = .false. return @@ -2896,9 +2898,9 @@ function grid_is_ok(grid) result(OK) do j = J1+1, J2 accurate_lat = -alpha + (j-1)*dalpha - if ( abs(accurate_lat - corner_lats(1,j-J1+1)) > 5.0*tolerance) then + if ( abs(accurate_lat - latRe(j-J1+1)) > 5.0*tolerance) then print*, "accurate_lat: ", accurate_lat - print*, "edge_lat : ", corner_lats(1,j-J1+1) + print*, "edge_lat : ", latRe(j-J1+1) print*, "edge point : ", j print*, "Error: It could be " print*, " 1)Grid is NOT gnomonic_ed;" @@ -3383,4 +3385,93 @@ module function MAPL_GetCorrectedPhase(gc,rc) result(phase) _RETURN(_SUCCESS) end function MAPL_GetCorrectedPhase + module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts, lon, lat, lonR8, latR8, lonRe, latRe, rc) + type(ESMF_Grid), intent(inout) :: Grid + logical, intent(out) :: stretched + integer, intent(in ) :: npts ! number of points in lat and lon arrays + real, optional, intent(in ) :: lon(npts) ! array of longitudes in radians + real, optional, intent(in ) :: lat(npts) ! array of latitudes in radians + real(ESMF_KIND_R8), optional, intent(in ) :: lonR8(npts) ! array of longitudes in radians + real(ESMF_KIND_R8), optional, intent(in ) :: latR8(npts) ! + real(ESMF_KIND_R8), optional, intent(out ) :: lonRe(npts) ! + real(ESMF_KIND_R8), optional, intent(out ) :: latRe(npts) ! + integer, optional, intent(out) :: rc + + logical :: factorPresent, lonPresent, latPresent + integer :: status + real(ESMF_KIND_R8) :: c2p1, c2m1, half_pi, two_pi, stretch_factor, target_lon, target_lat + real(ESMF_KIND_R8), dimension(npts) :: x,y,z, Xx, Yy, Zz + logical, dimension(npts) :: n_s + + _RETURN_IF( npts == 0 ) + + call ESMF_AttributeGet(grid, name='STRETCH_FACTOR', isPresent= factorPresent, _RC) + call ESMF_AttributeGet(grid, name='TARGET_LON', isPresent= lonPresent, _RC) + call ESMF_AttributeGet(grid, name='TARGET_LAT', isPresent= latPresent, _RC) + + if ( factorPresent .and. lonPresent .and. latPresent) then + stretched = .true. + else + stretched = .false. + endif + + if (present(lonRe) .and. present(latRe)) then + if (present(lonR8) .and. present(latR8)) then + lonRe = lonR8 + latRe = latR8 + else if (present(lon) .and. present(lat)) then + lonRe = lon + latRe = lat + else + _FAIL("Need input to get the output lonRe, latRe") + endif + else + _RETURN(_SUCCESS) + endif + + if (.not. stretched) then + _RETURN(_SUCCESS) + endif + + call ESMF_AttributeGet(grid, name='STRETCH_FACTOR', value=stretch_factor, _RC) + call ESMF_AttributeGet(grid, name='TARGET_LON', value=target_lon, _RC) + call ESMF_AttributeGet(grid, name='TARGET_LAT', value=target_lat, _RC) + + c2p1 = 1 + stretch_factor*stretch_factor + c2m1 = 1 - stretch_factor*stretch_factor + + half_pi = MAPL_PI/2 + two_pi = MAPL_PI*2 + + x = cos(latRe)*cos(lonRe - target_lon) + y = cos(latRe)*sin(lonRe - target_lon) + z = sin(latRe) + + Xx = sin(target_lat)*x - cos(target_lat)*z + Yy = -y + Zz = -cos(target_lat)*x - sin(target_lat)*z + + n_s = (1. - abs(Zz)) < 10**(-7) + + where(n_s) + lonRe = 0.0d0 + latRe = half_pi*sign(1.0, Zz) + elsewhere + lonRe = atan2(Yy,Xx) + latRe = asin(Zz) + endwhere + + if (abs(c2m1) > 10**(-7)) then !# unstretch + latRe = asin( (c2m1-c2p1*sin(latRe))/(c2m1*sin(latRe)-c2p1)) + endif + + where ( lonRe < 0) + lonRe = lonRe + two_pi + elsewhere (lonRe >= two_pi) + lonRe = lonRe - two_pi + endwhere + + _RETURN(_SUCCESS) + end subroutine + end submodule Base_Implementation From d3a0810121474b19a80d4c748302e4f35f511e6b Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Mon, 30 Sep 2024 16:22:03 -0400 Subject: [PATCH 4/7] clean up error info --- base/Base/Base_Base_implementation.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index 4fecd89782e8..ca2c1bc8d822 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -2882,7 +2882,8 @@ function grid_is_ok(grid) result(OK) 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) + 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 @@ -2905,7 +2906,6 @@ function grid_is_ok(grid) result(OK) 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)" - print*, " 3)This is a stretched grid which is not yet supported" OK = .false. return endif From 9462f11fab338bd7e27c90671ee97d29975b3f7d Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang <52509753+weiyuan-jiang@users.noreply.github.com> Date: Mon, 30 Sep 2024 16:38:15 -0400 Subject: [PATCH 5/7] Update Base_Base_implementation.F90 keep the types of "sign" arguments consistent --- base/Base/Base_Base_implementation.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index ca2c1bc8d822..66c1b551f1ff 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -3455,7 +3455,7 @@ module subroutine MAPL_Reverse_Schmidt(Grid, stretched, npts, lon, lat, lonR8, l where(n_s) lonRe = 0.0d0 - latRe = half_pi*sign(1.0, Zz) + latRe = half_pi*sign(1.0d0, Zz) elsewhere lonRe = atan2(Yy,Xx) latRe = asin(Zz) From 7637bd1ff1b3d6aad4f88d997a9973cbabab09e8 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Thu, 3 Oct 2024 10:54:16 -0400 Subject: [PATCH 6/7] fixed bug on units --- base/Base/Base_Base_implementation.F90 | 68 ++++++++++++++------------ 1 file changed, 37 insertions(+), 31 deletions(-) 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) From 7db8a1ef2e3d5dfe3b97144a2dbd3e983f027fa3 Mon Sep 17 00:00:00 2001 From: Weiyuan Jiang Date: Fri, 4 Oct 2024 10:17:47 -0400 Subject: [PATCH 7/7] added more error information --- base/Base/Base_Base_implementation.F90 | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/base/Base/Base_Base_implementation.F90 b/base/Base/Base_Base_implementation.F90 index 70b33640f724..6b4ff3cbce7a 100644 --- a/base/Base/Base_Base_implementation.F90 +++ b/base/Base/Base_Base_implementation.F90 @@ -2896,20 +2896,12 @@ function grid_is_ok(grid) result(OK) 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 - print*, "Error: Grid should have pi/18 shift" - OK = .false. - return - endif - - if ( any(abs(accurate_lat - latRe) > 2.0*tolerance)) then - print*, "accurate_lat: ", accurate_lat - print*, "edge_lat : ", latRe + if (any(abs(accurate_lon - lonRe) > 2.0* tolerance) .or. any(abs(accurate_lat - latRe) > 2.0*tolerance)) then 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)" + print*, " 1) grid may not have pi/18 Japan mountain shift" + print*, " 2) grid is NOT gnomonic_ed;" + print*, " 3) lats lons from MAPL_GridGetCorners are NOT accurate (single precision from ESMF)" + print*, " 4) strtech grid rotates north pole" OK = .false. return endif