57 USE constants_mod,
ONLY: deg_to_rad, rad_to_deg, radius
67 #include<file_version.h>
73 REAL,
allocatable,
DIMENSION(:,:) :: glo_lat
74 REAL,
allocatable,
DIMENSION(:,:) :: glo_lon
75 REAL,
allocatable,
DIMENSION(:,:) :: aglo_lat
78 REAL,
allocatable,
DIMENSION(:,:) :: aglo_lon
89 INTEGER :: tile_number
93 CHARACTER(len=128) :: grid_type
131 TYPE(
domain2d),
INTENT(in) :: domain
132 CLASS(*),
INTENT(in),
DIMENSION(:,:) :: glo_lat
133 CLASS(*),
INTENT(in),
DIMENSION(:,:) :: glo_lon
134 CLASS(*),
INTENT(in),
DIMENSION(:,:) :: aglo_lat
135 CLASS(*),
INTENT(in),
DIMENSION(:,:) :: aglo_lon
137 INTEGER,
DIMENSION(1) :: tile
140 INTEGER :: i_dim, j_dim
141 INTEGER :: ai_dim, aj_dim
142 INTEGER,
DIMENSION(2) :: latdim, londim
143 INTEGER,
DIMENSION(2) :: alatdim, alondim
144 INTEGER :: mype, npes, npespertile
145 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: xbegin, xend, ybegin, yend
152 IF (
ALLOCATED(xbegin) )
DEALLOCATE(xbegin)
153 IF (
ALLOCATED(ybegin) )
DEALLOCATE(ybegin)
154 IF (
ALLOCATED(xend) )
DEALLOCATE(xend)
155 IF (
ALLOCATED(yend) )
DEALLOCATE(yend)
158 mype =
mpp_pe() -mpp_root_pe() + 1
162 ALLOCATE(xbegin(npes), &
165 & yend(npes), stat=stat)
166 IF ( stat .NE. 0 )
THEN
167 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
168 &
'Could not allocate memory for the compute grid indices&
177 npespertile = npes / ntiles
183 & xbegin=xbegin, xend=xend,&
184 & ybegin=ybegin, yend=yend)
190 latdim = shape(glo_lat)
191 londim = shape(glo_lon)
192 IF ( (latdim(1) == londim(1)) .AND.&
193 &(latdim(2) == londim(2)) )
THEN
197 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
198 &
'glo_lat and glo_lon must be the same shape.', fatal)
202 alatdim = shape(aglo_lat)
203 alondim = shape(aglo_lon)
204 IF ( (alatdim(1) == alondim(1)) .AND. &
205 &(alatdim(2) == alondim(2)) )
THEN
206 IF ( tile(1) == 4 .OR. tile(1) == 5 )
THEN
215 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
216 &
"a-grid's glo_lat and glo_lon must be the same shape.", fatal)
222 IF (
mpp_pe() == mpp_root_pe() ) &
223 &
CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
224 &
'The global grid has already been initialized', warning)
228 IF ( stat .NE. 0 )
THEN
229 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
230 &
'Could not allocate memory for the global grid.', fatal)
237 IF (
mpp_pe() == mpp_root_pe() ) &
238 &
CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
239 &
'The global a-grid has already been initialized', warning)
243 IF ( stat .NE. 0 )
THEN
244 CALL error_mesg(
'diag_global_mod::diag_grid_init',&
245 &
'Could not allocate memory for the global a-grid', fatal)
253 IF ( tile(1) == 4 .OR. tile(1) == 5 )
THEN
254 SELECT TYPE (aglo_lat)
255 TYPE IS (real(kind=r4_kind))
257 TYPE IS (real(kind=r8_kind))
260 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
261 &
'The a-grid latitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
264 SELECT TYPE (aglo_lon)
265 TYPE IS (real(kind=r4_kind))
267 TYPE IS (real(kind=r8_kind))
270 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
271 &
'The a-grid longitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
274 SELECT TYPE (aglo_lat)
275 TYPE IS (real(kind=r4_kind))
277 TYPE IS (real(kind=r8_kind))
280 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
281 &
'The a-grid latitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
284 SELECT TYPE (aglo_lon)
285 TYPE IS (real(kind=r4_kind))
287 TYPE IS (real(kind=r8_kind))
290 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
291 &
'The a-grid longitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
295 SELECT TYPE (glo_lat)
296 TYPE IS (real(kind=r4_kind))
298 TYPE IS (real(kind=r8_kind))
301 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
302 &
'The grid latitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
305 SELECT TYPE (glo_lon)
306 TYPE IS (real(kind=r4_kind))
308 TYPE IS (real(kind=r8_kind))
311 CALL error_mesg(
'diag_grid_mod::diag_grid_init',&
312 &
'The grid longitude data is not one of the supported types of real(kind=4) or real(kind=8)', fatal)
354 ELSE IF (
mpp_pe() == mpp_root_pe() )
THEN
355 CALL error_mesg(
'diag_grid_mod::diag_grid_end',&
356 &
'diag_global_grid%glo_lat was not allocated.', warning)
361 ELSE IF (
mpp_pe() == mpp_root_pe() )
THEN
362 CALL error_mesg(
'diag_grid_mod::diag_grid_end',&
363 &
'diag_global_grid%glo_lon was not allocated.', warning)
368 ELSE IF (
mpp_pe() == mpp_root_pe() )
THEN
369 CALL error_mesg(
'diag_grid_mod::diag_grid_end',&
370 &
'diag_global_grid%aglo_lat was not allocated.', warning)
375 ELSE IF (
mpp_pe() == mpp_root_pe() )
THEN
376 CALL error_mesg(
'diag_grid_mod::diag_grid_end',&
377 &
'diag_global_grid%aglo_lon was not allocated.', warning)
390 & istart, iend, jstart, jend)
391 REAL,
INTENT(in) :: latstart
392 REAL,
INTENT(in) :: lonstart
393 REAL,
INTENT(in) :: latend
394 REAL,
INTENT(in) :: lonend
395 INTEGER,
INTENT(out) :: istart
396 INTEGER,
INTENT(out) :: jstart
397 INTEGER,
INTENT(out) :: iend
398 INTEGER,
INTENT(out) :: jend
400 REAL,
ALLOCATABLE,
DIMENSION(:,:) :: delta_lat, delta_lon, grid_lon
402 REAL,
DIMENSION(4) :: dists_lon, dists_lat
403 REAL :: lonendadj, my_lonstart, my_lonend
405 INTEGER,
ALLOCATABLE,
DIMENSION(:,:) :: ijmin, ijmax
406 INTEGER :: mytile, ntiles, i, j, k, dimi, dimj, istat
414 REAL :: minimum_distance
415 REAL :: global_min_distance
419 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes',&
420 &
'Module not initialized, first initialze module with a call &
421 &to diag_grid_init', fatal)
424 if ( lonstart < 0. )
then
425 my_lonstart = lonstart + 360.
427 my_lonstart = lonstart
429 if ( lonend < 0. )
then
430 my_lonend = lonend + 360.
435 IF (latstart .EQ. latend .AND. my_lonstart .EQ. my_lonend)
THEN
443 allocate(ijmin(ntiles,2))
455 global_min_distance = minimum_distance
456 CALL mpp_min(global_min_distance)
460 IF (global_min_distance .EQ. minimum_distance)
THEN
468 IF (rank_buf .EQ. -1)
THEN
470 "No rank has minimum distance.", &
474 IF (rank_buf .EQ.
mpp_pe())
THEN
484 istart = ijmin(mytile,1)
485 jstart = ijmin(mytile,2)
496 ALLOCATE(ijmin(ntiles,2), stat=istat)
498 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes',&
499 &
'Cannot allocate ijMin index array', fatal)
500 ALLOCATE(ijmax(ntiles,2), stat=istat)
502 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes',&
503 &
'Cannot allocate ijMax index array', fatal)
514 ALLOCATE(delta_lat(dimi,dimj), stat=istat)
516 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes',&
517 &
'Cannot allocate latitude delta array', fatal)
518 ALLOCATE(delta_lon(dimi,dimj), stat=istat)
520 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes',&
521 &
'Cannot allocate longitude delta array', fatal)
550 IF ( dists_lon(k) > 180.0 )
THEN
551 dists_lon(k) = 360.0 - dists_lon(k)
554 delta_lon(i,j) = sum(dists_lon)/real(count)
555 delta_lat(i,j) = sum(dists_lat)/real(count)
563 ALLOCATE(grid_lon(dimi,dimj), stat=istat)
565 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes',&
566 &
'Cannot allocate temporary longitude array', fatal)
570 IF ( my_lonstart > my_lonend )
THEN
571 WHERE ( grid_lon < my_lonstart )
572 grid_lon = grid_lon + 360.0
574 lonendadj = my_lonend + 360.0
576 lonendadj = my_lonend
585 IF ( (abs(latstart)-delta_lat(i,j) <= 90.0 .AND.&
586 & 90.0 <= abs(latend)+delta_lat(i,j)) .AND.&
589 ELSE IF ( (my_lonstart-delta_lon(i,j) <= grid_lon(i,j) .AND.&
590 & grid_lon(i,j) < lonendadj+delta_lon(i,j)) )
THEN
604 DEALLOCATE(delta_lon)
605 DEALLOCATE(delta_lat)
616 IF ( ijmin(mytile,1) == huge(1) .OR. ijmax(mytile,1) == -huge(1) )
THEN
620 IF ( ijmin(mytile,2) == huge(1) .OR. ijmax(mytile,2) == -huge(1) )
THEN
625 istart = ijmin(mytile,1)
626 jstart = ijmin(mytile,2)
627 iend = ijmax(mytile,1)
628 jend = ijmax(mytile,2)
640 REAL,
INTENT(in) :: lat
641 REAL,
INTENT(in) :: lon
642 INTEGER,
INTENT(out) :: iindex
643 INTEGER,
INTENT(out) :: jindex
645 INTEGER :: indexes(2)
648 &
CALL error_mesg(
'diag_grid_mod::get_local_indexes2',&
649 &
'Module not initialized, first initialze module with a call &
650 &to diag_grid_init', fatal)
683 REAL,
INTENT(in) :: angle
696 REAL,
INTENT(in) :: angle
711 REAL,
INTENT(in) :: lat
712 REAL,
INTENT(in) :: lon
720 INTEGER :: nearestcorner
721 INTEGER,
DIMENSION(4,2) :: ijarray
724 REAL,
DIMENSION(4) :: pntdistances
725 TYPE(
point) :: origpt
726 REAL,
DIMENSION(4,2) :: cornerpts
727 TYPE(
point),
DIMENSION(9) :: points
728 REAL,
DIMENSION(9) :: distsqrd
738 IF ( lat == 90.0 )
THEN
740 ELSE IF ( lat == -90.0 )
THEN
749 iloop:
DO i=0, dimi-2
750 jloop:
DO j = 0, dimj-2
755 & (/ 4, 2 /), order=(/2,1/) )
757 maxctrdist = max(
gcirdistance(cornerpts(1,1),cornerpts(1,2), cornerpts(2,1),cornerpts(2,2)),&
758 &
gcirdistance(cornerpts(3,1),cornerpts(3,2), cornerpts(4,1),cornerpts(4,2)))
761 pntdistances =
gcirdistance(cornerpts(:,1),cornerpts(:,2), llat,llon)
763 IF ( (minval(pntdistances) <= maxctrdist) .AND. (i*j.NE.0) )
THEN
765 ijarray = reshape( (/ i, j, i+1, j+1, i+1, j, i, j+1 /), (/ 4, 2 /), order=(/2,1/) )
768 nearestcorner = minloc(pntdistances,1)
770 indxi = ijarray(nearestcorner,1)
771 indxj = ijarray(nearestcorner,2)
780 valid:
IF ( (indxi <= 0 .OR. dimi-1 <= indxi) .OR. &
781 & (indxj <= 0 .OR. dimj-1 <= indxj) )
THEN
827 SELECT CASE (minloc(distsqrd,1))
875 REAL,
INTENT(in) :: lat
876 REAL,
INTENT(in) :: lon
878 INTEGER :: indxi, indxj
880 INTEGER :: dimi, dimj
882 INTEGER :: jstart, jend, nextj
883 TYPE(
point) :: origpt
884 TYPE(
point),
DIMENSION(4) :: points
885 REAL,
DIMENSION(4) :: distsqrd
908 iloop:
DO i=0, dimi-2
925 IF ( indxi > 0 )
THEN
926 jloop:
DO j=jstart, jend, nextj
936 valid:
IF ( (indxi <= 0 .OR. dimi-1 < indxi) .OR. &
937 & (indxj <= 0 .OR. dimj-1 < indxj) )
THEN
976 SELECT CASE (minloc(distsqrd,1))
1018 REAL,
INTENT(in) :: lat
1022 REAL,
INTENT(in) :: lon
1039 latlon2xyz%x = radius * sin(theta) * cos(phi)
1040 latlon2xyz%y = radius * sin(theta) * sin(phi)
1054 TYPE(
point),
INTENT(in) :: pt1, pt2
1057 & (pt1%y-pt2%y)**2 +&
1067 REAL,
INTENT(in) :: lat1, lat2, lon1, lon2
1069 REAL :: theta1, theta2
1075 deltalambda =
deg2rad(lon2-lon1)
1076 deltatheta =
deg2rad(lat2-lat1)
1078 gcirdistance = radius * 2. * asin(sqrt((sin(deltatheta/2.))**2 + cos(theta1)*cos(theta2)*(sin(deltalambda/2.))**2))
1090 REAL,
INTENT(IN) :: lat
1091 REAL,
INTENT(IN) :: lon
1092 INTEGER,
INTENT(OUT) :: minI
1093 INTEGER,
INTENT(OUT) :: minJ
1094 REAL,
INTENT(OUT) :: minimum_distance
1105 IF (lat .EQ. 90.0)
THEN
1107 ELSEIF (lat .EQ. -90.0)
THEN
1120 minimum_distance = 2.0*radius*3.141592653
1127 IF (dist .LT. minimum_distance)
THEN
1142 minimum_distance = dist
1148 IF (mini .EQ. 0 .OR. minj .EQ. 0)
THEN
1149 call error_mesg(
"find_nearest_agrid_index", &
1150 "A minimum distance was not found.", &
1156 END MODULE diag_grid_mod
pure elemental real function distancesqrd(pt1, pt2)
Find the distance between two points in the Cartesian coordinate space.
pure integer function, dimension(2) find_equator_index_agrid(lat, lon)
Return the closest index (i,j) to the given (lat,lon) point.
subroutine, public diag_grid_init(domain, glo_lat, glo_lon, aglo_lat, aglo_lon)
Send the global grid to the diag_manager_mod for regional output.
pure elemental type(point) function latlon2xyz(lat, lon)
Return the (x,y,z) position of a given (lat,lon) point.
pure elemental real function gcirdistance(lat1, lon1, lat2, lon2)
Find the distance, along the geodesic, between two points.
type(diag_global_grid_type) diag_global_grid
Variable to hold the global grid data.
logical diag_grid_initialized
Indicates if the diag_grid_mod has been initialized.
subroutine, public diag_grid_end()
Unallocate the diag_global_grid variable.
pure elemental real function deg2rad(angle)
Convert an angle in degrees to radians.
subroutine, public get_local_indexes2(lat, lon, iindex, jindex)
Find the indices of the nearest grid point of the a-grid to the specified (lon,lat) location on the l...
subroutine, public get_local_indexes(latStart, latEnd, lonStart, lonEnd, istart, iend, jstart, jend)
Find the local start and local end indexes on the local PE for regional output.
pure integer function, dimension(2) find_pole_index_agrid(lat, lon)
Return the closest index (i,j) to the given (lat,lon) point.
subroutine find_nearest_agrid_index(lat, lon, minI, minJ, minimum_distance)
Find the i,j indices and distance of the a-grid point nearest to the inputted lat,...
pure elemental real function rad2deg(angle)
Private type to hold the model's global grid data, and other grid information for use in this module.
Private type to hold the corresponding (x,y,z) location for a (lat,lon) location.
subroutine, public write_version_number(version, tag, unit)
Prints to the log file (or a specified unit) the version id string and tag name.
subroutine, public error_mesg(routine, message, level)
Print notes, warnings and error messages; terminates program for warning and error messages....
integer function, dimension(size(domain%tile_id(:))) mpp_get_tile_id(domain)
Returns the tile_id on current pe.
integer function mpp_get_ntile_count(domain)
Returns number of tiles in mosaic.
Retrieve the entire array of compute domain extents associated with a decomposition.
The domain2D type contains all the necessary information to define the global, compute and data domai...
integer function mpp_npes()
Returns processor count for current pelist.
integer function mpp_pe()
Returns processor ID.
Reduction operations. Find the max of scalar a from the PEs in pelist result is also automatically br...
Reduction operations. Find the min of scalar a from the PEs in pelist result is also automatically br...