WAVEWATCH III  beta 0.0.1
Functions/Subroutines | Variables
scrip_remap_conservative Module Reference

Functions/Subroutines

subroutine remap_conserv (l_master, l_test)
 
subroutine cellblock_integrate (ibegin, iend, grid_num, phi_or_theta)
 
subroutine cell_integrate (cell_add, grid_num, phi_or_theta)
 
subroutine modify_polar_cell (ncorners, nalloc, cell_corner_lat, cell_corner_lon)
 
subroutine intersection (seg_cell_id, seg_grid_id, beglat, beglon, endlat, endlon, begseg, begedge, cell_id, ncorners, cell_corner_lat, cell_corner_lon, cell_grid_id, intrsct_lat, intrsct_lon, intedge, sinang2, lcoinc, lthresh)
 
subroutine pole_intersection (location, ncorners, cell_corners_lat, cell_corners_lon, cell_grid_id, beglat, beglon, endlat, endlon, begseg, begedge, intedge, intrsct_lat, intrsct_lon, sinang2, lcoinc, lthresh)
 
subroutine line_integral (phi_or_theta, weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
 
subroutine line_integral_phi (weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
 
subroutine line_integral_theta (weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
 
subroutine store_link_cnsrv (add1, add2, weights)
 
subroutine locate_segstart (cell_grid_num, cell, beglat, beglon, endlat, endlon, offset, srch_grid_num, cont_cell, lboundary, edgeid)
 
subroutine locate_point (ptlat, ptlon, cell, cell_grid_num, srch_grid_num, cont_cell, lboundary, edgeid)
 
subroutine ptincell (ptlat, ptlon, cell_add, ncorners, cell_corner_lat, cell_corner_lon, cell_center_lat, cell_center_lon, cell_grid_id, inpoly, lboundary, edgeid)
 
subroutine ptinpoly (ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, inpoly, lboundary, edgeid)
 
subroutine ptinpolarpoly (ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, whichpole, inpoly, lboundary, edgeid)
 
subroutine ptinpolygen (ptx, pty, ncorners, cell_corner_x, cell_corner_y, cell_center_x, cell_center_y, latlon, inpoly, lboundary, edgeid)
 
subroutine ptinpolygen2 (ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, inpoly, lboundary, edgeid)
 
subroutine get_srch_cells (cell_add, cell_grid_num, srch_grid_num, num_srch_cells, srch_add, srch_corners, srch_corner_lat, srch_corner_lon, srch_center_lat, srch_center_lon)
 
subroutine find_adj_cell (cell_add, edge_id, cell_grid_num, adj_add)
 
subroutine converge_to_bdry (cell_add, cell_grid_num, ncorners, cell_corner_lat, cell_corner_lon, cell_center_lat, cell_center_lon, inpt_lat, inpt_lon, outpt_lat, outpt_lon, bpt_lat, bpt_lon, bedgeid)
 

Variables

integer(scrip_i4) nthreads =2
 
integer(scrip_i4), save avoid_pole_count = 0
 
real(scrip_r8), save avoid_pole_offset = tiny
 
integer(scrip_i4), dimension(:,:), allocatable, save link_add1
 
integer(scrip_i4), dimension(:,:), allocatable, save link_add2
 
logical(scrip_logical), save first_call_store_link_cnsrv = .true.
 
logical(scrip_logical), save first_call_locate_segstart = .true.
 
integer(scrip_i4), save last_cell_locate_segstart =0
 
integer(scrip_i4), save last_cell_grid_num_locate_segstart =0
 
integer(scrip_i4), save last_srch_grid_num_locate_segstart =0
 
integer(scrip_i4), save num_srch_cells_locate_segstart =0
 
integer(scrip_i4), save srch_corners_locate_segstart
 
integer(scrip_i4), dimension(:), allocatable, save srch_add_locate_segstart
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_locate_segstart
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_locate_segstart
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_locate_segstart
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_locate_segstart
 
logical(scrip_logical), save first_call_locate_point = .true.
 
integer(scrip_i4), save last_cell_locate_point =0
 
integer(scrip_i4), save last_cell_grid_num_locate_point =0
 
integer(scrip_i4), save last_srch_grid_num_locate_point =0
 
integer(scrip_i4), save num_srch_cell_locate_points =0
 
integer(scrip_i4), save srch_corners_locate_point
 
integer(scrip_i4), dimension(:), allocatable, save srch_add_locate_point
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_locate_point
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_locate_point
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_locate_point
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_locate_point
 
integer(scrip_i4), save num_srch_cells_loc_get_srch_cells
 
integer(scrip_i4), save srch_corners_loc_get_srch_cells
 
integer(scrip_i4), dimension(:), allocatable, save srch_add_loc_get_srch_cells
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_loc_get_srch_cells
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_loc_get_srch_cells
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_loc_get_srch_cells
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_loc_get_srch_cells
 
integer(scrip_i4), save last_cell_add_get_srch_cells
 
integer(scrip_i4), save last_cell_grid_num_get_srch_cells
 
integer(scrip_i4), save last_srch_grid_num_get_srch_cells
 
logical(scrip_logical), save first_call_get_srch_cells =.true.
 
logical(scrip_logical), save first_call_find_adj_cell =.true.
 
integer(scrip_i4), save last_cell_find_adj_cell
 
integer(scrip_i4), save last_cell_grid_num_find_adj_cell
 
integer(scrip_i4), save num_srch_cells_find_adj_cell
 
integer(scrip_i4), save srch_corners_find_adj_cell
 
integer(scrip_i4), dimension(:), allocatable, save srch_add_find_adj_cell
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_find_adj_cell
 
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_find_adj_cell
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_find_adj_cell
 
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_find_adj_cell
 

Function/Subroutine Documentation

◆ cell_integrate()

subroutine scrip_remap_conservative::cell_integrate ( integer (scrip_i4)  cell_add,
integer (scrip_i4)  grid_num,
integer (scrip_i4)  phi_or_theta 
)

Definition at line 1005 of file scrip_remap_conservative.f.

1005 
1006 !-----------------------------------------------------------------------
1007 !
1008 ! Integrate around cell while finding intersecting with opposite
1009 ! grid cells and finding segments of cell boundary lying in cells
1010 ! of opposite grid
1011 !
1012 !-----------------------------------------------------------------------
1013 
1014 !-----------------------------------------------------------------------
1015 !
1016 ! Input variables
1017 !
1018 !-----------------------------------------------------------------------
1019 
1020  integer (SCRIP_i4) ::
1021  & cell_add, ! cell to be processed
1022  & grid_num, ! grid that the cell belongs to
1023  & phi_or_theta ! Integration var :
1024  ! phi (lon) or theta (lat)
1025 
1026 
1027 !-----------------------------------------------------------------------
1028 !
1029 ! local variables
1030 !
1031 !-----------------------------------------------------------------------
1032 
1033  integer (SCRIP_i4), parameter ::
1034  & max_subseg = 500 ! max number of subsegments per segment
1035  ! to prevent infinite loop
1036 
1037 
1038  integer (SCRIP_i4) ::
1039  & i, inext, !
1040  & j, jnext, ! generic counters
1041  & ic, k, ns, !
1042  & n, next_n, !
1043  & nwgt, it, !
1044  & oppcell_add, ! Cell from opposite grid we are
1045  ! intersecting
1046  & opp_grid_num, ! Index of opposite grid (2,1)
1047  & min_add, ! addresses for restricting search of
1048  & max_add, ! destination grid
1049  & corner, ! corner of cell that segment starts
1050  ! from
1051  & next_corn, ! corner of cell that segment ends on
1052  & nseg, ! number of segments to use to represent
1053  ! edges near the pole
1054  & num_subseg, ! number of subsegments
1055  & bedgeid1, !
1056  & bedgeid2, ! ID of edge that a point is on
1057  & bedgeid3, !
1058  & intedge, ! ID of intersected edge
1059  & last_add, ! Address of last cell we were in
1060  & next_add, ! Address of next cell we will go into
1061  & adj_add ! Address of cell adjacent to current
1062  ! one
1063 
1064  logical (SCRIP_logical) ::
1065  & lcoinc, ! Are segments coincident?
1066  & lrevers, ! Are we integrating segment in reverse?
1067  & lboundary1,
1068  & lboundary2, ! Is point is on cell boundary?
1069  & lboundary3,
1070  & last_lboundary, ! Is last point is on cell bdry?
1071  & loutside, ! Is point outside the grid?
1072  & lthresh, ! Has segment crossed threshold?
1073  & srch_success, ! Was search for segment start
1074  ! successful?
1075  & intrsct_success, ! Was intersection of segment with
1076  ! opposite grid successful?
1077  & inpoly, ! Is point is in polygon
1078  & last_endpt_inpoly, ! Was end point of last segment in cell
1079  & last_endpt_onedge, ! Was end point of last segment on edge
1080  ! of cell
1081  & lstuck, ! Is the walk stuck inside a cell
1082  & seg_outside, ! Is segment completely outside the grid
1083  & bndedge, ! Is segment on the boundary of the grid
1084  & search, ! Do we have to search to locate point
1085  ! in grid
1086  & inpolar, ! Are we in the polar region?
1087  & special_cell, ! Is this a special cell
1088  ! (only 1 vtx at pole)
1089  & L_exit_do ! Do we need to escape from infinite
1090  ! loop? (NRL)
1091 
1092  real (SCRIP_r8) ::
1093  & intrsct_lat, ! lat of next intersection point
1094  & intrsct_lon, ! lon of next intersection point
1095  & beglat, beglon, ! start point of current sub seg
1096  & endlat, endlon, ! endpoint of current seg
1097  ! (chg to endseg?)
1098  & endlat1, endlon1, ! endpoint of current subseg
1099  & norm_factor ! factor for normalizing wts
1100 
1101  real (SCRIP_r8), dimension(2) ::
1102  & begseg ! begin lat/lon for full segment
1103 
1104  real (SCRIP_r8), dimension(6) ::
1105  & weights, ! local wgt array
1106  & rev_weights ! Weights for grid1 and grid2 flipped
1107 
1108  real (SCRIP_r8) ::
1109  & vec1_lat, vec1_lon, ! vectors, products
1110  & vec2_lat, vec2_lon, ! used in grid search
1111  & vec1_len, dp,
1112  & midlat, midlon, ! Midpoint of segment
1113  & tmplat, tmplon,
1114  & srchpt_lat, ! Search point (offset from seg. start)
1115  & srchpt_lon,
1116  & offset, delta, ! Offset and offset increase for search
1117  & sinang2, ! Square of sine of angle b/w two
1118  ! segments
1119  & dist2, ! Square of distance b/w two points
1120  & fullseg_len2, ! Square of full segment length
1121  & partseg_len2, ! Square of length of segment integrated
1122  ! so far
1123  & fullseg_dlat, ! Lat diff of full segment endpoints
1124  & fullseg_dlon, ! Lon diff of full segment endpoints
1125  & prevlon,
1126  & nextlon,
1127  & pole_lat,
1128  & cell_center_lat,
1129  & cell_center_lon,
1130  & oppcell_center_lat,
1131  & oppcell_center_lon
1132 
1133  real (SCRIP_r8), dimension(:), allocatable ::
1134  & cell_corner_lat, ! Local copies of cell coordinates
1135  & cell_corner_lon, ! May be augmented for computational
1136  ! reasons
1137  & oppcell_corner_lat,
1138  & oppcell_corner_lon
1139 
1140  integer (SCRIP_i4) ::
1141  & ncorners, ! Number of corners in local copy of cell
1142  & ncorners_opp, ! Number of corners in local copy of oppcell
1143  & nalloc, ! Allocation for the cell_corner_* array
1144  & nalloc_opp ! Allocation for the oppcell_corner_* array
1145 
1146  real (SCRIP_r8) ::
1147  & tmpwt1, tmpwt2
1148 
1149  integer (SCRIP_i4) ::
1150  & ncorners_at_pole,
1151  & previdx,
1152  & nextidx
1153 
1154 
1155  if (grid_num .eq. 1) then
1156 
1157  !***
1158  !*** Set up a local copy of the cell with room to add
1159  !*** degenerate edges
1160  !***
1161 
1162  ncorners = grid1_corners
1163  nalloc = min(ncorners + 2,
1164  & size(grid1_corner_lat(:,1)))
1165  allocate (cell_corner_lat(nalloc),
1166  & cell_corner_lon(nalloc))
1167 
1168  do corner = 1, ncorners
1169  cell_corner_lat(corner) = grid1_corner_lat(corner,cell_add)
1170  cell_corner_lon(corner) = grid1_corner_lon(corner,cell_add)
1171  enddo
1172 
1173  cell_center_lat = grid1_center_lat(cell_add)
1174  cell_center_lon = grid1_center_lon(cell_add)
1175 
1176  special_cell = special_polar_cell1(cell_add)
1177 
1178 
1179  !***
1180  !*** Also, allocate storage for the cell from the opposite grid
1181  !***
1182 
1183  opp_grid_num = 2
1184  ncorners_opp = grid2_corners
1185  nalloc_opp = ncorners_opp+2
1186  allocate (oppcell_corner_lat(nalloc_opp),
1187  & oppcell_corner_lon(nalloc_opp))
1188 
1189  else
1190 
1191  !***
1192  !*** Set up the cell info with room to add degenerate edges
1193  !***
1194 
1195  ncorners = grid2_corners
1196  nalloc = min(ncorners + 2,
1197  & size(grid2_corner_lat(:,1)))
1198  allocate (cell_corner_lat(nalloc),
1199  & cell_corner_lon(nalloc))
1200 
1201  do corner = 1, ncorners
1202  cell_corner_lat(corner) = grid2_corner_lat(corner,cell_add)
1203  cell_corner_lon(corner) = grid2_corner_lon(corner,cell_add)
1204  enddo
1205 
1206  cell_center_lat = grid2_center_lat(cell_add)
1207  cell_center_lon = grid2_center_lon(cell_add)
1208 
1209  special_cell = special_polar_cell2(cell_add)
1210 
1211  !***
1212  !*** Also, allocate storage for the cell from the opposite grid
1213  !***
1214 
1215  opp_grid_num = 1
1216  ncorners_opp = grid1_corners
1217  nalloc_opp = ncorners_opp + 2
1218  allocate (oppcell_corner_lat(nalloc_opp),
1219  & oppcell_corner_lon(nalloc_opp))
1220 
1221  endif
1222 
1223  if (special_cell) then
1224 
1225  !***
1226  !*** Special cell with only one corner at the pole Such cells
1227  !*** can have an artificially extreme distortion of the edges
1228  !*** when mapped to the Lambert projection because of the span
1229  !*** of longitudes on the edges So we will augment such cells
1230  !*** with degenerate edges at the pole so that the edges coming
1231  !*** off the pole will actually have the same longitude values
1232  !*** at both ends
1233  !***
1234  !*** lon_p lon_p+ lon_p lon_p-
1235  !*** pi/2 pi/2 pi/2 pi/2
1236  !*** * *--------*------*
1237  !*** / \ | |
1238  !*** / \ | |
1239  !*** / \ | |
1240  !*** * * * *
1241  !*** lon_p+ lon_p- lon_p+ lon_p-
1242  !*** lat_p+ lat_p- lat_p+ lat_p-
1243  !***
1244 
1245  call modify_polar_cell(ncorners,nalloc,cell_corner_lat,
1246  & cell_corner_lon)
1247 
1248  endif
1249 
1250  !***
1251  !*** Cell info set up - Now process the cell
1252  !***
1253 
1254  do corner = 1, ncorners
1255  next_corn = mod(corner,ncorners) + 1
1256 
1257  !***
1258  !*** define endpoints of the current segment
1259  !***
1260 
1261  beglat = cell_corner_lat(corner)
1262  beglon = cell_corner_lon(corner)
1263  endlat = cell_corner_lat(next_corn)
1264  endlon = cell_corner_lon(next_corn)
1265  lrevers = .false.
1266 
1267  !***
1268  !*** if this is a constant-longitude segment, skip the rest
1269  !*** since the line integral contribution will be zero.
1270  !***
1271 
1272  if ((phi_or_theta == 1 .and. endlon == beglon) .or.
1273  & (phi_or_theta == 2 .and. endlat == beglat)) cycle
1274 
1275  !***
1276  !*** to ensure exact path taken during both
1277  !*** sweeps, always integrate segments in the same
1278  !*** direction (SW to NE).
1279  !***
1280 
1281  if ((endlat < beglat) .or.
1282  & (endlat == beglat .and. endlon < beglon)) then
1283  tmplat = beglat
1284  beglat = endlat
1285  endlat = tmplat
1286  tmplon = beglon
1287  beglon = endlon
1288  endlon = tmplon
1289  lrevers = .not. lrevers
1290  endif
1291 
1292  !*** But if one of the segment ends is in the polar region,
1293  !*** we want to start from that (makes some logic easier)
1294 
1295  if ((beglat < north_thresh .and. endlat > north_thresh) .or.
1296  & (beglat > south_thresh .and. endlat < south_thresh))
1297  & then
1298  tmplat = beglat
1299  beglat = endlat
1300  endlat = tmplat
1301  tmplon = beglon
1302  beglon = endlon
1303  endlon = tmplon
1304  lrevers = .not. lrevers
1305  endif
1306 
1307  begseg(1) = beglat
1308  begseg(2) = beglon
1309 
1310  fullseg_dlat = endlat-beglat
1311  fullseg_dlon = endlon-beglon
1312  if (fullseg_dlon > pi) fullseg_dlon = fullseg_dlon - pi2
1313  if (fullseg_dlon < -pi) fullseg_dlon = fullseg_dlon + pi2
1314  fullseg_len2 = fullseg_dlat*fullseg_dlat +
1315  & fullseg_dlon*fullseg_dlon
1316 
1317  partseg_len2 = 0.0
1318 
1319  !***
1320  !*** Is this an edge on the boundary of the grid or
1321  !*** on the boundary of the active cells
1322  !***
1323 
1324 ! Commented out by MD
1325 ! call find_adj_cell(cell_add, corner, grid_num, adj_add)
1326 ! if (grid_num .eq. 1) then
1327 ! if (adj_add .eq. 0 .or. .not. grid1_mask(adj_add)) then
1328 ! bndedge = .true.
1329 ! else
1330 ! bndedge = .false.
1331 ! endif
1332 ! else
1333 ! if (adj_add .eq. 0 .or. .not. grid2_mask(adj_add)) then
1334 ! bndedge = .true.
1335 ! else
1336 ! bndedge = .false.
1337 ! endif
1338 ! endif
1339 
1340  call find_adj_cell(cell_add, corner, grid_num, adj_add)
1341  bndedge = .false.
1342  if (grid_num .eq. 1) then
1343  if (adj_add .eq. 0) then
1344  bndedge = .true.
1345  else
1346  if (.not. grid1_mask(adj_add)) then
1347  bndedge = .true.
1348  endif
1349  endif
1350  else
1351  if (adj_add .eq. 0) then
1352  bndedge = .true.
1353  else
1354  if (.not. grid2_mask(adj_add)) then
1355  bndedge = .true.
1356  endif
1357  endif
1358  endif
1359 
1360  !***
1361  !*** integrate along this segment, detecting intersections
1362  !*** and computing the line integral for each sub-segment
1363  !***
1364 
1365  if (beglat .gt. north_thresh .or. beglat .lt. south_thresh)
1366  & then
1367  nseg = npseg ! Use multiple subsegments near the pole
1368  inpolar = .true.
1369  else
1370  nseg = 1
1371  inpolar = .false.
1372  endif
1373 
1374 
1375  last_add = 0
1376  last_lboundary = .false.
1377  last_endpt_inpoly = .false.
1378  last_endpt_onedge = .false.
1379  next_add = 0
1380  search = .true.
1381  ns = 1
1382 
1383 ! outer "do while"
1384 
1385  do while (beglat /= endlat .or. beglon /= endlon)
1386 
1387  l_exit_do=.false. !NRL
1388 
1389  if ((ns .eq. nseg) .or. (inpolar .eqv. .false.)) then
1390  !
1391  ! Last subseg or out of the polar region
1392  ! Go directly to end of segment
1393  !
1394  endlat1 = endlat
1395  endlon1 = endlon
1396  else
1397  endlat1 = begseg(1) + ns*(fullseg_dlat)/nseg
1398  endlon1 = begseg(2) + ns*(fullseg_dlon)/nseg
1399  endif
1400 
1401  num_subseg = 0
1402 
1403 ! inner "do while"
1404 
1405  do while (beglat /= endlat1 .or. beglon /= endlon1)
1406 
1407  !***
1408  !*** If we integrated to the end or just past it (due to
1409  !*** numerical errors), we are done with this segment
1410  !***
1411 
1412 !NRL see notes below re: infinite "do while" loop
1413  l_exit_do=.false. !NRL
1414  if (partseg_len2 .ge. fullseg_len2) then
1415  write(*,*).ge.'partseg_len2 fullseg_len2'
1416  write(*,*)'beglat,beglon = ',beglat,beglon
1417  write(*,*)'endlat,endlon = ',endlat,endlon
1418  write(*,*)'endlat1,endlon1 = ',endlat1,endlon1
1419  write(*,*)'exiting inner do while loop'
1420  l_exit_do=.true. !NRL
1421  exit
1422  end if
1423 
1424  !******************************************************
1425  !*** Try to find which cell of the opposite grid this
1426  !*** segment is starting in and where it is exiting this
1427  !*** cell
1428  !******************************************************
1429 
1430  vec1_lat = endlat1-beglat
1431  vec1_lon = endlon1-beglon
1432  if (vec1_lon > pi) vec1_lon = vec1_lon - pi2
1433  if (vec1_lon < -pi) vec1_lon = vec1_lon + pi2
1434  vec1_len = sqrt(vec1_lat*vec1_lat+vec1_lon*vec1_lon)
1435  vec1_lat = vec1_lat/vec1_len
1436  vec1_lon = vec1_lon/vec1_len
1437 
1438  offset = 100.0*tiny
1439  oppcell_add = 0
1440  delta = 10*tiny
1441  intrsct_success = .false.
1442  loutside = .false.
1443  lstuck = .false.
1444  lboundary1 = .false.
1445  lboundary2 = .false.
1446  lboundary3 = .false.
1447 
1448  do while (.not. intrsct_success)
1449 
1450  !*************************************************
1451  !*** Find out which cell the segment starts in
1452  !*************************************************
1453 
1454  srch_success = .false.
1455  if (search) then
1456 
1457  !***
1458  !*** Offset the start point in ever increasing
1459  !*** amounts until we are able to reliably locate
1460  !*** the point in a cell of grid2. Inability to locate
1461  !*** the point causes the offset amount to increase
1462 
1463  it = 0
1464  do while (.not. srch_success)
1465 
1466  srchpt_lat = beglat + offset*vec1_lat
1467  srchpt_lon = beglon + offset*vec1_lon
1468 
1469  call locate_point(srchpt_lat, srchpt_lon,
1470  & cell_add, grid_num, opp_grid_num,
1471  & oppcell_add, lboundary1, bedgeid1)
1472 
1473  if (oppcell_add .eq. 0) then
1474  loutside = .true.
1475 ! lcoinc added by MD
1476  lcoinc = .false.
1477  exit ! exit the search loop
1478  else
1479  if (oppcell_add .ne. last_add .or. lthresh)
1480  & then
1481  srch_success = .true.
1482  else
1483  offset = offset + delta
1484  if (offset .ge. vec1_len) then
1485  exit
1486  endif
1487  if (it .gt. 3) then
1488  delta = 2.0*delta
1489  it = 0
1490  endif
1491  endif
1492  endif
1493 
1494  it = it + 1
1495  enddo ! do while (.not. srch_success)
1496 
1497  else
1498  if (last_endpt_inpoly) then
1499 
1500  !*** We know the grid cell the end of the last
1501  !*** segment (which is the beginning of this
1502  !*** segment)
1503 
1504  oppcell_add = last_add
1505  lboundary1 = last_lboundary
1506 
1507  else if (next_add .ne. 0) then
1508 
1509  !*** We know the edge of the grid2 cell that the
1510  !*** last segment intersected, so we move into
1511  !*** the adjacent cell
1512 
1513  oppcell_add = next_add
1514  lboundary1 = .true.
1515 
1516  endif
1517 
1518  srch_success = .true.
1519 
1520  endif
1521 
1522  !*****************************************************
1523  !*** Find where the segment exits this cell, if at all
1524  !*****************************************************
1525 
1526  if (srch_success) then
1527 
1528  !***
1529  !*** First setup local copy of oppcell with room for
1530  !*** adding degenerate edges
1531  !***
1532 
1533  if (grid_num .eq. 1) then
1534  ncorners_opp = grid2_corners
1535  do i = 1, ncorners_opp
1536  oppcell_corner_lat(i) =
1537  & grid2_corner_lat(i,oppcell_add)
1538  oppcell_corner_lon(i) =
1539  & grid2_corner_lon(i,oppcell_add)
1540  enddo
1541  oppcell_center_lat =
1542  & grid2_center_lat(oppcell_add)
1543  oppcell_center_lon =
1544  & grid2_center_lon(oppcell_add)
1545 
1546  special_cell = special_polar_cell2(oppcell_add)
1547  else
1548  ncorners_opp = grid1_corners
1549  do i = 1, ncorners_opp
1550  oppcell_corner_lat(i) =
1551  & grid1_corner_lat(i,oppcell_add)
1552  oppcell_corner_lon(i) =
1553  & grid1_corner_lon(i,oppcell_add)
1554  enddo
1555  oppcell_center_lat =
1556  & grid1_center_lat(oppcell_add)
1557  oppcell_center_lon =
1558  & grid1_center_lon(oppcell_add)
1559 
1560  special_cell = special_polar_cell1(oppcell_add)
1561  endif
1562 
1563  if (special_cell) then
1564  call modify_polar_cell(ncorners_opp, nalloc_opp,
1565  & oppcell_corner_lat, oppcell_corner_lon)
1566  endif
1567 
1568  !***
1569  !*** First see if the segment end is
1570  !*** in the same cell
1571  !***
1572 
1573  call ptincell(endlat1,endlon1, oppcell_add,
1574  & ncorners_opp,
1575  & oppcell_corner_lat,oppcell_corner_lon,
1576  & oppcell_center_lat,oppcell_center_lon,
1577  & opp_grid_num,inpoly,
1578  & lboundary2,bedgeid2)
1579 
1580  if (inpoly) then
1581  intrsct_lat = endlat1
1582  intrsct_lon = endlon1
1583  intrsct_success = .true.
1584  search = .false.
1585  next_add = 0
1586  last_add = oppcell_add ! for next subseg
1587  last_lboundary = lboundary2
1588  last_endpt_inpoly = .true.
1589 
1590  if (lboundary1 .and. lboundary2) then
1591 
1592  !*** This is a edge on the boundary of the
1593  !*** active mesh and both of its endpoints
1594  !*** are on the boundary of the containing
1595  !*** cell. Check if the the segment is also
1596  !*** on the boundary
1597 
1598  midlat = (beglat+endlat1)/2.0
1599  if (abs(beglon-endlon1) .ge. pi) then
1600  midlon = (beglon+endlon1)/2.0 - pi
1601  else
1602  midlon = (beglon+endlon1)/2.0
1603  endif
1604 
1605  call ptincell(midlat,midlon, oppcell_add,
1606  & ncorners_opp,
1607  & oppcell_corner_lat, oppcell_corner_lon,
1608  & oppcell_center_lat, oppcell_center_lon,
1609  & opp_grid_num, inpoly, lboundary3,
1610  & bedgeid3)
1611 
1612  if (inpoly .and. lboundary3) then
1613  lcoinc = .true.
1614  intedge = bedgeid3
1615  endif
1616 
1617  else
1618  lcoinc = .false.
1619  endif
1620 
1621  else
1622 
1623  !***
1624  !*** Do an intersection to find out where the
1625  !*** segment exits the cell
1626  !***
1627 
1628  call intersection(cell_add,grid_num,
1629  & beglat, beglon, endlat1, endlon1,
1630  & begseg,
1631  & bedgeid1,
1632  & oppcell_add, ncorners_opp,
1633  & oppcell_corner_lat, oppcell_corner_lon,
1634  & opp_grid_num,
1635  & intrsct_lat, intrsct_lon, intedge,
1636  & sinang2, lcoinc, lthresh)
1637 
1638  if (intedge /= 0) then
1639  intrsct_success = .true.
1640  last_add = oppcell_add ! for next subseg
1641  last_endpt_onedge = .true.
1642  last_endpt_inpoly = .false.
1643  last_lboundary = .true.
1644 
1645  if (.not. lthresh) then
1646  call find_adj_cell(oppcell_add,intedge,
1647  & opp_grid_num,next_add)
1648  if (next_add .ne. 0) then
1649  search = .false.
1650  else
1651  search = .true.
1652  endif
1653  else
1654  search = .true.
1655  endif
1656  endif
1657 
1658  endif
1659 
1660  if (.not. intrsct_success) then
1661 
1662  !*** Offset point and try again
1663 
1664  search = .true.
1665  delta = 2.0*delta
1666  offset = offset + delta
1667  if (offset .gt. vec1_len) then
1668 
1669  ! Punt - exit the intersection loop
1670 
1671  intrsct_lat = endlat1
1672  intrsct_lon = endlon1
1673  last_add = 0
1674  last_lboundary = .false.
1675  exit
1676 
1677  endif
1678  endif
1679 
1680 !NRL if (lcoinc .and. .not. bndedge) then
1681 
1682  if (lcoinc .and. .not. bndedge !NRL
1683  & .and. intedge /= 0) then !NRL
1684 
1685  !***
1686  !*** Segment is coincident with edge of other grid
1687  !*** which means it could belong to one of 2 cells
1688  !*** Choose the cell such that edge that is
1689  !*** coincident with the segment is in the same
1690  !*** dir as the segment
1691 
1692  i = intedge
1693  inext = mod(i,ncorners_opp)+1
1694  vec2_lat = oppcell_corner_lat(inext) -
1695  & oppcell_corner_lat(i)
1696  vec2_lon = oppcell_corner_lon(inext) -
1697  & oppcell_corner_lon(i)
1698 
1699  if (vec2_lon > pi) vec2_lon = vec2_lon - pi2
1700  if (vec2_lon < -pi) vec2_lon = vec2_lon + pi2
1701 
1702  dp = vec1_lat*vec2_lat + vec1_lon*vec2_lon
1703 
1704  if ((.not. lrevers .and. dp .lt. 0) .or.
1705  & (lrevers .and. dp .gt. 0)) then
1706 
1707  !*** Integrals from this segment must be
1708  !*** assigned to the adjacent cell of
1709  !*** opcell_add but only if such an adjacent
1710  !*** cell exists
1711 
1712  call find_adj_cell(oppcell_add, intedge,
1713  & opp_grid_num, adj_add)
1714 
1715  if (adj_add .gt. 0) then
1716  oppcell_add = adj_add
1717 
1718  if (grid_num .eq. 1) then
1719  ncorners_opp = grid2_corners
1720  do i = 1, ncorners_opp
1721  oppcell_corner_lat(i) =
1722  & grid2_corner_lat(i,oppcell_add)
1723  oppcell_corner_lon(i) =
1724  & grid2_corner_lon(i,oppcell_add)
1725  enddo
1726  oppcell_center_lat =
1727  & grid2_center_lat(oppcell_add)
1728  oppcell_center_lon =
1729  & grid2_center_lon(oppcell_add)
1730 
1731  special_cell =
1732  & special_polar_cell2(oppcell_add)
1733  else
1734  ncorners_opp = grid1_corners
1735  do i = 1, ncorners_opp
1736  oppcell_corner_lat(i) =
1737  & grid1_corner_lat(i,oppcell_add)
1738  oppcell_corner_lon(i) =
1739  & grid1_corner_lon(i,oppcell_add)
1740  enddo
1741  oppcell_center_lat =
1742  & grid1_center_lat(oppcell_add)
1743  oppcell_center_lon =
1744  & grid1_center_lon(oppcell_add)
1745 
1746  special_cell =
1747  & special_polar_cell1(oppcell_add)
1748  endif
1749 
1750  if (special_cell) then
1751  call modify_polar_cell(ncorners_opp,
1752  & nalloc_opp, oppcell_corner_lat,
1753  & oppcell_corner_lon)
1754  endif
1755 
1756  endif
1757 
1758  endif
1759 
1760  endif
1761 
1762  else
1763 
1764  !***
1765  !*** Could not locate a viable cell for the segment
1766  !*** start
1767  !***
1768 
1769  if (oppcell_add .eq. 0) then
1770  loutside = .true.
1771 ! lcoinc added by MD
1772  lcoinc = .false.
1773 
1774  !***
1775  !*** Take baby steps to see if any part of the
1776  !*** segment is inside a cell of the other grid
1777  !***
1778 
1779  seg_outside = .false.
1780  delta = vec1_len/100.00
1781  offset = delta
1782  do while (.not. srch_success)
1783 
1784  srchpt_lat = beglat + offset*vec1_lat
1785  srchpt_lon = beglon + offset*vec1_lon
1786 
1787  call locate_point(srchpt_lat, srchpt_lon,
1788  & cell_add, grid_num, opp_grid_num,
1789  & oppcell_add, lboundary1, bedgeid1)
1790 
1791  if (oppcell_add /= 0) then
1792  srch_success = .true.
1793 
1794  !***
1795  !*** Found a point of the segment in the
1796  !*** cell. Do a bisection method to find
1797  !*** the starting point of the segment
1798  !*** in the cell
1799  !***
1800 
1801  call converge_to_bdry(oppcell_add,
1802  & opp_grid_num, ncorners_opp,
1803  & oppcell_corner_lat,
1804  & oppcell_corner_lon,
1805  & oppcell_center_lat,
1806  & oppcell_center_lon,
1807  & srchpt_lat, srchpt_lon,
1808  & beglat, beglon,
1809  & intrsct_lat, intrsct_lon,
1810  & bedgeid1)
1811 
1812  search = .false.
1813  last_endpt_onedge = .true.
1814  next_add = oppcell_add
1815  last_lboundary = .true.
1816 
1817  oppcell_add = 0
1818 
1819  else
1820 
1821  offset = offset + delta
1822 
1823  if (offset .ge. vec1_len) then
1824 ! print *,
1825 ! & 'Segment fully outside grid2'
1826 ! print *, 'Segment of grid1_add',
1827 ! & grid1_add
1828 ! print *, beglat,beglon
1829 ! print *, endlat1,endlon1
1830 
1831  seg_outside = .true.
1832 
1833  intrsct_lat = endlat1
1834  intrsct_lon = endlon1
1835 
1836  search = .true.
1837  last_add = 0
1838  last_lboundary = .false.
1839 
1840  exit ! leave search loop
1841  endif
1842  endif
1843 
1844  enddo
1845 
1846  ! int. loop
1847  if (srch_success .or. seg_outside) exit
1848 
1849  else
1850 
1851  if(is_master)then
1852  print *, 'Unable to move out of last cell'
1853  print *, 'Segment of edge ',corner,
1854  & ' of grid cell ',cell_add
1855  print *, 'Stuck in opposite grid cell ',
1856  & oppcell_add
1857  dist2 =
1858  & (endlat1-begseg(1))*(endlat1-begseg(1)) +
1859  & (endlon1-begseg(2))*(endlon1-begseg(2))
1860  print *, 'Fraction of segment left ',
1861  & vec1_len/sqrt(dist2)
1862  endif
1863  lstuck = .true.
1864 
1865  !***
1866  !*** Punt - just assign the rest of the segment
1867  !*** to the current cell it is stuck in by
1868  !*** tagging the segment endpoint as the
1869  !*** intersection point
1870  !***
1871 
1872  intrsct_lat = endlat1
1873  intrsct_lon = endlon1
1874 
1875  search = .true.
1876  last_add = 0
1877  last_lboundary = .false.
1878 
1879  endif
1880 
1881  exit ! exit the intersection loop
1882 
1883  endif ! if (srch_success) then ... else ....
1884 
1885  end do ! do while (.not. intrsct_success)
1886 
1887  !********************************************************
1888  !*** Compute the line integrals for this subsegment
1889  !********************************************************
1890 
1891  if (oppcell_add /= 0) then
1892  call line_integral(phi_or_theta, weights, num_wts,
1893  & beglon, intrsct_lon, beglat, intrsct_lat,
1894  & cell_center_lat, cell_center_lon,
1895  & oppcell_center_lat, oppcell_center_lon)
1896  else
1897  call line_integral(phi_or_theta, weights, num_wts,
1898  & beglon, intrsct_lon, beglat, intrsct_lat,
1899  & cell_center_lat, cell_center_lon,
1900  & cell_center_lat, cell_center_lon)
1901  endif
1902 
1903  !***
1904  !*** if integrating in reverse order, change
1905  !*** sign of weights
1906  !***
1907 
1908  if (lrevers) then
1909  weights = -weights
1910  endif
1911 
1912  !***
1913  !*** store the appropriate addresses and weights.
1914  !*** also add contributions to cell areas and centroids.
1915  !***
1916 
1917  if (grid_num .eq. 1) then
1918 
1919  if (oppcell_add /= 0) then
1920  if (grid1_mask(cell_add)) then
1921  call store_link_cnsrv(cell_add, oppcell_add,
1922  & weights)
1923 
1924 C$OMP CRITICAL(block1)
1925 !
1926 ! Could have another thread that found an intersection between that
1927 ! cell address and oppcell_add in which case it will try to write
1928 ! into this address - we have to block that until we are finished
1929 !
1930  grid1_frac(cell_add) =
1931  & grid1_frac(cell_add) + weights(1)
1932 
1933  grid2_frac(oppcell_add) =
1934  & grid2_frac(oppcell_add) + weights(num_wts+1)
1935 C$OMP END CRITICAL(block1)
1936  endif
1937 
1938  endif
1939 
1940 C$OMP CRITICAL(block2)
1941  grid1_area(cell_add) = grid1_area(cell_add) +
1942  & weights(1)
1943  grid1_centroid_lat(cell_add) =
1944  & grid1_centroid_lat(cell_add) + weights(2)
1945  grid1_centroid_lon(cell_add) =
1946  & grid1_centroid_lon(cell_add) + weights(3)
1947 C$OMP END CRITICAL(block2)
1948 
1949  else
1950 
1951  !*** swap weights because in store_link_cnsrv
1952  !*** we are always sending in grid1 weights first
1953  !*** and then grid2 weights
1954 
1955  do i = 1, num_wts
1956  rev_weights(num_wts+i) = weights(i)
1957  rev_weights(i) = weights(num_wts+i)
1958  enddo
1959 
1960  if (.not. lcoinc .and. oppcell_add /= 0) then
1961  if (grid1_mask(oppcell_add)) then
1962  call store_link_cnsrv(oppcell_add, cell_add,
1963  & rev_weights)
1964 
1965 C$OMP CRITICAL(block3)
1966 !
1967 ! Could have another thread that found an intersection between that
1968 ! cell address and oppcell_add in which case it will try to write
1969 ! into this address - we have to block that until we are finished
1970 !
1971  grid2_frac(cell_add) =
1972  & grid2_frac(cell_add) + weights(1)
1973 
1974  grid1_frac(oppcell_add) =
1975  & grid1_frac(oppcell_add) + weights(num_wts+1)
1976 C$OMP END CRITICAL(block3)
1977 
1978  endif
1979 
1980  endif
1981 
1982 C$OMP CRITICAL(block4)
1983  grid2_area(cell_add) = grid2_area(cell_add) +
1984  & weights(1)
1985  grid2_centroid_lat(cell_add) =
1986  & grid2_centroid_lat(cell_add) + weights(2)
1987  grid2_centroid_lon(cell_add) =
1988  & grid2_centroid_lon(cell_add) + weights(3)
1989 C$OMP END CRITICAL(block4)
1990  endif
1991 
1992  !***
1993  !*** reset beglat and beglon for next subsegment.
1994  !***
1995 
1996  beglat = intrsct_lat
1997  beglon = intrsct_lon
1998 
1999  !***
2000  !*** How far have we come from the start of the segment
2001  !***
2002 
2003  vec2_lat = intrsct_lat-begseg(1)
2004  vec2_lon = intrsct_lon-begseg(2)
2005  if (vec2_lon > pi) vec2_lon = vec2_lon - pi2
2006  if (vec2_lon < -pi) vec2_lon = vec2_lon + pi2
2007 
2008  partseg_len2 = vec2_lat*vec2_lat + vec2_lon*vec2_lon
2009 
2010  !***
2011  !*** prevent infinite loops if integration gets stuck
2012  !*** near cell or threshold boundary
2013  !***
2014 
2015  num_subseg = num_subseg + 1
2016  if (num_subseg > max_subseg) then
2017  print *,
2018  & 'integration stalled: num_subseg exceeded limit'
2019  print *, 'Cell ',cell_add
2020  print *, 'Edge ',corner
2021  print *, 'Grid ',1
2022  dist2 = (endlat1-begseg(1))*(endlat1-begseg(1)) +
2023  & (endlon1-begseg(2))*(endlon1-begseg(2))
2024  print *, 'Fraction of segment left ',
2025  & vec1_len/sqrt(dist2)
2026 ! exit ! Give up and exit
2027  stop ! Give up and stop
2028  endif
2029 
2030 ! inner "do while"
2031  end do ! do while (beglat /= endlat1 ...
2032 
2033 !NRL We add an exit to outer do similar to exit of inner do:
2034 !NRL This was an apparent bug: exit statement would escape
2035 !NRL inner do but then computation could not get out of
2036 !NRL outer do since beglat, beglon controlling outer do
2037 !NRL never changed b/c it never gets to the part of the
2038 !NRL code that changes beglat, beglon, b/c it keeps
2039 !NRL exiting inner do.
2040 
2041 !NRL This should happen very rarely, so we have a print
2042 !NRL statement to notify user.
2043 
2044  if (l_exit_do)then ! NRL
2045  write(*,*)'partseg_len2,fullseg_len2 = ', ! NRL
2046  & partseg_len2,fullseg_len2 ! NRL
2047  write(*,*)'exiting outer do while loop' ! NRL
2048  exit ! NRL
2049  endif ! NRL
2050 
2051  ns = ns + 1
2052  if ((beglat > 0 .and. beglat < north_thresh) .or.
2053  & (beglat < 0 .and. beglat > south_thresh))
2054  & then
2055  inpolar = .false.
2056  endif
2057 
2058 ! outer "do while"
2059  end do ! do while (beglat /= endlat ....
2060 
2061  call line_integral(phi_or_theta, weights, num_wts,
2062  & begseg(2), endlon, begseg(1), endlat,
2063  & cell_center_lat,
2064  & cell_center_lon,
2065  & cell_center_lat,
2066  & cell_center_lon)
2067 
2068  !***
2069  !*** end of segment
2070  !***
2071 
2072  end do ! do corner=....
2073 

References converge_to_bdry(), find_adj_cell(), scrip_grids::grid1_area, scrip_grids::grid1_center_lat, scrip_grids::grid1_center_lon, scrip_grids::grid1_centroid_lat, scrip_grids::grid1_centroid_lon, scrip_grids::grid1_corner_lat, scrip_grids::grid1_corner_lon, scrip_grids::grid1_corners, scrip_grids::grid1_frac, scrip_grids::grid1_mask, scrip_grids::grid2_area, scrip_grids::grid2_center_lat, scrip_grids::grid2_center_lon, scrip_grids::grid2_centroid_lat, scrip_grids::grid2_centroid_lon, scrip_grids::grid2_corner_lat, scrip_grids::grid2_corner_lon, scrip_grids::grid2_corners, scrip_grids::grid2_frac, scrip_grids::grid2_mask, intersection(), line_integral(), locate_point(), modify_polar_cell(), scrip_grids::north_thresh, scrip_grids::npseg, scrip_remap_vars::num_wts, scrip_constants::pi, scrip_constants::pi2, ptincell(), scrip_grids::south_thresh, scrip_grids::special_polar_cell1, scrip_grids::special_polar_cell2, store_link_cnsrv(), and scrip_constants::tiny.

Referenced by cellblock_integrate(), and remap_conserv().

◆ cellblock_integrate()

subroutine scrip_remap_conservative::cellblock_integrate ( integer (scrip_i4)  ibegin,
integer (scrip_i4)  iend,
integer (scrip_i4)  grid_num,
integer (scrip_i4)  phi_or_theta 
)

Definition at line 982 of file scrip_remap_conservative.f.

982 
983  integer (SCRIP_i4) :: ibegin, iend, grid_num, phi_or_theta
984 
985  integer (SCRIP_i4) :: cell_add
986 
987 
988  do cell_add = ibegin, iend
989 
990  call cell_integrate(cell_add, grid_num, phi_or_theta)
991 
992  enddo
993 
994 
995 

References cell_integrate().

◆ converge_to_bdry()

subroutine scrip_remap_conservative::converge_to_bdry ( integer (scrip_i4), intent(in)  cell_add,
integer (scrip_i4), intent(in)  cell_grid_num,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_lat,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_lon,
real (scrip_r8), intent(in)  cell_center_lat,
real (scrip_r8), intent(in)  cell_center_lon,
real (scrip_r8), intent(in)  inpt_lat,
real (scrip_r8), intent(in)  inpt_lon,
real (scrip_r8), intent(in)  outpt_lat,
real (scrip_r8), intent(in)  outpt_lon,
real (scrip_r8), intent(out)  bpt_lat,
real (scrip_r8), intent(out)  bpt_lon,
integer (scrip_i4), intent(out)  bedgeid 
)

Definition at line 6217 of file scrip_remap_conservative.f.

6217 
6218 !----------------------------------------------------------------------
6219 !
6220 ! Input variables
6221 !
6222 !----------------------------------------------------------------------
6223 
6224  integer (SCRIP_i4), intent(in) ::
6225  & cell_add, ! Cell in which we are operating
6226  & cell_grid_num, ! Grid to which cell belongs
6227  & ncorners ! Number of corners in cell
6228 
6229  real (SCRIP_r8), dimension(ncorners), intent(in) ::
6230  & cell_corner_lat, ! Latitude values of cell corners
6231  & cell_corner_lon ! Longitude values of cell corners
6232 
6233  real (SCRIP_r8), intent(in) ::
6234  & cell_center_lat, ! Latitude of cell center
6235  & cell_center_lon, ! Longitude of cell center,
6236  & inpt_lat, ! Latitude of inside point
6237  & inpt_lon, ! Longitude of inside point
6238  & outpt_lat, ! Latitude of outside point
6239  & outpt_lon ! Longitude of outside point
6240 
6241 
6242 !----------------------------------------------------------------------
6243 !
6244 ! Output variables
6245 !
6246 !----------------------------------------------------------------------
6247 
6248  real (SCRIP_r8), intent(out) ::
6249  & bpt_lat, ! Latitude of boundary point
6250  & bpt_lon ! Longitude of boundary point
6251 
6252  integer (SCRIP_i4), intent(out) ::
6253  & bedgeid ! ID of edge that point converged to
6254 
6255 !----------------------------------------------------------------------
6256 !
6257 ! Local variables
6258 !
6259 !----------------------------------------------------------------------
6260 
6261  logical (SCRIP_logical) ::
6262  & converged,
6263  & lboundary,
6264  & inpoly
6265 
6266  integer (SCRIP_i4) ::
6267  & it
6268 
6269  real (SCRIP_r8) ::
6270  & lat1, lon1,
6271  & lat2, lon2,
6272  & midlat, midlon
6273 
6274  bedgeid = 0
6275 
6276  lat1 = inpt_lat
6277  lon1 = inpt_lon
6278  lat2 = outpt_lat
6279  lon2 = outpt_lon
6280 
6281 
6282  converged = .false.
6283  it = 0
6284  do while (.not. converged)
6285 
6286  midlat = (lat1+lat2)/2.0
6287  if (abs(lon1-lon2) < pi) then
6288  midlon = (lon1+lon2)/2.0
6289  else
6290  midlon = (lon1+lon2)/2.0 - pi2
6291  endif
6292 
6293 
6294  call ptincell(midlat, midlon,
6295  & cell_add, ncorners,
6296  & cell_corner_lat, cell_corner_lon,
6297  & cell_center_lat, cell_center_lon,
6298  & cell_grid_num,
6299  & inpoly, lboundary, bedgeid)
6300 
6301  if (inpoly) then
6302  lat1 = midlat
6303  lon1 = midlon
6304  else
6305  lat2 = midlat
6306  lon2 = midlon
6307  endif
6308 
6309  if (abs(lat1-lat2) < tiny .and.
6310  & abs(lon1-lon2) < tiny .and. lboundary) then
6311  converged = .true.
6312  endif
6313 
6314  if (it > 100) then
6315  exit
6316  endif
6317 
6318  it = it + 1
6319  enddo ! do while (not converged)
6320 
6321 
6322  bpt_lat = midlat
6323  bpt_lon = midlon
6324 

References scrip_constants::pi, scrip_constants::pi2, ptincell(), and scrip_constants::tiny.

Referenced by cell_integrate().

◆ find_adj_cell()

subroutine scrip_remap_conservative::find_adj_cell ( integer (scrip_i4), intent(in)  cell_add,
integer (scrip_i4), intent(in)  edge_id,
integer (scrip_i4), intent(in)  cell_grid_num,
integer (scrip_i4), intent(out)  adj_add 
)

Definition at line 5987 of file scrip_remap_conservative.f.

5987 
5988 !----------------------------------------------------------------------
5989 !
5990 ! Input variables
5991 !
5992 !----------------------------------------------------------------------
5993 
5994  integer (SCRIP_i4), intent(in) ::
5995  & cell_add, ! cell whose edge we are checking
5996  & edge_id, ! index of edge that we are check
5997  & cell_grid_num ! grid to which cell belongs
5998 
5999 !----------------------------------------------------------------------
6000 !
6001 ! Output variables
6002 !
6003 !----------------------------------------------------------------------
6004 
6005  integer (SCRIP_i4), intent(out) :: adj_add
6006 
6007 !----------------------------------------------------------------------
6008 !
6009 ! Local variables
6010 !
6011 !----------------------------------------------------------------------
6012 
6013  integer (SCRIP_i4) :: i, inx, n, global_add
6014  logical (SCRIP_logical) :: found
6015  real (SCRIP_r8) :: lat1, lon1, lat2, lon2
6016 
6017  adj_add = 0
6018 
6019  if (cell_grid_num .eq. 1) then
6020 
6021  i = edge_id
6022  inx = 1 + mod(edge_id,grid1_corners)
6023 
6024  lat1 = grid1_corner_lat(i,cell_add)
6025  lon1 = grid1_corner_lon(i,cell_add)
6026  lat2 = grid1_corner_lat(inx,cell_add)
6027  lon2 = grid1_corner_lon(inx,cell_add)
6028 
6029  !***
6030  !*** Often the cell with the next or previous index is
6031  !*** the adjacent cell. Check that first
6032  !***
6033 
6034  if (cell_add .lt. grid1_size) then
6035 
6036  global_add = cell_add + 1
6037 
6038  do i = 1, grid1_corners
6039  inx = mod(i,grid1_corners)+1
6040  if (abs(grid1_corner_lat(inx,global_add)-lat1) .le. tiny
6041  & .and.
6042  & abs(grid1_corner_lat(i,global_add)-lat2) .le. tiny
6043  & .and.
6044  & abs(grid1_corner_lon(inx,global_add)-lon1) .le. tiny
6045  & .and.
6046  & abs(grid1_corner_lon(i,global_add)-lon2) .le. tiny)
6047  & then
6048 
6049  adj_add = global_add
6050  return
6051  endif
6052  enddo
6053 
6054  endif
6055 
6056  if (cell_add .gt. 1) then
6057 
6058  global_add = cell_add - 1
6059 
6060  do i = 1, grid1_corners
6061  inx = mod(i,grid1_corners)+1
6062  if (abs(grid1_corner_lat(inx,global_add)-lat1) .le. tiny
6063  & .and.
6064  & abs(grid1_corner_lat(i,global_add)-lat2) .le. tiny
6065  & .and.
6066  & abs(grid1_corner_lon(inx,global_add)-lon1) .le. tiny
6067  & .and.
6068  & abs(grid1_corner_lon(i,global_add)-lon2) .le. tiny)
6069  & then
6070 
6071  adj_add = global_add
6072  return
6073  endif
6074  enddo
6075 
6076  endif
6077 
6078 
6079 
6080  else
6081 
6082  i = edge_id
6083  inx = 1 + mod(edge_id,grid2_corners)
6084 
6085  lat1 = grid2_corner_lat(i,cell_add)
6086  lon1 = grid2_corner_lon(i,cell_add)
6087  lat2 = grid2_corner_lat(inx,cell_add)
6088  lon2 = grid2_corner_lon(inx,cell_add)
6089 
6090 
6091  !***
6092  !*** Often the cell with the next or previous index is
6093  !*** the adjacent cell. Check that first
6094  !***
6095 
6096  if (cell_add .lt. grid2_size) then
6097 
6098  global_add = cell_add + 1
6099 
6100  do i = 1, grid2_corners
6101  inx = mod(i,grid2_corners)+1
6102  if (abs(grid2_corner_lat(inx,global_add)-lat1) .le. tiny
6103  & .and.
6104  & abs(grid2_corner_lat(i,global_add)-lat2) .le. tiny
6105  & .and.
6106  & abs(grid2_corner_lon(inx,global_add)-lon1) .le. tiny
6107  & .and.
6108  & abs(grid2_corner_lon(i,global_add)-lon2) .le. tiny)
6109  & then
6110 
6111  adj_add = global_add
6112  return
6113  endif
6114  enddo
6115 
6116  endif
6117 
6118  if (cell_add .gt. 1) then
6119 
6120  global_add = cell_add - 1
6121 
6122  do i = 1, grid2_corners
6123  inx = mod(i,grid2_corners)+1
6124  if (abs(grid2_corner_lat(inx,global_add)-lat1) .le. tiny
6125  & .and.
6126  & abs(grid2_corner_lat(i,global_add)-lat2) .le. tiny
6127  & .and.
6128  & abs(grid2_corner_lon(inx,global_add)-lon1) .le. tiny
6129  & .and.
6130  & abs(grid2_corner_lon(i,global_add)-lon2) .le. tiny)
6131  & then
6132 
6133  adj_add = global_add
6134  return
6135  endif
6136  enddo
6137 
6138  endif
6139 
6140 
6141  endif
6142 
6143 
6144  if (cell_add /= last_cell_find_adj_cell .or.
6145  & cell_grid_num /= last_cell_grid_num_find_adj_cell) then
6146 
6147  last_cell_find_adj_cell = cell_add
6148  last_cell_grid_num_find_adj_cell = cell_grid_num
6149 
6150  if (first_call_find_adj_cell) then
6151  first_call_find_adj_cell = .false.
6152  last_cell_find_adj_cell = 0
6153  last_cell_grid_num_find_adj_cell = 0
6154  else
6155  if (num_srch_cells_find_adj_cell .gt. 0) then
6156  deallocate(srch_add_find_adj_cell,
6157  & srch_corner_lat_find_adj_cell,
6158  & srch_corner_lon_find_adj_cell,
6159  & srch_center_lat_find_adj_cell,
6160  & srch_center_lon_find_adj_cell)
6161  endif
6162  endif
6163 
6164  call get_srch_cells(cell_add, cell_grid_num, cell_grid_num,
6165  & num_srch_cells_find_adj_cell, srch_add_find_adj_cell,
6166  & srch_corners_find_adj_cell, srch_corner_lat_find_adj_cell,
6167  & srch_corner_lon_find_adj_cell,
6168  & srch_center_lat_find_adj_cell,
6169  & srch_center_lon_find_adj_cell)
6170 
6171  endif
6172 
6173 
6174  found = .false.
6175  do n = 1, num_srch_cells_find_adj_cell
6176 
6177  global_add = srch_add_find_adj_cell(n)
6178 
6179  do i = 1, srch_corners_find_adj_cell
6180  inx = mod(i,srch_corners_find_adj_cell)+1
6181  if (abs(srch_corner_lat_find_adj_cell(inx,n)-lat1) .le. tiny
6182  & .and.
6183  & abs(srch_corner_lat_find_adj_cell(i,n)-lat2) .le. tiny
6184  & .and.
6185  & abs(srch_corner_lon_find_adj_cell(inx,n)-lon1) .le.tiny
6186  & .and.
6187  & abs(srch_corner_lon_find_adj_cell(i,n)-lon2) .le. tiny)
6188  & then
6189 
6190  adj_add = global_add
6191  found = .true.
6192 
6193  exit
6194  endif
6195  enddo
6196 
6197  if (found) exit
6198 
6199  enddo
6200 
6201  return

References first_call_find_adj_cell, get_srch_cells(), scrip_grids::grid1_corner_lat, scrip_grids::grid1_corner_lon, scrip_grids::grid1_corners, scrip_grids::grid1_size, scrip_grids::grid2_corner_lat, scrip_grids::grid2_corner_lon, scrip_grids::grid2_corners, scrip_grids::grid2_size, last_cell_find_adj_cell, last_cell_grid_num_find_adj_cell, num_srch_cells_find_adj_cell, srch_add_find_adj_cell, srch_center_lat_find_adj_cell, srch_center_lon_find_adj_cell, srch_corner_lat_find_adj_cell, srch_corner_lon_find_adj_cell, srch_corners_find_adj_cell, and scrip_constants::tiny.

Referenced by cell_integrate().

◆ get_srch_cells()

subroutine scrip_remap_conservative::get_srch_cells ( integer (scrip_i4), intent(in)  cell_add,
integer (scrip_i4), intent(in)  cell_grid_num,
integer (scrip_i4), intent(in)  srch_grid_num,
integer (scrip_i4), intent(out)  num_srch_cells,
integer (scrip_i4), dimension(:), intent(out), allocatable  srch_add,
integer (scrip_i4), intent(out)  srch_corners,
real (scrip_r8), dimension(:,:), intent(out), allocatable  srch_corner_lat,
real (scrip_r8), dimension(:,:), intent(out), allocatable  srch_corner_lon,
real (scrip_r8), dimension(:), intent(out), allocatable  srch_center_lat,
real (scrip_r8), dimension(:), intent(out), allocatable  srch_center_lon 
)

Definition at line 5557 of file scrip_remap_conservative.f.

5557 
5558 !----------------------------------------------------------------------
5559 !
5560 ! Input arguments
5561 !
5562 !----------------------------------------------------------------------
5563 
5564  integer (SCRIP_i4), intent(in) ::
5565  & cell_add, ! cell in whose nbrhood we must find other
5566  & cell_grid_num, ! cells grid number from which 'cell_add'
5567  & srch_grid_num ! is grid number in which we must find
5568  ! search cells
5569 
5570 !----------------------------------------------------------------------
5571 !
5572 ! Output arguments
5573 !
5574 !----------------------------------------------------------------------
5575 
5576  integer (SCRIP_i4), intent(out) ::
5577  & num_srch_cells,
5578  & srch_corners ! Number of corners for search cells
5579 
5580  integer (SCRIP_i4), dimension(:), allocatable, intent(out) ::
5581  & srch_add ! Global addresses of search cells
5582 
5583  real (SCRIP_r8), dimension(:,:), allocatable, intent(out) ::
5584  & srch_corner_lat, srch_corner_lon
5585 
5586  real (SCRIP_r8), dimension(:), allocatable, intent(out) ::
5587  & srch_center_lat, srch_center_lon
5588 
5589 
5590 !-----------------------------------------------------------------------
5591 !
5592 ! Local arguments
5593 !
5594 !-----------------------------------------------------------------------
5595 
5596  logical (SCRIP_logical), dimension(:), allocatable ::
5597  & srch_mask
5598 
5599  integer (SCRIP_i4) :: grid1_add, grid2_add, max_add, min_add,
5600  & n
5601 
5602 !-----------------------------------------------------------------------
5603 
5604  num_srch_cells = 0
5605 
5606  !***
5607  !*** restrict searches first using search bins
5608  !***
5609 
5610  if (last_cell_add_get_srch_cells /= cell_add .or.
5611  & last_cell_grid_num_get_srch_cells /= cell_grid_num .or.
5612  & last_srch_grid_num_get_srch_cells /= srch_grid_num) then
5613 
5614  if (first_call_get_srch_cells) then
5615  first_call_get_srch_cells = .false.
5616  num_srch_cells_loc_get_srch_cells = 0
5617  srch_corners_loc_get_srch_cells = 0
5618  last_cell_add_get_srch_cells = 0
5619  last_cell_grid_num_get_srch_cells = 0
5620  last_srch_grid_num_get_srch_cells = 0
5621  else
5622  if (num_srch_cells_loc_get_srch_cells .gt. 0) then
5623  deallocate(srch_add_loc_get_srch_cells,
5624  & srch_corner_lat_loc_get_srch_cells,
5625  & srch_corner_lon_loc_get_srch_cells,
5626  & srch_center_lat_loc_get_srch_cells,
5627  & srch_center_lon_loc_get_srch_cells)
5628  endif
5629 
5630  endif
5631 
5632 
5633  last_cell_add_get_srch_cells = cell_add
5634  last_cell_grid_num_get_srch_cells = cell_grid_num
5635  last_srch_grid_num_get_srch_cells = srch_grid_num
5636 
5637 
5638  if (cell_grid_num == 1) then
5639 
5640  if (srch_grid_num == 1) then
5641 
5642  !*** Grid 1 neighbors of grid 1 cell
5643 
5644  allocate(srch_mask(grid1_size))
5645 
5646  min_add = grid1_size
5647  max_add = 1
5648  do n=1,num_srch_bins
5649  if (cell_add >= bin_addr1(1,n) .and.
5650  & cell_add <= bin_addr1(2,n)) then
5651  min_add = min(min_add, bin_addr1(1,n))
5652  max_add = max(max_add, bin_addr1(2,n))
5653  endif
5654  end do
5655 
5656  !***
5657  !*** further restrict searches using bounding boxes
5658  !***
5659 
5660  num_srch_cells_loc_get_srch_cells = 0
5661  do grid1_add = min_add,max_add
5662  srch_mask(grid1_add) =
5663  & (grid1_bound_box(1,grid1_add) <=
5664  & grid1_bound_box(2,cell_add)) .and.
5665  & (grid1_bound_box(2,grid1_add) >=
5666  & grid1_bound_box(1,cell_add)) .and.
5667  & (grid1_bound_box(3,grid1_add) <=
5668  & grid1_bound_box(4,cell_add)) .and.
5669  & (grid1_bound_box(4,grid1_add) >=
5670  & grid1_bound_box(3,cell_add))
5671 
5672  if (srch_mask(grid1_add))
5673  & num_srch_cells_loc_get_srch_cells =
5674  & num_srch_cells_loc_get_srch_cells+1
5675  end do
5676 
5677  if (num_srch_cells_loc_get_srch_cells /= 0) then
5678 
5679  !***
5680  !*** create search arrays
5681  !***
5682 
5683  allocate(srch_add_loc_get_srch_cells
5684  & (num_srch_cells_loc_get_srch_cells),
5685  & srch_corner_lat_loc_get_srch_cells
5686  & (grid1_corners,num_srch_cells_loc_get_srch_cells),
5687  & srch_corner_lon_loc_get_srch_cells
5688  & (grid1_corners,num_srch_cells_loc_get_srch_cells),
5689  & srch_center_lat_loc_get_srch_cells
5690  & (num_srch_cells_loc_get_srch_cells),
5691  & srch_center_lon_loc_get_srch_cells
5692  & (num_srch_cells_loc_get_srch_cells))
5693 
5694  n = 0
5695  do grid1_add = min_add,max_add
5696  if (srch_mask(grid1_add)) then
5697  n = n+1
5698  srch_add_loc_get_srch_cells(n) = grid1_add
5699  srch_corner_lat_loc_get_srch_cells(:,n) =
5700  & grid1_corner_lat(:,grid1_add)
5701  srch_corner_lon_loc_get_srch_cells(:,n) =
5702  & grid1_corner_lon(:,grid1_add)
5703  srch_center_lat_loc_get_srch_cells(n) =
5704  & grid1_center_lat(grid1_add)
5705  srch_center_lon_loc_get_srch_cells(n) =
5706  & grid1_center_lon(grid1_add)
5707  endif
5708  end do
5709 
5710  srch_corners_loc_get_srch_cells = grid1_corners
5711  endif
5712 
5713  deallocate(srch_mask)
5714 
5715  else
5716 
5717  !*** Grid 2 neighbors of grid 1 cell
5718 
5719  allocate(srch_mask(grid2_size))
5720 
5721  min_add = grid2_size
5722  max_add = 1
5723  do n=1,num_srch_bins
5724  if (cell_add >= bin_addr1(1,n) .and.
5725  & cell_add <= bin_addr1(2,n)) then
5726  min_add = min(min_add, bin_addr2(1,n))
5727  max_add = max(max_add, bin_addr2(2,n))
5728  endif
5729  end do
5730 
5731  !***
5732  !*** further restrict searches using bounding boxes
5733  !***
5734 
5735  num_srch_cells_loc_get_srch_cells = 0
5736  do grid2_add = min_add,max_add
5737  srch_mask(grid2_add) =
5738  & (grid2_bound_box(1,grid2_add) <=
5739  & grid1_bound_box(2,cell_add)) .and.
5740  & (grid2_bound_box(2,grid2_add) >=
5741  & grid1_bound_box(1,cell_add)) .and.
5742  & (grid2_bound_box(3,grid2_add) <=
5743  & grid1_bound_box(4,cell_add)) .and.
5744  & (grid2_bound_box(4,grid2_add) >=
5745  & grid1_bound_box(3,cell_add))
5746 
5747 
5748 
5749  if (srch_mask(grid2_add))
5750  & num_srch_cells_loc_get_srch_cells =
5751  & num_srch_cells_loc_get_srch_cells+1
5752  end do
5753 
5754 
5755  if (num_srch_cells_loc_get_srch_cells /= 0) then
5756 
5757  !***
5758  !*** create search arrays
5759  !***
5760 
5761  allocate(srch_add_loc_get_srch_cells(
5762  & num_srch_cells_loc_get_srch_cells),
5763  & srch_corner_lat_loc_get_srch_cells(
5764  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5765  & srch_corner_lon_loc_get_srch_cells(
5766  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5767  & srch_center_lat_loc_get_srch_cells(
5768  & num_srch_cells_loc_get_srch_cells),
5769  & srch_center_lon_loc_get_srch_cells(
5770  & num_srch_cells_loc_get_srch_cells))
5771 
5772  n = 0
5773  do grid2_add = min_add,max_add
5774  if (srch_mask(grid2_add)) then
5775  n = n+1
5776  srch_add_loc_get_srch_cells(n) = grid2_add
5777  srch_corner_lat_loc_get_srch_cells(:,n) =
5778  & grid2_corner_lat(:,grid2_add)
5779  srch_corner_lon_loc_get_srch_cells(:,n) =
5780  & grid2_corner_lon(:,grid2_add)
5781  srch_center_lat_loc_get_srch_cells(n) =
5782  & grid2_center_lat(grid2_add)
5783  srch_center_lon_loc_get_srch_cells(n) =
5784  & grid2_center_lon(grid2_add)
5785  endif
5786  end do
5787 
5788  srch_corners_loc_get_srch_cells = grid2_corners
5789  endif
5790 
5791  deallocate(srch_mask)
5792  endif
5793 
5794  else
5795 
5796  if (srch_grid_num == 1) then
5797 
5798  !*** Grid 1 neighbors of grid 2 cell
5799 
5800  allocate(srch_mask(grid1_size))
5801 
5802  min_add = grid1_size
5803  max_add = 1
5804  do n=1,num_srch_bins
5805  if (cell_add >= bin_addr2(1,n) .and.
5806  & cell_add <= bin_addr2(2,n)) then
5807  min_add = min(min_add, bin_addr1(1,n))
5808  max_add = max(max_add, bin_addr1(2,n))
5809  endif
5810  end do
5811 
5812  !***
5813  !*** further restrict searches using bounding boxes
5814  !***
5815 
5816  num_srch_cells_loc_get_srch_cells = 0
5817  do grid1_add = min_add,max_add
5818  srch_mask(grid1_add) =
5819  & (grid1_bound_box(1,grid1_add) <=
5820  & grid2_bound_box(2,cell_add)) .and.
5821  & (grid1_bound_box(2,grid1_add) >=
5822  & grid2_bound_box(1,cell_add)) .and.
5823  & (grid1_bound_box(3,grid1_add) <=
5824  & grid2_bound_box(4,cell_add)) .and.
5825  & (grid1_bound_box(4,grid1_add) >=
5826  & grid2_bound_box(3,cell_add))
5827 
5828  if (srch_mask(grid1_add))
5829  & num_srch_cells_loc_get_srch_cells =
5830  & num_srch_cells_loc_get_srch_cells+1
5831  end do
5832 
5833 
5834  if (num_srch_cells_loc_get_srch_cells /= 0) then
5835 
5836  !***
5837  !*** create search arrays
5838  !***
5839 
5840  allocate(srch_add_loc_get_srch_cells(
5841  & num_srch_cells_loc_get_srch_cells),
5842  & srch_corner_lat_loc_get_srch_cells(
5843  & grid1_corners,num_srch_cells_loc_get_srch_cells),
5844  & srch_corner_lon_loc_get_srch_cells(
5845  & grid1_corners,num_srch_cells_loc_get_srch_cells),
5846  & srch_center_lat_loc_get_srch_cells(
5847  & num_srch_cells_loc_get_srch_cells),
5848  & srch_center_lon_loc_get_srch_cells(
5849  & num_srch_cells_loc_get_srch_cells))
5850 
5851  n = 0
5852  do grid1_add = min_add,max_add
5853  if (srch_mask(grid1_add)) then
5854  n = n+1
5855  srch_add_loc_get_srch_cells(n) = grid1_add
5856  srch_corner_lat_loc_get_srch_cells(:,n) =
5857  & grid1_corner_lat(:,grid1_add)
5858  srch_corner_lon_loc_get_srch_cells(:,n) =
5859  & grid1_corner_lon(:,grid1_add)
5860  srch_center_lat_loc_get_srch_cells(n) =
5861  & grid1_center_lat(grid1_add)
5862  srch_center_lon_loc_get_srch_cells(n) =
5863  & grid1_center_lon(grid1_add)
5864  endif
5865  end do
5866 
5867  srch_corners_loc_get_srch_cells = grid1_corners
5868  endif
5869 
5870  deallocate(srch_mask)
5871 
5872  else
5873 
5874  !*** Grid 2 neighbors of grid 2 cell
5875 
5876  allocate(srch_mask(grid2_size))
5877 
5878  min_add = grid2_size
5879  max_add = 1
5880  do n=1,num_srch_bins
5881  if (cell_add >= bin_addr2(1,n) .and.
5882  & cell_add <= bin_addr2(2,n)) then
5883  min_add = min(min_add, bin_addr2(1,n))
5884  max_add = max(max_add, bin_addr2(2,n))
5885  endif
5886  end do
5887 
5888  !***
5889  !*** further restrict searches using bounding boxes
5890  !***
5891 
5892  num_srch_cells_loc_get_srch_cells = 0
5893  do grid2_add = min_add,max_add
5894  srch_mask(grid2_add) =
5895  & (grid2_bound_box(1,grid2_add) <=
5896  & grid2_bound_box(2,cell_add)) .and.
5897  & (grid2_bound_box(2,grid2_add) >=
5898  & grid2_bound_box(1,cell_add)) .and.
5899  & (grid2_bound_box(3,grid2_add) <=
5900  & grid2_bound_box(4,cell_add)) .and.
5901  & (grid2_bound_box(4,grid2_add) >=
5902  & grid2_bound_box(3,cell_add))
5903 
5904  if (srch_mask(grid2_add))
5905  & num_srch_cells_loc_get_srch_cells =
5906  & num_srch_cells_loc_get_srch_cells+1
5907  end do
5908 
5909 
5910  if (num_srch_cells_loc_get_srch_cells /= 0) then
5911 
5912  !***
5913  !*** create search arrays
5914  !***
5915 
5916  allocate(srch_add_loc_get_srch_cells(
5917  & num_srch_cells_loc_get_srch_cells),
5918  & srch_corner_lat_loc_get_srch_cells(
5919  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5920  & srch_corner_lon_loc_get_srch_cells(
5921  & grid2_corners,num_srch_cells_loc_get_srch_cells),
5922  & srch_center_lat_loc_get_srch_cells(
5923  & num_srch_cells_loc_get_srch_cells),
5924  & srch_center_lon_loc_get_srch_cells(
5925  & num_srch_cells_loc_get_srch_cells))
5926 
5927  n = 0
5928  do grid2_add = min_add,max_add
5929  if (srch_mask(grid2_add)) then
5930  n = n+1
5931  srch_add_loc_get_srch_cells(n) = grid2_add
5932  srch_corner_lat_loc_get_srch_cells(:,n) =
5933  & grid2_corner_lat(:,grid2_add)
5934  srch_corner_lon_loc_get_srch_cells(:,n) =
5935  & grid2_corner_lon(:,grid2_add)
5936  srch_center_lat_loc_get_srch_cells(n) =
5937  & grid2_center_lat(grid2_add)
5938  srch_center_lon_loc_get_srch_cells(n) =
5939  & grid2_center_lon(grid2_add)
5940  endif
5941  end do
5942 
5943  srch_corners_loc_get_srch_cells = grid2_corners
5944  endif
5945 
5946  deallocate(srch_mask)
5947 
5948  endif
5949 
5950  endif
5951 
5952  endif
5953 
5954 
5955  num_srch_cells = num_srch_cells_loc_get_srch_cells
5956 
5957  if (num_srch_cells .eq. 0) then
5958  return
5959  endif
5960 
5961  srch_corners = srch_corners_loc_get_srch_cells
5962  allocate(srch_add(num_srch_cells),
5963  & srch_corner_lat(srch_corners,num_srch_cells),
5964  & srch_corner_lon(srch_corners,num_srch_cells),
5965  & srch_center_lat(num_srch_cells),
5966  & srch_center_lon(num_srch_cells))
5967  srch_add = srch_add_loc_get_srch_cells
5968  srch_corner_lat = srch_corner_lat_loc_get_srch_cells
5969  srch_corner_lon = srch_corner_lon_loc_get_srch_cells
5970  srch_center_lat = srch_center_lat_loc_get_srch_cells
5971  srch_center_lon = srch_center_lon_loc_get_srch_cells
5972 

References scrip_grids::bin_addr1, scrip_grids::bin_addr2, first_call_get_srch_cells, scrip_grids::grid1_bound_box, scrip_grids::grid1_center_lat, scrip_grids::grid1_center_lon, scrip_grids::grid1_corner_lat, scrip_grids::grid1_corner_lon, scrip_grids::grid1_corners, scrip_grids::grid1_size, scrip_grids::grid2_bound_box, scrip_grids::grid2_center_lat, scrip_grids::grid2_center_lon, scrip_grids::grid2_corner_lat, scrip_grids::grid2_corner_lon, scrip_grids::grid2_corners, scrip_grids::grid2_size, last_cell_add_get_srch_cells, last_cell_grid_num_get_srch_cells, last_srch_grid_num_get_srch_cells, scrip_grids::num_srch_bins, num_srch_cells_loc_get_srch_cells, srch_add_loc_get_srch_cells, srch_center_lat_loc_get_srch_cells, srch_center_lon_loc_get_srch_cells, srch_corner_lat_loc_get_srch_cells, srch_corner_lon_loc_get_srch_cells, and srch_corners_loc_get_srch_cells.

Referenced by find_adj_cell(), locate_point(), and locate_segstart().

◆ intersection()

subroutine scrip_remap_conservative::intersection ( integer (scrip_i4), intent(in)  seg_cell_id,
integer (scrip_i4), intent(in)  seg_grid_id,
real (scrip_r8), intent(in)  beglat,
real (scrip_r8), intent(in)  beglon,
real (scrip_r8), intent(in)  endlat,
real (scrip_r8), intent(in)  endlon,
real (scrip_r8), dimension(2), intent(inout)  begseg,
integer (scrip_i4), intent(in)  begedge,
integer (scrip_i4), intent(in)  cell_id,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_lat,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_lon,
integer (scrip_i4), intent(in)  cell_grid_id,
real (scrip_r8), intent(out)  intrsct_lat,
real (scrip_r8), intent(out)  intrsct_lon,
integer (scrip_i4), intent(out)  intedge,
real (scrip_r8), intent(out)  sinang2,
logical (scrip_logical), intent(out)  lcoinc,
logical (scrip_logical), intent(out)  lthresh 
)

Definition at line 2196 of file scrip_remap_conservative.f.

2196 
2197 !-----------------------------------------------------------------------
2198 !
2199 ! this routine finds the intersection of a line segment given by
2200 ! beglon, endlon, etc. with a cell from another grid
2201 ! A coincidence flag is returned if the segment is entirely
2202 ! coincident with an edge of the opposite.
2203 !
2204 !-----------------------------------------------------------------------
2205 
2206 !-----------------------------------------------------------------------
2207 !
2208 ! intent(in):
2209 !
2210 !-----------------------------------------------------------------------
2211 
2212  integer (SCRIP_i4), intent(in) ::
2213  & seg_cell_id ! ID of cell that intersecting segment is from
2214 
2215  integer (SCRIP_i4), intent(in) ::
2216  & seg_grid_id ! ID of grid that intersecting segment is from
2217 
2218  real (SCRIP_r8), intent(in) ::
2219  & beglat, beglon,! beginning lat/lon endpoints for segment
2220  & endlat, endlon ! ending lat/lon endpoints for segment
2221 
2222  real (SCRIP_r8), dimension(2), intent(inout) ::
2223  & begseg ! begin lat/lon of full segment
2224 
2225  integer (SCRIP_i4), intent(in) ::
2226  & begedge ! edge that beginning point is on (can be 0)
2227 
2228  integer (SCRIP_i4), intent(in) ::
2229  & cell_id ! cell to intersect with
2230 
2231  integer (SCRIP_i4), intent(in) ::
2232  & ncorners ! number of corners of cell
2233 
2234  real (SCRIP_r8), dimension(ncorners), intent(in) ::
2235  & cell_corner_lat, ! coordinates of cell corners
2236  & cell_corner_lon
2237 
2238  integer (SCRIP_i4), intent(in) ::
2239  & cell_grid_id ! which grid is the cell from?
2240 
2241 
2242 !-----------------------------------------------------------------------
2243 !
2244 ! intent(out):
2245 !
2246 !-----------------------------------------------------------------------
2247 
2248  real (SCRIP_r8), intent(out) ::
2249  & intrsct_lat,
2250  & intrsct_lon ! lat/lon coords of intersection
2251 
2252  real (SCRIP_r8), intent(out) ::
2253  & sinang2 ! square of sine of angle between
2254  ! intersecting lines
2255 
2256  integer (SCRIP_i4), intent(out) ::
2257  & intedge ! edge that is intersected
2258 
2259  logical (SCRIP_logical), intent(out) ::
2260  & lcoinc ! True if segment is coincident with
2261  ! a cell edge
2262 
2263  logical (SCRIP_logical), intent(out) ::
2264  & lthresh ! True if segment crosses threshold
2265 
2266 !-----------------------------------------------------------------------
2267 !
2268 ! local variables
2269 !
2270 !-----------------------------------------------------------------------
2271 
2272  integer (SCRIP_i4) ::
2273  & n, next_n
2274 
2275  logical (SCRIP_logical) ::
2276  & found, first
2277 
2278  real (SCRIP_r8) ::
2279  & lon1, lon2, ! local longitude variables for segment
2280  & lat1, lat2, ! local latitude variables for segment
2281  & grdlon1, grdlon2, ! local longitude variables for grid cell
2282  & grdlat1, grdlat2, ! local latitude variables for grid cell
2283  & vec1_lat, vec1_lon,
2284  & vec2_lat, vec2_lon, !
2285  & vec3_lat, vec3_lon, ! vectors and vector products used
2286  & cross_product, ! during grid search
2287  & dot_product, !
2288  & lensqr1, lensqr2, !
2289  & lensqr3, !
2290  & s1, s2, determ,
2291  & mat1, mat2, ! variables used for linear solve to
2292  & mat3, mat4, ! find intersection
2293  & rhs1, rhs2, !
2294  & denom,
2295  & begsegloc(2), ! local copy of full segment start
2296  & dist2, ! distance from start pt to intersection
2297  ! pt
2298  & maxdist2, ! max dist from start pt to any
2299  ! intersection pt
2300  & max_intrsct_lat, ! latitude of farthest intersection point
2301  & max_intrsct_lon, ! longitude of farthest intersection
2302  ! point
2303  & minlat, maxlat, ! min and max latitudes of segment
2304  & minlon, maxlon, ! min and max longitudes of segment
2305  & tmplat, tmplon
2306 
2307 
2308 !-----------------------------------------------------------------------
2309 !
2310 ! initialize defaults, flags, etc.
2311 !
2312 !-----------------------------------------------------------------------
2313 
2314  lcoinc = .false.
2315  lthresh = .false.
2316  intedge = 0
2317  first = .true.
2318 
2319  lat1 = beglat
2320  lon1 = beglon
2321  lat2 = endlat
2322  lon2 = endlon
2323 
2324  ! No edge is allowed to span more than pi radians
2325  ! Accordingly transform one or the other end point
2326 
2327  if ((lon2-lon1) > pi) then
2328  lon2 = lon2 - pi2
2329  else if ((lon2-lon1) < -pi) then
2330  lon1 = lon1 - pi2
2331  endif
2332  s1 = zero
2333 
2334 !-----------------------------------------------------------------------
2335 !
2336 ! loop over sides of the cell to find intersection with side
2337 ! must check all sides for coincidences or intersections
2338 !
2339 !-----------------------------------------------------------------------
2340 
2341  if (beglat > north_thresh .or. beglat < south_thresh) then
2342 
2343  !*** Special intersection routine for cells near the pole
2344  !*** Intersection is done in a transformed space using
2345  !*** multi-segmented representation of the cell
2346 
2347  call pole_intersection(cell_id,ncorners,
2348  & cell_corner_lat,cell_corner_lon,cell_grid_id,
2349  & beglat, beglon, endlat,
2350  & endlon, begseg, begedge,
2351  & intedge,intrsct_lat,intrsct_lon,
2352  & sinang2,lcoinc,lthresh)
2353 
2354  return
2355 
2356  endif
2357 
2358 
2359  maxdist2 = -9999999.0
2360 
2361  begsegloc(1) = begseg(1)
2362  begsegloc(2) = begseg(2)
2363 
2364  lthresh = .false.
2365  intrsct_loop: do n=1,ncorners
2366  next_n = mod(n,ncorners) + 1
2367 
2368  grdlat1 = cell_corner_lat(n)
2369  grdlon1 = cell_corner_lon(n)
2370  grdlat2 = cell_corner_lat(next_n)
2371  grdlon2 = cell_corner_lon(next_n)
2372 
2373  lensqr2 = (grdlat1-grdlat2)*(grdlat1-grdlat2) +
2374  & (grdlon1-grdlon2)*(grdlon1-grdlon2)
2375 
2376  if (lensqr2 .le. tiny*tiny) cycle ! degenerate edge
2377 
2378  ! No edge can span more than pi radians
2379 
2380  if (grdlon2-grdlon1 > pi) then
2381  grdlon2 = grdlon2 - pi2
2382  else if (grdlon2-grdlon1 < -pi) then
2383  grdlon1 = grdlon1 - pi2
2384  endif
2385 
2386  ! Also the two intersecting segments together
2387  ! cannot span more than 2*pi radians
2388 
2389  minlon = min(lon1,lon2)
2390  maxlon = max(grdlon1,grdlon2)
2391  if (maxlon-minlon > pi2) then
2392  grdlon1 = grdlon1 - pi2
2393  grdlon2 = grdlon2 - pi2
2394  else
2395  minlon = min(grdlon1,grdlon2)
2396  maxlon = max(lon1,lon2)
2397  if (maxlon-minlon > pi2) then
2398  grdlon1 = grdlon1 + pi2
2399  grdlon2 = grdlon2 + pi2
2400  endif
2401  endif
2402 
2403 
2404  !***
2405  !*** set up linear system to solve for intersection
2406  !***
2407 
2408  mat1 = lat2 - lat1
2409  mat2 = grdlat1 - grdlat2
2410  mat3 = lon2 - lon1
2411  mat4 = grdlon1 - grdlon2
2412  rhs1 = grdlat1 - lat1
2413  rhs2 = grdlon1 - lon1
2414 
2415  determ = mat1*mat4 - mat2*mat3
2416 
2417  !***
2418  !*** if the determinant is zero, the segments are either
2419  !*** parallel or coincident. coincidences were detected
2420  !*** above so do nothing.
2421 
2422  if (abs(determ) > tiny*tiny) then
2423 
2424  !*** if the determinant is non-zero, solve for the linear
2425  !*** parameters s for the intersection point on each line
2426  !*** segment.
2427  !*** if 0<s1,s2<1 then the segment intersects with this side.
2428  !*** return the point of intersection (adding a small
2429  !*** number so the intersection is off the grid line).
2430  !***
2431 
2432  s1 = (rhs1*mat4 - mat2*rhs2)/determ
2433  s2 = (mat1*rhs2 - rhs1*mat3)/determ
2434 
2435  if (s2 >= zero .and. s2 <= one .and.
2436  & s1 > zero .and. s1 <= one) then
2437 
2438  !***
2439  !*** recompute intersection based on full segment
2440  !*** so intersections are consistent for both sweeps
2441  !***
2442 
2443  if (lon2-begsegloc(2) > pi) then
2444  lon2 = lon2 - pi2
2445  else if (lon2-begsegloc(2) < -pi) then
2446  begsegloc(2) = begsegloc(2) - pi2
2447  endif
2448 
2449 
2450  ! Also the two intersecting segments together
2451  ! cannot span more than 2*pi radians
2452 
2453  minlon = min(begsegloc(2),lon2)
2454  maxlon = max(grdlon1,grdlon2)
2455  if (maxlon-minlon > pi2) then
2456  grdlon1 = grdlon1 - pi2
2457  grdlon2 = grdlon2 - pi2
2458  else
2459  minlon = min(grdlon1,grdlon2)
2460  maxlon = max(begsegloc(2),lon2)
2461  if (maxlon-minlon > pi2) then
2462  grdlon1 = grdlon1 + pi2
2463  grdlon2 = grdlon2 + pi2
2464  endif
2465  endif
2466 
2467 
2468  mat1 = lat2 - begsegloc(1)
2469  mat3 = lon2 - begsegloc(2)
2470  rhs1 = grdlat1 - begsegloc(1)
2471  rhs2 = grdlon1 - begsegloc(2)
2472 
2473  determ = mat1*mat4 - mat2*mat3
2474 
2475  !***
2476  !*** sometimes due to roundoff, the previous
2477  !*** determinant is non-zero, but the lines
2478  !*** are actually coincident. if this is the
2479  !*** case, skip the rest.
2480  !***
2481 
2482  if (determ /= zero) then
2483  s1 = (rhs1*mat4 - mat2*rhs2)/determ
2484  s2 = (mat1*rhs2 - rhs1*mat3)/determ
2485 
2486  intrsct_lat = begsegloc(1) + mat1*s1
2487  intrsct_lon = begsegloc(2) + mat3*s1
2488 
2489  if (intrsct_lon < 0.0) then
2490  intrsct_lon = intrsct_lon + pi2
2491  else if (intrsct_lon > pi2) then
2492  intrsct_lon = intrsct_lon - pi2
2493  endif
2494 
2495  !***
2496  !*** Make sure the intersection point is not within
2497  !*** tolerance of the starting point
2498  !***
2499 
2500  if (first) then
2501  max_intrsct_lat = intrsct_lat
2502  max_intrsct_lon = intrsct_lon
2503 
2504  vec1_lat = intrsct_lat-beglat
2505  vec1_lon = intrsct_lon-beglon
2506  if (vec1_lon > pi) then
2507  vec1_lon = vec1_lon - pi2
2508  else if (vec1_lon < -pi) then
2509  vec1_lon = vec1_lon + pi2
2510  endif
2511 
2512  maxdist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
2513  dist2 = maxdist2
2514 
2515  denom = (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
2516  sinang2 = determ*determ/denom
2517  intedge = n
2518  first = .false.
2519  else
2520  vec1_lat = intrsct_lat-beglat
2521  vec1_lon = intrsct_lon-beglon
2522  if (vec1_lon > pi) then
2523  vec1_lon = vec1_lon - pi2
2524  else if (vec1_lon < -pi) then
2525  vec1_lon = vec1_lon + pi2
2526  endif
2527 
2528  dist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
2529 
2530  if (dist2 > maxdist2) then
2531  if (begedge .eq. 0 .or. begedge .ne. n) then
2532  max_intrsct_lat = intrsct_lat
2533  max_intrsct_lon = intrsct_lon
2534  maxdist2 = dist2
2535 
2536  denom =
2537  & (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
2538  sinang2 = determ*determ/denom
2539  intedge = n
2540  endif
2541  endif
2542  endif
2543 
2544  else
2545  print *, 'DEBUG: zero determ'
2546  stop
2547  endif
2548 
2549  endif
2550 
2551  else
2552 
2553  !***
2554  !*** Coincident lines or parallel lines
2555  !***
2556 
2557  cross_product = mat2*rhs2 - mat4*rhs1
2558 
2559  !***
2560  !*** If area of triangle formed by endlat,endlon and
2561  !*** the gridline is negligible then the lines are coincident
2562  !***
2563 
2564 
2565  if (abs(cross_product) < tiny) then
2566 
2567  dot_product = mat1*(-mat2) + mat3*(-mat4)
2568 
2569  lensqr1 = mat1*mat1 + mat3*mat3 ! length sqrd of input
2570  ! segment
2571 
2572  if (dot_product < zero) then
2573 
2574  !***
2575  !*** Segments oriented in the same direction
2576  !***
2577 
2578 
2579  tmplat = grdlat2
2580  tmplon = grdlon2
2581  grdlat2 = grdlat1
2582  grdlon2 = grdlon1
2583  grdlat1 = tmplat
2584  grdlon1 = tmplon
2585 
2586  endif
2587 
2588 
2589  vec2_lat = grdlat1 - lat1
2590  vec2_lon = grdlon1 - lon1
2591  if (vec2_lon > pi) vec2_lon = vec2_lon - pi2
2592  if (vec2_lon < -pi) vec2_lon = vec2_lon + pi2
2593 
2594  lensqr2 = vec2_lat*vec2_lat + vec2_lon*vec2_lon
2595 
2596  if (vec2_lat*mat1 + vec2_lon*mat3 < 0) then
2597  lensqr2 = -lensqr2
2598  endif
2599 
2600  vec3_lat = grdlat2 - lat1
2601  vec3_lon = grdlon2 - lon1
2602  if (vec3_lon > pi) vec3_lon = vec3_lon - pi2
2603  if (vec3_lon < -pi) vec3_lon = vec3_lon + pi2
2604 
2605  lensqr3 = (vec3_lat*vec3_lat+vec3_lon*vec3_lon)
2606 
2607  if (vec3_lat*mat1 + vec3_lon*mat3 < 0) then
2608  lensqr3 = -lensqr3
2609  endif
2610 
2611  found = .false.
2612 
2613  if (lensqr2 > 0) then
2614  if (lensqr2 <= lensqr1) then
2615  intrsct_lat = grdlat1
2616  intrsct_lon = grdlon1
2617  found = .true.
2618  endif
2619  else
2620  if (lensqr3 > 0) then
2621  if (lensqr3 > lensqr1) then
2622  intrsct_lat = lat2
2623  intrsct_lon = lon2
2624  found = .true.
2625  else
2626  intrsct_lat = grdlat2
2627  intrsct_lon = grdlon2
2628  found = .true.
2629  endif
2630  endif
2631  endif
2632 
2633  if (found) then
2634 
2635  dist2 = (intrsct_lat-beglat)*(intrsct_lat-beglat)+
2636  & (intrsct_lon-beglon)*(intrsct_lon-beglon)
2637 
2638  !*** Coincidence intersection always wins
2639 
2640  max_intrsct_lat = intrsct_lat
2641  max_intrsct_lon = intrsct_lon
2642  maxdist2 = dist2
2643  sinang2 = 0
2644  intedge = n
2645  lcoinc = .true.
2646 
2647  exit intrsct_loop
2648  endif
2649 
2650  endif
2651 
2652  endif
2653 
2654  !*** restore lon1 and lon2 in case it got modified
2655 
2656  lon1 = beglon
2657  lon2 = endlon
2658  begsegloc(2) = begseg(2)
2659  if ((lon2-lon1) > pi) then
2660  lon2 = lon2 - pi2
2661  else if ((lon2-lon1) < -pi) then
2662  lon1 = lon1 - pi2
2663  endif
2664 
2665  end do intrsct_loop
2666 
2667  if (intedge .eq. 0) then
2668  return
2669  else
2670  if (maxdist2 < 1e6*tiny*tiny) then
2671  intedge = 0
2672  return
2673  else
2674  intrsct_lat = max_intrsct_lat
2675  intrsct_lon = max_intrsct_lon
2676  endif
2677  endif
2678 
2679 !-----------------------------------------------------------------------
2680 !
2681 ! if the segment crosses a pole threshold, reset the intersection
2682 ! to be the threshold latitude. only check if this was not a
2683 ! threshold segment since sometimes coordinate transform can end
2684 ! up on other side of threshold again.
2685 !
2686 !-----------------------------------------------------------------------
2687 
2688  if (lthresh) then
2689  if (intrsct_lat < north_thresh .or. intrsct_lat > south_thresh)
2690  & lthresh = .false.
2691  else if (lat1 > zero .and. intrsct_lat > north_thresh) then
2692 ! intrsct_lat = north_thresh + tiny
2693  intrsct_lat = north_thresh
2694  mat1 = lat2 - begsegloc(1)
2695  mat3 = lon2 - begsegloc(2)
2696  s1 = (intrsct_lat - begsegloc(1))/mat1
2697  intrsct_lon = begsegloc(2) + s1*mat3
2698  lthresh = .true.
2699  else if (lat1 < zero .and. intrsct_lat < south_thresh) then
2700 ! intrsct_lat = south_thresh - tiny
2701  intrsct_lat = south_thresh
2702  mat1 = lat2 - begsegloc(1)
2703  mat3 = lon2 - begsegloc(2)
2704  s1 = (intrsct_lat - begsegloc(1))/mat1
2705  intrsct_lon = begsegloc(2) + s1*mat3
2706  lthresh = .true.
2707  endif
2708 
2709  if (intrsct_lon < 0.0) then
2710  intrsct_lon = intrsct_lon + pi2
2711  else if (intrsct_lon > pi2) then
2712  intrsct_lon = intrsct_lon - pi2
2713  endif
2714 
2715 
2716 !-----------------------------------------------------------------------
2717 

References scrip_grids::north_thresh, scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, pole_intersection(), scrip_grids::south_thresh, scrip_constants::tiny, and scrip_constants::zero.

Referenced by cell_integrate().

◆ line_integral()

subroutine scrip_remap_conservative::line_integral ( integer (scrip_i4), intent(in)  phi_or_theta,
real (scrip_r8), dimension(2*num_wts), intent(out)  weights,
integer (scrip_i4), intent(in)  num_wts,
real (scrip_r8), intent(in)  in_phi1,
real (scrip_r8), intent(in)  in_phi2,
real (scrip_r8), intent(in)  theta1,
real (scrip_r8), intent(in)  theta2,
real (scrip_r8), intent(in)  grid1_lat,
real (scrip_r8), intent(in)  grid1_lon,
real (scrip_r8), intent(in)  grid2_lat,
real (scrip_r8), intent(in)  grid2_lon 
)

Definition at line 3460 of file scrip_remap_conservative.f.

3460 
3461 !-----------------------------------------------------------------------
3462 !
3463 ! this routine computes the line integral of the flux function
3464 ! that results in the interpolation weights. the line is defined
3465 ! by the input lat/lon of the endpoints.
3466 !
3467 !-----------------------------------------------------------------------
3468 
3469 !-----------------------------------------------------------------------
3470 !
3471 ! intent(in):
3472 !
3473 !-----------------------------------------------------------------------
3474 
3475  integer (SCRIP_i4), intent(in) ::
3476  & phi_or_theta ! Integration variable (lat or lon)
3477 
3478  integer (SCRIP_i4), intent(in) ::
3479  & num_wts ! number of weights to compute
3480 
3481  real (SCRIP_r8), intent(in) ::
3482  & in_phi1, in_phi2, ! longitude endpoints for the segment
3483  & theta1, theta2, ! latitude endpoints for the segment
3484  & grid1_lat, grid1_lon, ! reference coordinates for each
3485  & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
3486 
3487 !-----------------------------------------------------------------------
3488 !
3489 ! intent(out):
3490 !
3491 !-----------------------------------------------------------------------
3492 
3493  real (SCRIP_r8), dimension(2*num_wts), intent(out) ::
3494  & weights ! line integral contribution to weights
3495 
3496 
3497 ! write(*,*)'subroutine line_integral'
3498  if (phi_or_theta .eq. 1) then
3499  call line_integral_phi(weights, num_wts, in_phi1, in_phi2,
3500  & theta1, theta2, grid1_lat, grid1_lon,
3501  & grid2_lat, grid2_lon)
3502  else
3503  call line_integral_theta(weights, num_wts,in_phi1,in_phi2,
3504  & theta1, theta2, grid1_lat, grid1_lon,
3505  & grid2_lat, grid2_lon)
3506  endif
3507 
3508 
3509  return
3510 
3511 !-----------------------------------------------------------------------
3512 

Referenced by cell_integrate(), and remap_conserv().

◆ line_integral_phi()

subroutine scrip_remap_conservative::line_integral_phi ( real (scrip_r8), dimension(2*num_wts), intent(out)  weights,
integer (scrip_i4), intent(in)  num_wts,
real (scrip_r8), intent(in)  in_phi1,
real (scrip_r8), intent(in)  in_phi2,
real (scrip_r8), intent(in)  theta1,
real (scrip_r8), intent(in)  theta2,
real (scrip_r8), intent(in)  grid1_lat,
real (scrip_r8), intent(in)  grid1_lon,
real (scrip_r8), intent(in)  grid2_lat,
real (scrip_r8), intent(in)  grid2_lon 
)

Definition at line 3522 of file scrip_remap_conservative.f.

3522 
3523 !-----------------------------------------------------------------------
3524 !
3525 ! this routine computes the line integral of the flux function
3526 ! that results in the interpolation weights. the line is defined
3527 ! by the input lat/lon of the endpoints. Integration is w.r.t. lon
3528 !
3529 !-----------------------------------------------------------------------
3530 
3531 !-----------------------------------------------------------------------
3532 !
3533 ! intent(in):
3534 !
3535 !-----------------------------------------------------------------------
3536 
3537  integer (SCRIP_i4), intent(in) ::
3538  & num_wts ! number of weights to compute
3539 
3540  real (SCRIP_r8), intent(in) ::
3541  & in_phi1, in_phi2, ! longitude endpoints for the segment
3542  & theta1, theta2, ! latitude endpoints for the segment
3543  & grid1_lat, grid1_lon, ! reference coordinates for each
3544  & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
3545 
3546 !-----------------------------------------------------------------------
3547 !
3548 ! intent(out):
3549 !
3550 !-----------------------------------------------------------------------
3551 
3552  real (SCRIP_r8), dimension(2*num_wts), intent(out) ::
3553  & weights ! line integral contribution to weights
3554 
3555 !-----------------------------------------------------------------------
3556 !
3557 ! local variables
3558 !
3559 !-----------------------------------------------------------------------
3560 
3561  real (SCRIP_r8) :: dphi, sinth1, sinth2, costh1, costh2, fac,
3562  & phi1, phi2
3563  real (SCRIP_r8) :: f1, f2, fint
3564 
3565 !-----------------------------------------------------------------------
3566 !
3567 ! weights for the general case based on a trapezoidal approx to
3568 ! the integrals.
3569 !
3570 !-----------------------------------------------------------------------
3571 
3572 
3573 ! write(*,*)'subroutine line_integral_phi'
3574 
3575  sinth1 = sin(theta1)
3576  sinth2 = sin(theta2)
3577  costh1 = cos(theta1)
3578  costh2 = cos(theta2)
3579 
3580  dphi = in_phi1 - in_phi2
3581  if (dphi > pi) then
3582  dphi = dphi - pi2
3583  else if (dphi < -pi) then
3584  dphi = dphi + pi2
3585  endif
3586  dphi = half*dphi
3587 
3588 !-----------------------------------------------------------------------
3589 !
3590 ! the first weight is the area overlap integral. the second and
3591 ! fourth are second-order latitude gradient weights.
3592 !
3593 !-----------------------------------------------------------------------
3594 
3595  weights( 1) = dphi*(sinth1 + sinth2)
3596  write(401,*)weights(1),' % A'
3597  weights(num_wts+1) = dphi*(sinth1 + sinth2)
3598  weights( 2) = dphi*(costh1 + costh2 + (theta1*sinth1 +
3599  & theta2*sinth2))
3600  weights(num_wts+2) = dphi*(costh1 + costh2 + (theta1*sinth1 +
3601  & theta2*sinth2))
3602 
3603 !-----------------------------------------------------------------------
3604 !
3605 ! the third and fifth weights are for the second-order phi gradient
3606 ! component. must be careful of longitude range.
3607 !
3608 !-----------------------------------------------------------------------
3609 
3610  f1 = half*(costh1*sinth1 + theta1)
3611  f2 = half*(costh2*sinth2 + theta2)
3612 
3613  phi1 = in_phi1 - grid1_lon
3614  if (phi1 > pi) then
3615  phi1 = phi1 - pi2
3616  else if (phi1 < -pi) then
3617  phi1 = phi1 + pi2
3618  endif
3619 
3620  phi2 = in_phi2 - grid1_lon
3621  if (phi2 > pi) then
3622  phi2 = phi2 - pi2
3623  else if (phi2 < -pi) then
3624  phi2 = phi2 + pi2
3625  endif
3626 
3627  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3628  weights(3) = dphi*(phi1*f1 + phi2*f2)
3629  else
3630  if (phi1 > zero) then
3631  fac = pi
3632  else
3633  fac = -pi
3634  endif
3635  fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
3636  weights(3) = half*phi1*(phi1-fac)*f1 -
3637  & half*phi2*(phi2+fac)*f2 +
3638  & half*fac*(phi1+phi2)*fint
3639  endif
3640 
3641  phi1 = in_phi1 - grid2_lon
3642  if (phi1 > pi) then
3643  phi1 = phi1 - pi2
3644  else if (phi1 < -pi) then
3645  phi1 = phi1 + pi2
3646  endif
3647 
3648  phi2 = in_phi2 - grid2_lon
3649  if (phi2 > pi) then
3650  phi2 = phi2 - pi2
3651  else if (phi2 < -pi) then
3652  phi2 = phi2 + pi2
3653  endif
3654 
3655  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3656  weights(num_wts+3) = dphi*(phi1*f1 + phi2*f2)
3657  else
3658  if (phi1 > zero) then
3659  fac = pi
3660  else
3661  fac = -pi
3662  endif
3663  fint = f1 + (f2-f1)*(fac-phi1)/abs(dphi)
3664  weights(num_wts+3) = half*phi1*(phi1-fac)*f1 -
3665  & half*phi2*(phi2+fac)*f2 +
3666  & half*fac*(phi1+phi2)*fint
3667  endif
3668 
3669 !-----------------------------------------------------------------------
3670 

References scrip_constants::half, scrip_constants::pi, scrip_constants::pi2, and scrip_constants::zero.

◆ line_integral_theta()

subroutine scrip_remap_conservative::line_integral_theta ( real (scrip_r8), dimension(2*num_wts), intent(out)  weights,
integer (scrip_i4), intent(in)  num_wts,
real (scrip_r8), intent(in)  in_phi1,
real (scrip_r8), intent(in)  in_phi2,
real (scrip_r8), intent(in)  theta1,
real (scrip_r8), intent(in)  theta2,
real (scrip_r8), intent(in)  grid1_lat,
real (scrip_r8), intent(in)  grid1_lon,
real (scrip_r8), intent(in)  grid2_lat,
real (scrip_r8), intent(in)  grid2_lon 
)

Definition at line 3682 of file scrip_remap_conservative.f.

3682 
3683 !-----------------------------------------------------------------------
3684 !
3685 ! this routine computes the line integral of the flux function
3686 ! that results in the interpolation weights. the line is defined
3687 ! by the input lat/lon of the endpoints. Integration is w.r.t. lat
3688 !
3689 ! Needed to use Simpson rule for this integration to get lower errors
3690 !-----------------------------------------------------------------------
3691 
3692 !-----------------------------------------------------------------------
3693 !
3694 ! intent(in):
3695 !
3696 !-----------------------------------------------------------------------
3697 
3698  integer (SCRIP_i4), intent(in) ::
3699  & num_wts ! number of weights to compute
3700 
3701  real (SCRIP_r8), intent(in) ::
3702  & in_phi1, in_phi2, ! longitude endpoints for the segment
3703  & theta1, theta2, ! latitude endpoints for the segment
3704  & grid1_lat, grid1_lon, ! reference coordinates for each
3705  & grid2_lat, grid2_lon ! grid (to ensure correct 0,2pi interv.
3706 
3707 !-----------------------------------------------------------------------
3708 !
3709 ! intent(out):
3710 !
3711 !-----------------------------------------------------------------------
3712 
3713  real (SCRIP_r8), dimension(2*num_wts), intent(out) ::
3714  & weights ! line integral contribution to weights
3715 
3716 !-----------------------------------------------------------------------
3717 !
3718 ! local variables
3719 !
3720 !-----------------------------------------------------------------------
3721 
3722  real (SCRIP_r8) :: dtheta, dtheta2, costh1, costh2, costhpi,
3723  & phi1, phi2, theta_pi, f1, f2, fpi,
3724  & fm, costhm, part1, part2
3725 
3726 !-----------------------------------------------------------------------
3727 !
3728 ! weights for the general case based on a trapezoidal approx to
3729 ! the integrals.
3730 !
3731 !-----------------------------------------------------------------------
3732 
3733  costh1 = cos(theta1)
3734  costh2 = cos(theta2)
3735  costhm = cos(half*(theta1+theta2))
3736 
3737  dtheta = theta2 - theta1
3738  dtheta2 = half*dtheta
3739 
3740 ! write(*,*)' subroutine line_integral_theta'
3741 
3742 
3743 !-----------------------------------------------------------------------
3744 !
3745 ! Need to account for double value of longitude in calculations of
3746 ! all the weights. First we transform all the phis to be relative
3747 ! to the grid center This takes care of a good number of cases where
3748 ! the the phis span the periodic boundary in the longitudinal
3749 ! direction. If we still have a line that spans the periodic
3750 ! boundary then we have to integrate along the line in two parts -
3751 ! from point 1 to the periodic boundary and from the periodic
3752 ! boundary to the second point
3753 !
3754 ! Example: Consider a line which has points at phi1 = -100 and phi2
3755 ! = 100 degrees and say the grid center is at phi_c = 0
3756 ! degrees. Then phi1-phi_c > -180 and phi2-phi_c < 180. But
3757 ! phi2-phi1 > 180.
3758 !
3759 ! *********************************************!!!!!!!!!!!
3760 ! If we are doing the second step anyway, why are we normalizing the
3761 ! coordinates with respect to the grid centers?
3762 !
3763 ! We need it particularly in this integration because phi figures
3764 ! explicitly in the expressions - so if a cell straddles the 0,2pi
3765 ! boundary, we integrate some edges with phi values close to zero
3766 ! and others with phi values close to 2pi leading to errors
3767 ! *********************************************!!!!!!!!!!!
3768 !
3769 !-----------------------------------------------------------------------
3770 
3771  phi1 = in_phi1 - grid1_lon
3772  if (phi1 > pi) then
3773  phi1 = phi1 - pi2
3774  else if (phi1 < -pi) then
3775  phi1 = phi1 + pi2
3776  endif
3777 
3778  phi2 = in_phi2 - grid1_lon
3779  if (phi2 > pi) then
3780  phi2 = phi2 - pi2
3781  else if (phi2 < -pi) then
3782  phi2 = phi2 + pi2
3783  endif
3784 
3785  f1 = phi1*costh1
3786  f2 = phi2*costh2
3787 
3788  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3789 
3790  fm = half*(phi1+phi2)*costhm
3791 
3792  weights(1) = dtheta*(f1 + 4*fm + f2)/6.0
3793 ! write(401,*)weights(1),' % A'
3794 
3795  weights(2) = dtheta2*(theta1*f1 + theta2*f2)
3796 
3797  weights(3) = half*dtheta2*(f1*f1 + f2*f2)
3798 
3799  else
3800  if (phi1 > zero) then ! Means phi2-phi1 < -pi
3801 
3802 ! theta at phi = pi
3803  theta_pi = theta1 + (pi - phi1)*dtheta/(phi2 + pi2 - phi1)
3804 ! print *, ''
3805 ! print *, 'phi1',phi1,' phi2',phi2
3806 ! print *, 'theta1',theta1,' theta2',theta2
3807 ! print *, 'theta_pi',theta_pi
3808 
3809  costhpi = cos(theta_pi)
3810  fpi = pi*costhpi
3811 
3812  fm = half*(phi1+pi)*cos(half*(theta1+theta_pi))
3813  part1 = (theta_pi - theta1)*(f1 + 4*fm + fpi)/6.0
3814 
3815  fm = half*(phi2-pi)*cos(half*(theta1+theta_pi))
3816  part2 = 0.5*(theta2 - theta_pi)*(-fpi + 4*fm + f2)/6.0
3817 
3818  weights(1) = part1 + part2
3819 ! write(401,*)weights(1),' % B'
3820 
3821  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 + theta_pi*fpi)
3822  part2 = 0.5*(theta2 - theta_pi)*(-theta_pi*fpi + theta2*f2)
3823  weights(2) = part1 + part2
3824 
3825 
3826  else ! Means phi2-phi1 > pi
3827 
3828 ! theta at phi = -pi
3829  theta_pi = theta1 + (-pi - phi1)*dtheta/(phi2 - pi2 - phi1)
3830 ! print *, ''
3831 ! print *, 'phi1',phi1,' phi2',phi2
3832 ! print *, 'theta1',theta1,' theta2',theta2
3833 ! print *, 'theta_pi',theta_pi
3834 
3835  costhpi = cos(theta_pi)
3836  fpi = pi*costhpi
3837 
3838  fm = half*(phi1-pi)*cos(half*(theta1+theta_pi))
3839  part1 = 0.5*(theta_pi - theta1)*(f1 + 4*fm - fpi)/6.0
3840 
3841  fm = half*(pi+phi2)*cos(half*(theta2+theta_pi))
3842  part2 = 0.5*(theta2 - theta_pi)*(fpi + 4*fm + f2)/6.0
3843  weights(1) = part1 + part2
3844 ! write(401,*)weights(1),' % C'
3845 
3846  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 - theta_pi*fpi)
3847  part2 = 0.5*(theta2 - theta_pi)*(theta_pi*fpi + theta2*f2)
3848  weights(2) = part1 + part2
3849 
3850 
3851  endif
3852 
3853  part1 = 0.25*(theta_pi - theta1)*(f1*f1 + fpi*fpi)
3854  part2 = 0.25*(theta2 - theta_pi)*(fpi*fpi + f2*f2)
3855  weights(3) = part1 + part2
3856 
3857  endif
3858 
3859 
3860  phi1 = in_phi1 - grid2_lon
3861  if (phi1 > pi) then
3862  phi1 = phi1 - pi2
3863  else if (phi1 < -pi) then
3864  phi1 = phi1 + pi2
3865  endif
3866 
3867  phi2 = in_phi2 - grid2_lon
3868  if (phi2 > pi) then
3869  phi2 = phi2 - pi2
3870  else if (phi2 < -pi) then
3871  phi2 = phi2 + pi2
3872  endif
3873 
3874 
3875  f1 = phi1*costh1
3876  f2 = phi2*costh2
3877 
3878  if ((phi2-phi1) < pi .and. (phi2-phi1) > -pi) then
3879 
3880  fm = half*(phi1+phi2)*costhm
3881 
3882  weights(num_wts+1) = dtheta2*(f1 + f2)
3883 
3884  weights(num_wts+2) = dtheta2*(theta1*f1 + theta2*f2)
3885 
3886  weights(num_wts+3) = half*dtheta2*(f1*f1 + f2*f2)
3887 
3888  else
3889  if (phi1 > zero) then
3890 
3891  theta_pi = theta1 + (pi - phi1)*dtheta/(phi2 + pi2 - phi1)
3892 ! print *, ''
3893 ! print *, 'phi1',phi1,' phi2',phi2
3894 ! print *, 'theta1',theta1,' theta2',theta2
3895 ! print *, 'theta_pi',theta_pi
3896 
3897  costhpi = cos(theta_pi)
3898  fpi = pi*costhpi
3899 
3900  fm = half*(phi1+pi)*cos(half*(theta1+theta_pi))
3901  part1 = (theta_pi - theta1)*(f1 + 4*fm + fpi)/6.0
3902 
3903  fm = half*(-pi+phi2)*cos(half*(theta2+theta_pi))
3904  part2 = (theta2 - theta_pi)*(-fpi + 4*fm + f2)/6.0
3905  weights(num_wts+1) = part1 + part2
3906 
3907  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 + theta_pi*fpi)
3908  part2 = 0.5*(theta2 - theta_pi)*(-theta_pi*fpi + theta2*f2)
3909  weights(num_wts+2) = part1 + part2
3910 
3911 
3912  else
3913 
3914  theta_pi = theta1 + (-pi - phi1)*dtheta/(phi2 - pi2 - phi1)
3915 ! print *, ''
3916 ! print *, 'phi1',phi1,' phi2',phi2
3917 ! print *, 'theta1',theta1,' theta2',theta2
3918 ! print *, 'theta_pi',theta_pi
3919 
3920  costhpi = cos(theta_pi)
3921  fpi = pi*costhpi
3922 
3923  fm = half*(phi1-pi)*cos(half*(theta1+theta_pi))
3924  part1 = (theta_pi - theta1)*(f1 +4*fm - fpi)/6.0
3925 
3926  fm = half*(phi2+pi)*cos(half*(theta2+theta_pi))
3927  part2 = 0.5*(theta2 - theta_pi)*(fpi + 4*fm + f2)/6.0
3928  weights(num_wts+1) = part1 + part2
3929 
3930  part1 = 0.5*(theta_pi - theta1)*(theta1*f1 - theta_pi*fpi)
3931  part2 = 0.5*(theta2 - theta_pi)*(theta_pi*fpi + theta2*f2)
3932  weights(num_wts+2) = part1 + part2
3933 
3934  endif
3935 
3936  part1 = 0.25*(theta_pi - theta1)*(f1*f1 + fpi*fpi)
3937  part2 = 0.25*(theta2 - theta_pi)*(fpi*fpi + f2*f2)
3938  weights(num_wts+3) = part1 + part2
3939 
3940  endif
3941 
3942 !-----------------------------------------------------------------------
3943 

References scrip_constants::half, scrip_constants::pi, scrip_constants::pi2, and scrip_constants::zero.

◆ locate_point()

subroutine scrip_remap_conservative::locate_point ( real (scrip_r8), intent(in)  ptlat,
real (scrip_r8), intent(in)  ptlon,
integer (scrip_i4), intent(in)  cell,
integer (scrip_i4), intent(in)  cell_grid_num,
integer (scrip_i4), intent(in)  srch_grid_num,
integer (scrip_i4), intent(out)  cont_cell,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 4423 of file scrip_remap_conservative.f.

4423 
4424 !-----------------------------------------------------------------------
4425 !
4426 ! Find the cell containing the given point
4427 !
4428 !-----------------------------------------------------------------------
4429 
4430  implicit none
4431 
4432 !-----------------------------------------------------------------------
4433 !
4434 ! intent(in):
4435 !
4436 !-----------------------------------------------------------------------
4437 
4438  real (SCRIP_r8), intent(in) ::
4439  & ptlat, ptlon ! Point to locate
4440 
4441  integer (SCRIP_i4), intent(in) ::
4442  & cell, ! Cell from which point originates
4443  ! Point will be on boundary of orig_cell
4444  & cell_grid_num ! Index of grid to which cell belongs
4445 
4446  integer (SCRIP_i4), intent(in) ::
4447  & srch_grid_num ! num indicating if we are locating a
4448  ! grid1 point in a cell of grid2 (num=2)
4449  ! or a grid2 point in a cell of grid1
4450  ! (num=1)
4451 
4452 !-----------------------------------------------------------------------
4453 !
4454 ! intent(out):
4455 !
4456 !-----------------------------------------------------------------------
4457 
4458  integer (SCRIP_i4), intent(out) ::
4459  & cont_cell ! grid cell containing this point
4460 
4461  logical (SCRIP_logical), intent(out) ::
4462  & lboundary ! flag points that lie on the boundary
4463  ! of the cell
4464 
4465  integer (SCRIP_i4), intent(out) ::
4466  & edgeid ! if point is on boundary, which local
4467  ! edge is it on? (0 otherwise)
4468 
4469 !-----------------------------------------------------------------------
4470 !
4471 ! local variables
4472 !
4473 !-----------------------------------------------------------------------
4474 
4475  integer (SCRIP_i4) :: i, j, n, ic
4476  integer (SCRIP_i4) :: whichpole, srch_cell_add,
4477  & grid1_add, grid2_add, min_add, max_add,
4478  & previdx, nextidx, pcorner, corner,
4479  & ncorners, nalloc
4480 
4481  real (SCRIP_r8), dimension(:), allocatable ::
4482  & cell_corner_lat,
4483  & cell_corner_lon
4484 
4485  real (SCRIP_r8) ::
4486  & prevlon,
4487  & nextlon,
4488  & polelat,
4489  & cell_center_lat,
4490  & cell_center_lon
4491 
4492 
4493  logical (SCRIP_logical) :: inpoly, latlon
4494  logical (SCRIP_logical) :: test
4495 
4496 !-----------------------------------------------------------------------
4497 !
4498 ! initialize defaults, flags, etc.
4499 !
4500 !-----------------------------------------------------------------------
4501 
4502  lboundary = .false.
4503  edgeid = 0
4504  cont_cell = 0
4505 
4506  if (cell /= last_cell_locate_point .or. cell_grid_num /=
4507  & last_cell_grid_num_locate_point
4508  & .or. srch_grid_num /= last_srch_grid_num_locate_point) then
4509 
4510  last_cell_locate_point = cell
4511  last_cell_grid_num_locate_point = cell_grid_num
4512  last_srch_grid_num_locate_point = srch_grid_num
4513 
4514  if (first_call_locate_point) then
4515  first_call_locate_point = .false.
4516  last_cell_locate_point = 0
4517  last_cell_grid_num_locate_point = 0
4518  last_srch_grid_num_locate_point = 0
4519  num_srch_cell_locate_points = 0
4520  else
4521  if (num_srch_cell_locate_points .gt. 0) then
4522  deallocate(srch_add_locate_point,
4523  & srch_corner_lat_locate_point,
4524  & srch_corner_lon_locate_point)
4525  endif
4526  endif
4527 
4528  call get_srch_cells(cell, cell_grid_num, srch_grid_num,
4529  & num_srch_cell_locate_points, srch_add_locate_point,
4530  & srch_corners_locate_point,
4531  & srch_corner_lat_locate_point,srch_corner_lon_locate_point,
4532  & srch_center_lat_locate_point,srch_center_lon_locate_point)
4533 
4534  endif
4535 
4536  if (num_srch_cell_locate_points == 0) return
4537 
4538 
4539  ncorners = srch_corners_locate_point
4540  nalloc = ncorners+2
4541  allocate(cell_corner_lat(nalloc),
4542  & cell_corner_lon(nalloc))
4543 
4544 
4545  do ic=1,num_srch_cell_locate_points
4546 
4547  srch_cell_add = srch_add_locate_point(ic)
4548 
4549  do i = 1, ncorners
4550  cell_corner_lat(i) = srch_corner_lat_locate_point(i,ic)
4551  cell_corner_lon(i) = srch_corner_lon_locate_point(i,ic)
4552  enddo
4553 
4554  cell_center_lat = srch_center_lat_locate_point(ic)
4555  cell_center_lon = srch_center_lon_locate_point(ic)
4556 
4557 ! if ((srch_grid_num .eq. 1 .and.
4558 ! & (special_polar_cell1(srch_cell_add))) .or.
4559 ! & (srch_grid_num .eq. 2 .and.
4560 ! & (special_polar_cell2(srch_cell_add)))) then
4561 !
4562 ! Modified by MD
4563  test=.false.
4564  if (srch_grid_num .eq. 1) then
4565  if (special_polar_cell1(srch_cell_add)) then
4566  test=.true.
4567  endif
4568  else
4569  if (special_polar_cell2(srch_cell_add)) then
4570  test=.true.
4571  endif
4572  endif
4573  if (test) then
4574  call modify_polar_cell(ncorners, nalloc, cell_corner_lat,
4575  & cell_corner_lon)
4576 
4577  endif
4578 
4579  call ptincell(ptlat, ptlon, srch_cell_add, ncorners,
4580  & cell_corner_lat, cell_corner_lon,
4581  & cell_center_lat, cell_center_lon,
4582  & srch_grid_num, inpoly, lboundary, edgeid)
4583 
4584 
4585  if (inpoly) then
4586  cont_cell = srch_cell_add
4587  exit
4588  endif
4589 
4590  ncorners = srch_corners_locate_point ! reset it for other srch
4591  !cells
4592  end do
4593 
4594 !----------------------------------------------------------------------
4595 

References first_call_locate_point, get_srch_cells(), last_cell_grid_num_locate_point, last_cell_locate_point, last_srch_grid_num_locate_point, modify_polar_cell(), num_srch_cell_locate_points, ptincell(), scrip_grids::special_polar_cell1, scrip_grids::special_polar_cell2, srch_add_locate_point, srch_center_lat_locate_point, srch_center_lon_locate_point, srch_corner_lat_locate_point, srch_corner_lon_locate_point, and srch_corners_locate_point.

Referenced by cell_integrate().

◆ locate_segstart()

subroutine scrip_remap_conservative::locate_segstart ( integer (scrip_i4), intent(in)  cell_grid_num,
integer (scrip_i4), intent(in)  cell,
real (scrip_r8), intent(in)  beglat,
real (scrip_r8), intent(in)  beglon,
real (scrip_r8), intent(in)  endlat,
real (scrip_r8), intent(in)  endlon,
real (scrip_r8), intent(in)  offset,
integer (scrip_i4), intent(in)  srch_grid_num,
integer (scrip_i4), intent(out)  cont_cell,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 4096 of file scrip_remap_conservative.f.

4096 
4097 !-----------------------------------------------------------------------
4098 !
4099 ! Find the cell containing the given point
4100 !
4101 !-----------------------------------------------------------------------
4102 
4103  implicit none
4104 
4105 !-----------------------------------------------------------------------
4106 !
4107 ! intent(in):
4108 !
4109 !-----------------------------------------------------------------------
4110 
4111  real (SCRIP_r8), intent(in) ::
4112  & beglat, beglon, ! beginning and end points of segment
4113  & endlat, endlon ! on which the point to be located lies
4114 
4115  real (SCRIP_r8), intent(in) ::
4116  & offset ! Offset to calculate the search point
4117 
4118  integer (SCRIP_i4), intent(in) ::
4119  & cell, ! Cell from which point originates
4120  ! Point will be on boundary of orig_cell
4121  & cell_grid_num ! Index of grid to which cell belongs
4122 
4123  integer (SCRIP_i4), intent(in) ::
4124  & srch_grid_num ! num indicating if we are locating a
4125  ! grid1 point in a cell of grid2 (num=2)
4126  ! or a grid2 point in a cell of grid1
4127  ! (num=1)
4128 
4129 !-----------------------------------------------------------------------
4130 !
4131 ! intent(out):
4132 !
4133 !-----------------------------------------------------------------------
4134 
4135  integer (SCRIP_i4), intent(out) ::
4136  & cont_cell ! grid cell containing this point
4137 
4138  logical (SCRIP_logical), intent(out) ::
4139  & lboundary ! flag points that lie on the boundary
4140  ! of the cell
4141 
4142  integer (SCRIP_i4), intent(out) ::
4143  & edgeid ! if point is on boundary, which local
4144  ! edge is it on? (0 otherwise)
4145 
4146 !-----------------------------------------------------------------------
4147 !
4148 ! local variables
4149 !
4150 !-----------------------------------------------------------------------
4151 
4152  integer (SCRIP_i4) :: i, j, k, n, ic
4153  integer (SCRIP_i4) :: whichpole, srch_cell_add,
4154  & grid1_add, grid2_add, min_add, max_add
4155 
4156  real (SCRIP_r8), dimension(:), allocatable ::
4157  & cell_corner_x, cell_corner_y
4158 
4159  logical (SCRIP_logical) :: inpoly, latlon
4160 
4161  real (SCRIP_r8) ::
4162  & vec1_x, vec1_y, vec1_lenx, vec1_lat, vec1_lon, vec1_len,
4163  & begx, begy, endx, endy, ptx, pty, rns, pi4, ptlat, ptlon,
4164  & lat, lon, cell_center_x, cell_center_y
4165 
4166 !-----------------------------------------------------------------------
4167 !
4168 ! initialize defaults, flags, etc.
4169 !
4170 !-----------------------------------------------------------------------
4171 
4172  lboundary = .false.
4173  edgeid = 0
4174  cont_cell = 0
4175 
4176 
4177  if (cell /= last_cell_locate_segstart .or.
4178  & cell_grid_num /= last_cell_grid_num_locate_segstart
4179  & .or. srch_grid_num /= last_srch_grid_num_locate_segstart)
4180  & then
4181 
4182  last_cell_locate_segstart = cell
4183  last_cell_grid_num_locate_segstart = cell_grid_num
4184  last_srch_grid_num_locate_segstart = srch_grid_num
4185 
4186  if (first_call_locate_segstart) then
4187  first_call_locate_segstart = .false.
4188  last_cell_locate_segstart = 0
4189  last_cell_grid_num_locate_segstart = 0
4190  last_srch_grid_num_locate_segstart = 0
4191  num_srch_cells_locate_segstart = 0
4192  else
4193  if (num_srch_cells_locate_segstart .gt. 0) then
4194  deallocate(srch_add_locate_segstart,
4195  & srch_corner_lat_locate_segstart,
4196  & srch_corner_lon_locate_segstart,
4197  & srch_center_lat_locate_segstart,
4198  & srch_center_lon_locate_segstart)
4199  endif
4200  endif
4201 
4202  call get_srch_cells(cell, cell_grid_num, srch_grid_num,
4203  & num_srch_cells_locate_segstart, srch_add_locate_segstart,
4204  & srch_corners_locate_segstart,
4205  & srch_corner_lat_locate_segstart,
4206  & srch_corner_lon_locate_segstart,
4207  & srch_center_lat_locate_segstart,
4208  & srch_center_lon_locate_segstart)
4209 
4210  endif
4211 
4212  if (num_srch_cells_locate_segstart == 0) return
4213 
4214 
4215  do ic=1,num_srch_cells_locate_segstart
4216 
4217  srch_cell_add = srch_add_locate_segstart(ic)
4218 
4219 
4220 
4221  !**** CAN WE ACCOMPLISH THE FOLLOWING THROUGH A SUBROUTINE
4222  !**** CALLED SEGSTART_INCELL ??
4223 
4224 
4225  !*** IF POINT IS IN POLAR REGION, CHECK IN A TRANSFORMED SPACE
4226  !*** HOWEVER, POINTS THAT ARE PRACTICALLY AT THE POLE CANNOT
4227  !*** BE CORRECTLY LOCATED THIS WAY BECAUSE THE POLE IS A SINGULARITY
4228  !*** AND CONTAINMENT IN ANY CELL INCIDENT ON THE POLE WILL GIVE US A
4229  !*** POSITIVE ANSWER. FOR THESE POINTS REVERT TO THE LATLON SPACE
4230  !***
4231 
4232 
4233 
4234  vec1_lat = endlat-beglat
4235  vec1_lon = endlon-beglon
4236  if (vec1_lon > pi) then
4237  vec1_lon = vec1_lon - pi2
4238  else if (vec1_lon < -pi) then
4239  vec1_lon = vec1_lon + pi2
4240  endif
4241  vec1_len = sqrt(vec1_lat*vec1_lat+vec1_lon*vec1_lon)
4242  vec1_lat = vec1_lat/vec1_len
4243  vec1_lon = vec1_lon/vec1_len
4244 
4245  ptlat = beglat + offset*vec1_lat
4246  ptlon = beglon + offset*vec1_lon
4247 
4248 
4249  if ((ptlat .gt. north_thresh .and. abs(ptlat-pih) .ge. 0.001)
4250  & .or.
4251  & (ptlat .lt. south_thresh .and. abs(ptlat+pih) .ge. 0.001))
4252  & then
4253 
4254  if (ptlat > zero) then
4255  pi4 = quart*pi
4256  rns = one
4257  else
4258  pi4 = -quart*pi
4259  rns = -one
4260  endif
4261 
4262 
4263  begx = rns*two*sin(pi4 - half*beglat)*cos(beglon)
4264  begy = two*sin(pi4 - half*beglat)*sin(beglon)
4265  endx = rns*two*sin(pi4 - half*endlat)*cos(endlon)
4266  endy = two*sin(pi4 - half*endlat)*sin(endlon)
4267 
4268  vec1_x = endx-begx
4269  vec1_y = endy-begy
4270 
4271  vec1_lenx = sqrt(vec1_x*vec1_x + vec1_y*vec1_y)
4272  vec1_x = vec1_x/vec1_lenx
4273  vec1_y = vec1_y/vec1_lenx
4274 
4275 
4276  !*** Must calculate ptx and pty as an offset on straight
4277  !*** line in polar space rather than calculating it on a
4278  !*** straight line in latlon space an offset point in latlon
4279  !*** space will be off the straight line in polar space
4280 
4281  ptx = begx + offset*vec1_x
4282  pty = begy + offset*vec1_y
4283 
4284  latlon = .false.
4285 
4286  ! Since we want greater fidelity for locating the points
4287  ! we send in the mid-points of the polygon edges too
4288  ! BUT THAT MAKES THE POLYGON NON-CONVEX SOMETIMES AND
4289  ! THE CROSS-PRODUCT CHECK FAILS. SO USE CODE TO CHECK GENERAL
4290  ! POLYGONS
4291 
4292 
4293  allocate(cell_corner_x(npseg*srch_corners_locate_segstart),
4294  & cell_corner_y(npseg*srch_corners_locate_segstart))
4295 
4296 
4297  k = 0
4298  do i = srch_corners_locate_segstart, 1, -1
4299  k = k+1
4300  lat = srch_corner_lat_locate_segstart(i,ic)
4301  lon = srch_corner_lon_locate_segstart(i,ic)
4302  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4303  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4304 
4305  j = i-1
4306  if (j .eq. 0) j = srch_corners_locate_segstart ! how do
4307  ! we do (j-1+n)%n in F90?
4308 
4309  vec1_lat = srch_corner_lat_locate_segstart(j,ic)
4310  & -srch_corner_lat_locate_segstart(i,ic)
4311  vec1_lon = srch_corner_lon_locate_segstart(j,ic)
4312  & -srch_corner_lon_locate_segstart(i,ic)
4313  if (vec1_lon > pi) then
4314  vec1_lon = vec1_lon - pi2
4315  else if (vec1_lon < -pi) then
4316  vec1_lon = vec1_lon + pi2
4317  endif
4318 
4319  do j = 1, npseg-1
4320  k = k+1
4321  lat = srch_corner_lat_locate_segstart(i,ic)
4322  & + j*vec1_lat/npseg
4323  lon = srch_corner_lon_locate_segstart(i,ic)
4324  & + j*vec1_lon/npseg
4325  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4326  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4327  enddo
4328  enddo
4329 
4330 
4331  call ptinpolygen2(ptx, pty, k, cell_corner_x,
4332  & cell_corner_y, latlon, inpoly, lboundary, edgeid)
4333 
4334  if (lboundary) then
4335  edgeid = (edgeid-1)/npseg + 1 ! convert from index in
4336  ! multi-segmented to regular cell
4337  endif
4338 
4339  deallocate(cell_corner_x, cell_corner_y)
4340 
4341  else
4342 
4343  latlon = .true.
4344 
4345  whichpole = 0
4346  if (srch_grid_num .eq. 1 .and.
4347  & srch_cell_add .eq. grid1_spole_cell) then
4348 
4349  whichpole = -1 ! S pole
4350  call ptinpolarpoly(ptlat, ptlon,
4351  & srch_corners_locate_segstart,
4352  & srch_corner_lat_locate_segstart(:,ic),
4353  & srch_corner_lon_locate_segstart(:,ic),
4354  & latlon, whichpole, inpoly, lboundary, edgeid)
4355 
4356  else if (srch_grid_num .eq. 1 .and.
4357  & srch_cell_add .eq. grid1_npole_cell) then
4358 
4359  whichpole = 1 ! N pole
4360  call ptinpolarpoly(ptlat, ptlon,
4361  & srch_corners_locate_segstart,
4362  & srch_corner_lat_locate_segstart(:,ic),
4363  & srch_corner_lon_locate_segstart(:,ic),
4364  & latlon, whichpole, inpoly, lboundary, edgeid)
4365 
4366  else if (srch_grid_num .eq. 2 .and.
4367  & srch_cell_add .eq. grid2_spole_cell) then
4368 
4369  whichpole = -1 ! S pole
4370  call ptinpolarpoly(ptlat, ptlon,
4371  & srch_corners_locate_segstart,
4372  & srch_corner_lat_locate_segstart(:,ic),
4373  & srch_corner_lon_locate_segstart(:,ic),
4374  & latlon, whichpole, inpoly, lboundary, edgeid)
4375 
4376  else if (srch_grid_num .eq. 2 .and.
4377  & srch_cell_add .eq. grid2_npole_cell) then
4378 
4379  whichpole = 1 ! N pole
4380  call ptinpolarpoly(ptlat, ptlon,
4381  & srch_corners_locate_segstart,
4382  & srch_corner_lat_locate_segstart(:,ic),
4383  & srch_corner_lon_locate_segstart(:,ic),
4384  & latlon, whichpole, inpoly, lboundary, edgeid)
4385 
4386  else
4387 
4388  !***
4389  !*** General cell
4390  !***
4391 
4392  call ptinpoly(ptlat, ptlon, srch_corners_locate_segstart,
4393  & srch_corner_lat_locate_segstart(:,ic),
4394  & srch_corner_lon_locate_segstart(:,ic),
4395  & latlon, inpoly, lboundary, edgeid)
4396 
4397  endif
4398 
4399  endif
4400 
4401  if (inpoly) then
4402  cont_cell = srch_cell_add
4403  exit
4404  endif
4405 
4406  end do
4407 
4408  return
4409 
4410 !----------------------------------------------------------------------
4411 

References first_call_locate_segstart, get_srch_cells(), scrip_grids::grid1_npole_cell, scrip_grids::grid1_spole_cell, scrip_grids::grid2_npole_cell, scrip_grids::grid2_spole_cell, scrip_constants::half, last_cell_grid_num_locate_segstart, last_cell_locate_segstart, last_srch_grid_num_locate_segstart, scrip_grids::north_thresh, scrip_grids::npseg, num_srch_cells_locate_segstart, scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, scrip_constants::pih, ptinpolarpoly(), ptinpoly(), ptinpolygen2(), scrip_constants::quart, scrip_grids::south_thresh, srch_add_locate_segstart, srch_center_lat_locate_segstart, srch_center_lon_locate_segstart, srch_corner_lat_locate_segstart, srch_corner_lon_locate_segstart, srch_corners_locate_segstart, scrip_constants::two, and scrip_constants::zero.

◆ modify_polar_cell()

subroutine scrip_remap_conservative::modify_polar_cell ( integer (scrip_i4), intent(inout)  ncorners,
integer (scrip_i4), intent(in)  nalloc,
real (scrip_r8), dimension(nalloc), intent(inout)  cell_corner_lat,
real (scrip_r8), dimension(nalloc), intent(inout)  cell_corner_lon 
)

Definition at line 2082 of file scrip_remap_conservative.f.

2082 
2083  !*** Input variables
2084 
2085  integer (SCRIP_i4), intent(in) ::
2086  & nalloc
2087 
2088  !*** In/Out Variables
2089 
2090  integer (SCRIP_i4), intent(inout) ::
2091  & ncorners
2092  real (SCRIP_r8), dimension(:), intent(inout) ::
2093  & cell_corner_lat(nalloc),
2094  & cell_corner_lon(nalloc)
2095 
2096  !*** Local variables
2097 
2098  integer (SCRIP_i4) ::
2099  & npcorners, ! Number of polar corners
2100  & pcorner, ! Index of the polar corner
2101  ! (if only 1 is found)
2102  & corner, ! Corner iterator variable
2103  & previdx, ! Index of previous corner to polar corner
2104  & nextidx ! Index of next corner to polar corner
2105 
2106  real (SCRIP_r8) ::
2107  & pole_lat, ! Latitude considered to be pole
2108  & prevlon, ! Latitude of previous corner to polar corner
2109  & nextlon ! Latitude of next corner to polar corner
2110 
2111 
2112  !***
2113  !*** Modify special cell with only one corner at the pole. Such
2114  !*** cells can have an artificially extreme distortion of the
2115  !*** edges when mapped to the Lambert projection because of the
2116  !*** span of longitudes on the edges So we will augment such
2117  !*** cells with degenerate edges at the pole so that the edges
2118  !*** coming off the pole will actually have the same longitude
2119  !*** values at both ends
2120  !***
2121  !*** lon_p lon_p+ lon_p lon_p-
2122  !*** pi/2 pi/2 pi/2 pi/2
2123  !*** * *--------*------*
2124  !*** / \ | |
2125  !*** / \ | |
2126  !*** / \ | |
2127  !*** * * * *
2128  !*** lon_p+ lon_p- lon_p+ lon_p-
2129  !*** lat_p+ lat_p- lat_p+ lat_p-
2130  !***
2131 
2132 
2133  !***
2134  !*** MAJOR ASSUMPTION HERE IS THAT CELL_CORNER_LAT AND
2135  !*** CELL_CORNER_LON HAVE ROOM TO GROW
2136  !***
2137  if (ncorners .ge. nalloc) return ! ** * No room to grow
2138 
2139  pcorner = 0
2140  npcorners = 0
2141  do corner = 1, ncorners
2142  if (abs(abs(cell_corner_lat(corner))-pih) .le. 1.0e-05) then
2143  pcorner = corner
2144  pole_lat = cell_corner_lat(corner)
2145  npcorners = npcorners + 1
2146  endif
2147  enddo
2148 
2149 
2150  if (npcorners .ne. 1) return !*** Not the kind of cell we want
2151 
2152  previdx = mod((pcorner-1)-1+ncorners,ncorners) + 1
2153  prevlon = cell_corner_lon(previdx)
2154 
2155  nextidx = mod(pcorner,ncorners) + 1
2156  nextlon = cell_corner_lon(nextidx)
2157 
2158  !*** Move entries from pcorner+1 on back by one
2159 
2160  do corner = ncorners, pcorner+1, -1
2161  cell_corner_lat(corner+1) = cell_corner_lat(corner)
2162  cell_corner_lon(corner+1) = cell_corner_lon(corner)
2163  enddo
2164 
2165  !*** Add a corner after pcorner
2166 
2167  cell_corner_lat(pcorner+1) = pole_lat
2168  cell_corner_lon(pcorner+1) = nextlon
2169 
2170  ncorners = ncorners+1
2171 
2172  !*** Move entries from pcorner on back by one
2173 
2174  do corner = ncorners, pcorner, -1
2175  cell_corner_lat(corner+1) = cell_corner_lat(corner)
2176  cell_corner_lon(corner+1) = cell_corner_lon(corner)
2177  enddo
2178 
2179  !*** Add a corner before pcorner
2180 
2181  cell_corner_lat(pcorner) = pole_lat
2182  cell_corner_lon(pcorner) = prevlon
2183 
2184  ncorners = ncorners+1
2185 

References scrip_constants::pih.

Referenced by cell_integrate(), and locate_point().

◆ pole_intersection()

subroutine scrip_remap_conservative::pole_intersection ( integer (scrip_i4), intent(in)  location,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corners_lat,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corners_lon,
integer (scrip_i4), intent(in)  cell_grid_id,
real (scrip_r8), intent(in)  beglat,
real (scrip_r8), intent(in)  beglon,
real (scrip_r8), intent(in)  endlat,
real (scrip_r8), intent(in)  endlon,
real (scrip_r8), dimension(2), intent(inout)  begseg,
integer (scrip_i4), intent(in)  begedge,
integer (scrip_i4), intent(out)  intedge,
real (scrip_r8), intent(out)  intrsct_lat,
real (scrip_r8), intent(out)  intrsct_lon,
real (scrip_r8), intent(out)  sinang2,
logical (scrip_logical), intent(out)  lcoinc,
logical (scrip_logical), intent(inout)  lthresh 
)

Definition at line 2728 of file scrip_remap_conservative.f.

2728 
2729 !-----------------------------------------------------------------------
2730 !
2731 ! Special intersection routine for line segment in cell close to
2732 ! poles
2733 ! A coordinate transformation (using a Lambert azimuthal
2734 ! equivalent projection) is performed to perform the intersection
2735 ! Also, since a straight line in lat-lon space is a curve in this
2736 ! transformed space, we represent each edge of the cell as having
2737 ! 'npseg' segments whose endpoints are mapped using the Lambert
2738 ! projection
2739 !
2740 !-----------------------------------------------------------------------
2741 
2742 !-----------------------------------------------------------------------
2743 !
2744 ! intent(in):
2745 !
2746 !-----------------------------------------------------------------------
2747 
2748  integer (SCRIP_i4), intent(in) ::
2749  & location ! cell to intersect segment with
2750 
2751  integer (SCRIP_i4), intent(in) ::
2752  & ncorners ! Number of cell corners
2753 
2754  real (SCRIP_r8), dimension(ncorners), intent(in) ::
2755  & cell_corners_lat, ! Cell corner coordinates
2756  & cell_corners_lon
2757 
2758  integer (SCRIP_i4), intent(in) ::
2759  & cell_grid_id ! which grid is the cell from?
2760 
2761  real (SCRIP_r8), intent(in) ::
2762  & beglat, beglon, ! beginning lat/lon coords for segment
2763  & endlat, endlon ! ending lat/lon coords for segment
2764 
2765  real (SCRIP_r8), dimension(2), intent(inout) ::
2766  & begseg ! begin lat/lon of full segment
2767 
2768  integer (SCRIP_i4), intent(in) ::
2769  & begedge ! edge on which segment start is on
2770  ! (can be 0)
2771 
2772 !-----------------------------------------------------------------------
2773 !
2774 ! intent(out):
2775 !
2776 !-----------------------------------------------------------------------
2777 
2778  integer (SCRIP_i4), intent(out) ::
2779  & intedge ! Edge that segment intersects
2780 
2781  real (SCRIP_r8), intent(out) ::
2782  & intrsct_lat, ! lat/lon coords of intersection
2783  & intrsct_lon
2784 
2785  real (SCRIP_r8), intent(out) ::
2786  & sinang2 ! square of sine of angle between
2787  ! intersecting line segments
2788 
2789  logical (SCRIP_logical), intent(out) ::
2790  & lcoinc ! True if segment is coincident with
2791  ! a cell edge
2792 
2793  logical (SCRIP_logical), intent(inout) ::
2794  & lthresh ! True if segment crosses threshold
2795 
2796 
2797 !-----------------------------------------------------------------------
2798 !
2799 ! local variables
2800 !
2801 !-----------------------------------------------------------------------
2802 
2803  integer (SCRIP_i4) ::
2804  & n, n1, next_n, prev_n,
2805  & it, i, j,
2806  & ncorners2,
2807  & intedge1
2808 
2809  logical (SCRIP_logical) ::
2810  & first,
2811  & found
2812 
2813  real (SCRIP_r8) ::
2814  & pi4, rns, ! north/south conversion
2815  & x1, x2, ! local x variables for segment
2816  & y1, y2, ! local y variables for segment
2817  & grdx1, grdx2, ! local x variables for grid cell
2818  & grdy1, grdy2, ! local y variables for grid cell
2819  & grdlat1, grdlat2, ! latitude vars for grid cell
2820  & grdlon1, grdlon2, ! longitude vars for grid cell
2821  & vec1_y, vec1_x, !
2822  & vec2_y, vec2_x, ! vectors and cross products used
2823  & vec3_y, vec3_x, !
2824  & vec1_lat, vec1_lon, !
2825  & vec2_lat, vec2_lon, !
2826  & vec3_lon, !
2827  & cross_product, !
2828  & dot_product, !
2829  & s1, s2, determ, ! variables used for linear solve to
2830  & mat1, mat2, !
2831  & mat3, mat4, ! find intersection
2832  & rhs1, rhs2, !
2833  & denom, !
2834  & intrsct_x, intrsct_y, ! intersection coordinates in
2835  ! transformed space
2836  & max_intrsct_lat, ! intersection point at max distance
2837  & max_intrsct_lon, ! from the start point
2838  & dist2, ! dist of intersection point from start
2839  ! point
2840  & maxdist2, ! max dist of intersection point from
2841  ! start pnt
2842  & lensqr1, lensqr2, ! various segment lengths
2843  & lensqr3,
2844  & tmpx, tmpy,
2845  & tmplat, tmplon,
2846  & ldummy
2847 
2848  !***
2849  !*** variables necessary if segment manages to hit pole
2850  !***
2851 
2852  real (SCRIP_r8), dimension(npseg*ncorners) ::
2853  & cell_corners_lat_loc,! Lat/Lon coordinates of multi-segmented
2854  & cell_corners_lon_loc ! version of cell
2855 
2856 
2857 
2858 !-----------------------------------------------------------------------
2859 !
2860 ! initialize defaults, flags, etc.
2861 !
2862 !-----------------------------------------------------------------------
2863 
2864  max_intrsct_lat = pi ! intersection point at max distance
2865  max_intrsct_lon = 4*pi ! from the start point
2866 
2867  intedge = 0
2868  first = .true.
2869  maxdist2 = -999999.00
2870 
2871  s1 = zero
2872 
2873 !-----------------------------------------------------------------------
2874 !
2875 ! convert coordinates
2876 !
2877 !-----------------------------------------------------------------------
2878 
2879  if (beglat > zero) then
2880  pi4 = quart*pi
2881  rns = one
2882  else
2883  pi4 = -quart*pi
2884  rns = -one
2885  endif
2886 
2887  x1 = rns*two*sin(pi4 - half*beglat)*cos(beglon)
2888  y1 = two*sin(pi4 - half*beglat)*sin(beglon)
2889  x2 = rns*two*sin(pi4 - half*endlat)*cos(endlon)
2890  y2 = two*sin(pi4 - half*endlat)*sin(endlon)
2891 
2892  intrsct_x = x2
2893  intrsct_y = y2
2894 
2895 
2896 !-----------------------------------------------------------------------
2897 !
2898 ! now that a cell is found, search for the next intersection.
2899 ! loop over sides of the cell to find intersection with side
2900 ! must check all sides for coincidences or intersections
2901 !
2902 !-----------------------------------------------------------------------
2903 
2904 
2905  if (abs(x1) .le. tiny .and. abs(y1) .le. tiny .and.
2906  & abs(x2) .le. tiny .and. abs(y2) .le. tiny) then
2907 
2908  !***
2909  !*** The segment is a polar segment which is degenerate
2910  !*** in the transformed Lambert space. Find out which
2911  !*** cell edge it is coincident with and find the
2912  !*** point where the segment exits this cell (if at all)
2913  !*** NOTE 1: THIS MUST BE DONE IN LAT-LON SPACE
2914  !*** NOTE 2: CODE RELEVANT ONLY FOR INTEGRATION W.R.T. phi
2915  !***
2916 
2917  intrsct_loop1: do n = 1, ncorners
2918  next_n = mod(n,ncorners) + 1
2919 
2920  grdlat1 = cell_corners_lat(n)
2921  grdlon1 = cell_corners_lon(n)
2922  grdlat2 = cell_corners_lat(next_n)
2923  grdlon2 = cell_corners_lon(next_n)
2924  grdx1 = rns*two*sin(pi4 - half*grdlat1)*cos(grdlon1)
2925  grdy1 = two*sin(pi4 - half*grdlat1)*sin(grdlon1)
2926  grdx2 = rns*two*sin(pi4 - half*grdlat2)*cos(grdlon2)
2927  grdy2 = two*sin(pi4 - half*grdlat2)*sin(grdlon2)
2928 
2929  if (abs(grdx1) .le. tiny .and. abs(grdy1) .le. tiny .and.
2930  & abs(grdx2) .le. tiny .and. abs(grdy2) .le. tiny) then
2931 
2932  !***
2933  !*** Found polar segment in cell
2934  !***
2935 
2936  vec1_lon = endlon-beglon
2937  if (vec1_lon .gt. pi) vec1_lon = vec1_lon - pi2
2938  if (vec1_lon .lt. -pi) vec1_lon = vec1_lon + pi2
2939 
2940  vec2_lon = grdlon2-grdlon1
2941  if (vec2_lon .gt. pi) vec2_lon = vec2_lon - pi2
2942  if (vec2_lon .lt. -pi) vec2_lon = vec2_lon + pi2
2943 
2944  if (vec1_lon*vec2_lon .lt. 0) then
2945 
2946  !*** switch coordinates to simplify logic below
2947 
2948  tmplat = grdlat2
2949  tmplon = grdlon2
2950  grdlat2 = grdlat1
2951  grdlon2 = grdlon1
2952  grdlat1 = tmplat
2953  grdlon1 = tmplon
2954  endif
2955 
2956  vec2_lon = grdlon1 - beglon
2957  if (vec2_lon .gt. pi) vec2_lon = vec2_lon - pi2
2958  if (vec2_lon .lt. -pi) vec2_lon = vec2_lon + pi2
2959 
2960  vec3_lon = grdlon2 - beglon
2961  if (vec3_lon .gt. pi) vec3_lon = vec3_lon - pi2
2962  if (vec3_lon .lt. -pi) vec3_lon = vec3_lon + pi2
2963 
2964  found = .false.
2965 
2966  if (vec2_lon*vec1_lon > 0) then
2967  if (abs(vec3_lon) < abs(vec1_lon)) then
2968  intrsct_lon = grdlon2
2969  found = .true.
2970  else if (abs(vec2_lon) < abs(vec1_lon)) then
2971  intrsct_lon = grdlon1 ! Shouldn't be here
2972  found = .true.
2973  endif
2974  else
2975  if (vec3_lon*vec1_lon > 0) then
2976  if (abs(vec3_lon) < abs(vec1_lon)) then
2977  intrsct_lon = grdlon2
2978  found = .true.
2979  endif
2980  endif
2981 
2982  endif
2983 
2984  if (found) then
2985  intrsct_lat = endlat
2986  lcoinc = .true.
2987  sinang2 = 0.0
2988  intedge = n
2989  exit intrsct_loop1
2990  endif
2991 
2992  endif
2993 
2994  end do intrsct_loop1
2995 
2996  return
2997  endif
2998 
2999 
3000 
3001 
3002  !****
3003  !**** General intersection
3004  !****
3005 
3006 
3007 
3008  !***
3009  !*** Construct multi-segmented version of the cell
3010  !***
3011 
3012  i = 0
3013  do n = ncorners, 1, -1
3014  i = i+1
3015  n1 = mod(n,ncorners)+1
3016  cell_corners_lat_loc(i) = cell_corners_lat(n1)
3017  cell_corners_lon_loc(i) = cell_corners_lon(n1)
3018 
3019  prev_n = n1-1
3020  if (prev_n .eq. 0) prev_n = ncorners ! how do we do (j-1+n)%n
3021  ! in F90 ?
3022 
3023  vec1_lat = cell_corners_lat(prev_n)-cell_corners_lat(n1)
3024  vec1_lon = cell_corners_lon(prev_n)-cell_corners_lon(n1)
3025  if (vec1_lon > pi) then
3026  vec1_lon = vec1_lon - pi2
3027  else if (vec1_lon < -pi) then
3028  vec1_lon = vec1_lon + pi2
3029  endif
3030 
3031  do j = 1, npseg-1
3032  i = i+1
3033  cell_corners_lat_loc(i) =
3034  & cell_corners_lat(n1) + j*vec1_lat/npseg
3035  cell_corners_lon_loc(i) =
3036  & cell_corners_lon(n1) + j*vec1_lon/npseg
3037  enddo
3038  enddo
3039 
3040  ncorners2 = npseg*ncorners
3041 
3042 
3043 
3044  !***
3045  !*** Now intersect segment with multi-segmented version of cell
3046  !***
3047 
3048 
3049  intrsct_loop2: do n= 1, ncorners2
3050 
3051  next_n = mod(n,ncorners2) + 1
3052 
3053  grdlat1 = cell_corners_lat_loc(n)
3054  grdlon1 = cell_corners_lon_loc(n)
3055  grdlat2 = cell_corners_lat_loc(next_n)
3056  grdlon2 = cell_corners_lon_loc(next_n)
3057  grdx1 = rns*two*sin(pi4 - half*grdlat1)*cos(grdlon1)
3058  grdy1 = two*sin(pi4 - half*grdlat1)*sin(grdlon1)
3059  grdx2 = rns*two*sin(pi4 - half*grdlat2)*cos(grdlon2)
3060  grdy2 = two*sin(pi4 - half*grdlat2)*sin(grdlon2)
3061 
3062  if ((grdx1-grdx2)*(grdx1-grdx2)+(grdy1-grdy2)*(grdy1-grdy2) .le.
3063  & tiny*tiny) cycle
3064 
3065 
3066  !***
3067  !*** set up linear system to solve for intersection
3068  !***
3069 
3070  mat1 = x2 - x1
3071  mat2 = grdx1 - grdx2
3072  mat3 = y2 - y1
3073  mat4 = grdy1 - grdy2
3074  rhs1 = grdx1 - x1
3075  rhs2 = grdy1 - y1
3076 
3077  determ = mat1*mat4 - mat2*mat3
3078 
3079  !***
3080  !*** if the determinant is zero, the segments are either
3081  !*** parallel or coincident or one segment has zero length.
3082 
3083  !*** if the determinant is non-zero, solve for the linear
3084  !*** parameters s for the intersection point on each line
3085  !*** segment.
3086  !*** if 0<s1,s2<1 then the segment intersects with this side.
3087  !*** return the point of intersection (adding a small
3088  !*** number so the intersection is off the grid line).
3089  !***
3090 
3091  if (abs(determ) > 1.e-30) then
3092 
3093  s1 = (rhs1*mat4 - mat2*rhs2)/determ
3094  s2 = (mat1*rhs2 - rhs1*mat3)/determ
3095 
3096  if (s2 >= zero .and. s2 <= one .and.
3097  & s1 > tiny .and. s1 <= one) then
3098 
3099  intrsct_x = x1 + s1*mat1
3100  intrsct_y = y1 + s1*mat3
3101 
3102  !***
3103  !*** convert back to lat/lon coordinates
3104  !***
3105 
3106  if (abs(intrsct_x) .gt. tiny .or.
3107  & abs(intrsct_y) .gt. tiny) then
3108 
3109  intrsct_lon = rns*atan2(intrsct_y,intrsct_x)
3110 
3111  else
3112 
3113  !*** Degenerate case - we don't have a good way of
3114  !*** finding out what the longitude corresponding
3115  !*** to a (0,0) intersection is. So we take the
3116  !*** the intersection as one of the two endpoints of
3117  !*** the grid segment
3118 
3119  if (abs(abs(grdlat1)-pih) .lt. 1e-5 .and.
3120  & abs(abs(grdlat2)-pih) .lt. 1e-5) then
3121 
3122  !*** Both endpoints of the grid segment are at the pole
3123  !*** but at different longitudes
3124 
3125  vec1_lat = grdlat1-beglat
3126  vec1_lon = grdlon1-beglon
3127  if (vec1_lon > pi) then
3128  vec1_lon = vec1_lon - pi2
3129  else if (vec1_lon < -pi) then
3130  vec1_lon = vec1_lon + pi2
3131  endif
3132  dist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
3133 
3134  vec2_lat = grdlat2-beglat
3135  vec2_lon = grdlon2-beglon
3136  if (vec2_lon > pi) then
3137  vec2_lon = vec2_lon - pi2
3138  else if (vec2_lon < -pi) then
3139  vec2_lon = vec2_lon + pi2
3140  endif
3141 
3142  !*** pick the endpoint of the grid segment that is
3143  !*** farthest from the beg point of the segment
3144 
3145  if ((vec1_lat*vec1_lat + vec1_lon*vec1_lon) .ge.
3146  & (vec2_lat*vec2_lat + vec2_lon*vec2_lon)) then
3147  intrsct_lon = grdlon1
3148  else
3149  intrsct_lon = grdlon2
3150  endif
3151 
3152  else if (abs(abs(grdlat1)-pih) .lt. 1e-5) then
3153  intrsct_lon = grdlon1
3154  else if (abs(abs(grdlat2)-pih) .lt. 1e-5) then
3155  intrsct_lon = grdlon2
3156  endif
3157 
3158  !*** Make sure this longitude is not outside the
3159  !*** beglon,endlon range
3160 
3161  vec1_lon = endlon-intrsct_lon
3162  if (vec1_lon > pi) then
3163  vec1_lon = vec1_lon - pi2
3164  else if (vec1_lon < -pi) then
3165  vec1_lon = vec1_lon + pi2
3166  endif
3167 
3168  vec2_lon = beglon-intrsct_lon
3169  if (vec2_lon > pi) then
3170  vec2_lon = vec2_lon - pi2
3171  else if (vec2_lon < -pi) then
3172  vec2_lon = vec2_lon + pi2
3173  endif
3174 
3175  !*** if vec1_lon and vec2_lon are of the same sign
3176  !*** then intrsct_lon is outside the beglon,endlon
3177  !*** range
3178 
3179  if (vec1_lon*vec2_lon > 0) cycle
3180 
3181  endif
3182 
3183  if (intrsct_lon < zero)
3184  & intrsct_lon = intrsct_lon + pi2
3185 
3186  if (abs(intrsct_x) > 1.d-10) then
3187  intrsct_lat = (pi4 -
3188  & asin(rns*half*intrsct_x/cos(intrsct_lon)))*two
3189  ldummy = two*(pi4 -
3190  & asin(sqrt(intrsct_x*intrsct_x+intrsct_y*intrsct_y)/2.))
3191  else if (abs(intrsct_y) > 1.d-10) then
3192  intrsct_lat = (pi4 -
3193  & asin(half*intrsct_y/sin(intrsct_lon)))*two
3194  ldummy = two*(pi4 -
3195  & asin(sqrt(intrsct_x*intrsct_x+intrsct_y*intrsct_y)/2.))
3196  else
3197  intrsct_lat = two*pi4
3198  endif
3199 
3200 
3201  !***
3202  !*** If there are multiple intersection points, accept the
3203  !*** one that is not on the edge we started from but is
3204  !*** closest to the start point - need this for
3205  !*** intersection to work for non-convex edges
3206  !***
3207 
3208  if (first) then
3209 
3210  intedge1 = (n-1)/npseg + 1
3211  intedge1 = ncorners - intedge1 + 1 ! dir of edges was
3212  ! reversed
3213  if (intedge1 .ne. begedge) then
3214 
3215  max_intrsct_lat = intrsct_lat
3216  max_intrsct_lon = intrsct_lon
3217 
3218  vec1_lat = intrsct_lat-beglat
3219  vec1_lon = intrsct_lon-beglon
3220  if (vec1_lon > pi) then
3221  vec1_lon = vec1_lon - pi2
3222  else if (vec1_lon < -pi) then
3223  vec1_lon = vec1_lon + pi2
3224  endif
3225  maxdist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
3226  dist2 = maxdist2
3227 
3228  denom = (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
3229  sinang2 = determ*determ/denom
3230  intedge = intedge1
3231 
3232  first = .false.
3233  endif
3234 
3235  else
3236  vec1_lat = intrsct_lat-beglat
3237  vec1_lon = intrsct_lon-beglon
3238  if (vec1_lon > pi) then
3239  vec1_lon = vec1_lon - pi2
3240  else if (vec1_lon < -pi) then
3241  vec1_lon = vec1_lon + pi2
3242  endif
3243  dist2 = vec1_lat*vec1_lat + vec1_lon*vec1_lon
3244 
3245  !*** if the first intersection was on the same edge
3246  !*** as the starting edge or
3247  !*** the current intersection point is not on the
3248  !*** starting edge and the distance to the beginning
3249  !*** point is less than that of the previous
3250  !*** intersection accept this intersection
3251 
3252  intedge1 = (n-1)/npseg + 1
3253  intedge1 = ncorners - intedge1 + 1 ! dir of edges was
3254  ! reversed
3255  if (dist2 > maxdist2) then
3256  if (begedge == 0 .or. intedge1 .ne. begedge) then
3257  max_intrsct_lat = intrsct_lat
3258  max_intrsct_lon = intrsct_lon
3259  maxdist2 = dist2
3260 
3261  denom =
3262  & (mat1*mat1+mat2*mat2)*(mat3*mat3+mat4*mat4)
3263  sinang2 = determ*determ/denom
3264  intedge = intedge1
3265  endif
3266  endif
3267  endif
3268  endif
3269 
3270  else
3271 
3272  !***
3273  !*** Coincident lines or parallel lines
3274  !***
3275 
3276  cross_product = mat2*rhs2 - mat4*rhs1
3277 
3278  if (abs(cross_product) < tiny) then
3279 
3280  dot_product = mat1*(-mat2) + mat3*(-mat4)
3281 
3282  !***
3283  !*** If area of triangle formed by x2,y2 and the gridline
3284  !*** is negligible then the lines are coincident
3285  !***
3286 
3287  lensqr1 = mat1*mat1 + mat3*mat3 ! length sqrd of input
3288  ! segment
3289 
3290  if (dot_product < zero) then
3291  tmpx = grdx2
3292  tmpy = grdy2
3293  tmplat = grdlat2
3294  tmplon = grdlon2
3295  grdx2 = grdx1
3296  grdy2 = grdy1
3297  grdlat2 = grdlat1
3298  grdlon2 = grdlon1
3299  grdx1 = tmpx
3300  grdy1 = tmpy
3301  grdlat1 = tmplat
3302  grdlon1 = tmplon
3303  endif
3304 
3305 
3306  vec2_x = grdx1 - x1
3307  vec2_y = grdy1 - y1
3308  lensqr2 = vec2_x*vec2_x + vec2_y*vec2_y
3309  if (vec2_x*mat1+vec2_y*mat3 < 0) then
3310  lensqr2 = -lensqr2
3311  endif
3312 
3313 
3314  vec3_x = grdx2 - x1
3315  vec3_y = grdy2 - y1
3316  lensqr3 = (vec3_x*vec3_x+vec3_y*vec3_y)
3317  if (vec3_x*mat1+vec3_y*mat3 < 0) then
3318  lensqr3 = -lensqr3
3319  endif
3320 
3321  found = .false.
3322 
3323  if (lensqr2 > 0) then
3324  if (lensqr2 <= lensqr1) then
3325  intrsct_x = grdx1
3326  intrsct_y = grdy1
3327  intrsct_lat = grdlat1
3328  intrsct_lon = grdlon1
3329  found = .true.
3330  endif
3331  else
3332  if (lensqr3 > 0) then
3333  if (lensqr3 > lensqr1) then
3334  intrsct_x = x2
3335  intrsct_y = y2
3336  intrsct_lat = endlat
3337  intrsct_lon = endlon
3338  found = .true.
3339  else
3340  intrsct_x = grdx2
3341  intrsct_y = grdy2
3342  intrsct_lat = grdlat2
3343  intrsct_lon = grdlon2
3344  found = .true.
3345  endif
3346  endif
3347  endif
3348 
3349  if (found) then
3350  dist2 = (intrsct_lat-beglat)*(intrsct_lat-beglat)+
3351  & (intrsct_lon-beglon)*(intrsct_lon-beglon)
3352 
3353  if (dist2 > tiny*tiny) then
3354 
3355  !*** Coincidence intersection always wins
3356 
3357  max_intrsct_lat = intrsct_lat
3358  max_intrsct_lon = intrsct_lon
3359  maxdist2 = dist2
3360  sinang2 = 0
3361  intedge = (n-1)/npseg + 1
3362  intedge = ncorners - intedge + 1
3363  lcoinc = .true.
3364 
3365  exit intrsct_loop2
3366  endif
3367  endif
3368 
3369 
3370  endif ! if (abs(cross_product) < tiny)
3371 
3372  endif ! if (abs(determ) > 1.e-30) .. else .. endif
3373 
3374  end do intrsct_loop2
3375 
3376  if (maxdist2 < 1e6*tiny*tiny) then
3377  intedge = 0
3378  return
3379  else
3380  intrsct_lat = max_intrsct_lat
3381  intrsct_lon = max_intrsct_lon
3382  endif
3383 
3384 !-----------------------------------------------------------------------
3385 !
3386 ! if segment manages to cross over pole, shift the beginning
3387 ! endpoint in order to avoid hitting pole directly
3388 ! (it is ok for endpoint to be pole point)
3389 !
3390 !-----------------------------------------------------------------------
3391 
3392  if (abs(intrsct_x) < 1.e-10 .and. abs(intrsct_y) < 1.e-10 .and.
3393  & (x2 /= zero .and. y2 /=0)) then
3394  if (avoid_pole_count > 2) then
3395  avoid_pole_count = 0
3396  avoid_pole_offset = 10.*avoid_pole_offset
3397  endif
3398 
3399  cross_product = x1*(y2-y1) - y1*(x2-x1)
3400  intrsct_lat = beglat
3401  if (cross_product*intrsct_lat > zero) then
3402  intrsct_lon = beglon + avoid_pole_offset
3403  else
3404  intrsct_lon = beglon - avoid_pole_offset
3405  endif
3406 
3407  avoid_pole_count = avoid_pole_count + 1
3408  else
3409  avoid_pole_count = 0
3410  avoid_pole_offset = tiny
3411  endif
3412 
3413 !-----------------------------------------------------------------------
3414 !
3415 ! if the segment crosses a pole threshold, reset the intersection
3416 ! to be the threshold latitude and do not reuse x,y intersect
3417 ! on next entry. only check if did not cross threshold last
3418 ! time - sometimes the coordinate transformation can place a
3419 ! segment on the other side of the threshold again
3420 !
3421 !-----------------------------------------------------------------------
3422 
3423  if (lthresh) then
3424  if (intrsct_lat > north_thresh .or. intrsct_lat < south_thresh)
3425  & lthresh = .false.
3426  else if (beglat > zero .and. intrsct_lat < north_thresh) then
3427  mat4 = endlat - begseg(1)
3428  mat3 = endlon - begseg(2)
3429  if (mat3 > pi) mat3 = mat3 - pi2
3430  if (mat3 < -pi) mat3 = mat3 + pi2
3431 ! intrsct_lat = north_thresh - tiny
3432  intrsct_lat = north_thresh
3433  s1 = (north_thresh - begseg(1))/mat4
3434  intrsct_lon = begseg(2) + s1*mat3
3435  lthresh = .true.
3436  else if (beglat < zero .and. intrsct_lat > south_thresh) then
3437  mat4 = endlat - begseg(1)
3438  mat3 = endlon - begseg(2)
3439  if (mat3 > pi) mat3 = mat3 - pi2
3440  if (mat3 < -pi) mat3 = mat3 + pi2
3441 ! intrsct_lat = south_thresh + tiny
3442  intrsct_lat = south_thresh
3443  s1 = (south_thresh - begseg(1))/mat4
3444  intrsct_lon = begseg(2) + s1*mat3
3445  lthresh = .true.
3446  endif
3447 
3448 
3449 !-----------------------------------------------------------------------
3450 

References avoid_pole_count, avoid_pole_offset, scrip_constants::half, scrip_grids::north_thresh, scrip_grids::npseg, scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, scrip_constants::pih, scrip_constants::quart, scrip_grids::south_thresh, scrip_constants::tiny, scrip_constants::two, and scrip_constants::zero.

Referenced by intersection().

◆ ptincell()

subroutine scrip_remap_conservative::ptincell ( real (scrip_r8), intent(in)  ptlat,
real (scrip_r8), intent(in)  ptlon,
integer (scrip_i4), intent(in)  cell_add,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_lat,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_lon,
real (scrip_r8), intent(in)  cell_center_lat,
real (scrip_r8), intent(in)  cell_center_lon,
integer (scrip_i4), intent(in)  cell_grid_id,
logical (scrip_logical), intent(out)  inpoly,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 4608 of file scrip_remap_conservative.f.

4608 
4609 !----------------------------------------------------------------------
4610 
4611  implicit none
4612 
4613 !-----------------------------------------------------------------------
4614 !
4615 ! intent(in):
4616 !
4617 !-----------------------------------------------------------------------
4618 
4619  real (SCRIP_r8), intent(in) ::
4620  & ptlat, ptlon ! Point to locate
4621 
4622  integer (SCRIP_i4), intent(in) ::
4623  & cell_add ! ID of cell
4624 
4625  integer (SCRIP_i4), intent(in) ::
4626  & ncorners
4627 
4628  real (SCRIP_r8), dimension(ncorners), intent(in) ::
4629  & cell_corner_lat, cell_corner_lon
4630 
4631  real (SCRIP_r8), intent(in) ::
4632  & cell_center_lat,
4633  & cell_center_lon
4634 
4635  integer (SCRIP_i4), intent(in) ::
4636  & cell_grid_id ! num indicating if we are locating a grid1
4637  ! point in a cell of grid2 (num = 2) or
4638  ! a grid2 point in a cell of grid1 (num = 1)
4639 
4640 
4641 !-----------------------------------------------------------------------
4642 !
4643 ! intent(out):
4644 !
4645 !-----------------------------------------------------------------------
4646 
4647  logical (SCRIP_logical), intent(out) ::
4648  & inpoly ! is point in polygon?
4649 
4650  logical (SCRIP_logical), intent(out) ::
4651  & lboundary ! flag points that lie on the boundary
4652  ! of the cell
4653 
4654  integer (SCRIP_i4), intent(out) ::
4655  & edgeid ! if point is on boundary, which local
4656  ! edge is it on? (0 otherwise)
4657 !-----------------------------------------------------------------------
4658 !
4659 ! local variables
4660 !
4661 !-----------------------------------------------------------------------
4662 
4663  integer (SCRIP_i4) :: i, j, k, ic
4664  integer (SCRIP_i4) :: whichpole
4665 
4666  real (SCRIP_r8) :: rns, pi4, ptx, pty, lat, lon,
4667  & cell_center_x, cell_center_y, vec1_lat, vec1_lon
4668 
4669  logical (SCRIP_logical) ::
4670  & latlon
4671 
4672  real (kind=scrip_r8), dimension(npseg*ncorners) ::
4673  & cell_corner_x, ! x of each corner of cell
4674  & cell_corner_y ! y of each corner of cell
4675 
4676 !----------------------------------------------------------------------
4677 
4678  edgeid = 0
4679 
4680 
4681  !*** IF POINTS ARE ABOVE THE THRESHOLD, CHECK THEM IN A TRANSFORMED
4682  !*** SPACE
4683  !*** HOWEVER, POINTS THAT ARE PRACTICALLY AT THE POLE CANNOT
4684  !*** BE CORRECTLY LOCATED THIS WAY BECAUSE THE POLE IS A SINGULARITY
4685  !*** AND CONTAINMENT IN ANY CELL INCIDENT ON THE POLE WILL GIVE US A
4686  !*** POSITIVE ANSWER. FOR THESE POINTS REVERT TO THE LATLON SPACE
4687  !***
4688 
4689  if ((ptlat .gt. north_thresh .and. abs(ptlat-pih) .ge. 0.001) .or.
4690  & (ptlat .lt. south_thresh .and. abs(ptlat+pih) .ge. 0.001))
4691  & then
4692 
4693  if (ptlat > zero) then
4694  pi4 = quart*pi
4695  rns = one
4696  else
4697  pi4 = -quart*pi
4698  rns = -one
4699  endif
4700 
4701  ptx = rns*two*sin(pi4 - half*ptlat)*cos(ptlon)
4702  pty = two*sin(pi4 - half*ptlat)*sin(ptlon)
4703 
4704  latlon = .false.
4705 
4706  ! Since we want greater fidelity for locating the points
4707  ! we send in the mid-points of the polygon edges too
4708  ! BUT THAT MAKES THE POLYGON NON-CONVEX SOMETIMES AND
4709  ! THE CROSS-PRODUCT CHECK FAILS. SO USE CODE TO CHECK GENERAL
4710  ! POLYGONS
4711 
4712 
4713  k = 0
4714  do i = ncorners, 1, -1
4715  k = k+1
4716  lat = cell_corner_lat(i)
4717  lon = cell_corner_lon(i)
4718  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4719  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4720 
4721  j = i-1
4722  if (j .eq. 0) j = ncorners ! how do we do (j-1+n)%n in F90?
4723 
4724  vec1_lat = cell_corner_lat(j)-cell_corner_lat(i)
4725  vec1_lon = cell_corner_lon(j)-cell_corner_lon(i)
4726  if (vec1_lon > pi) then
4727  vec1_lon = vec1_lon - pi2
4728  else if (vec1_lon < -pi) then
4729  vec1_lon = vec1_lon + pi2
4730  endif
4731 
4732  do j = 1, npseg-1
4733  k = k+1
4734  lat = cell_corner_lat(i) + j*vec1_lat/npseg
4735  lon = cell_corner_lon(i) + j*vec1_lon/npseg
4736  cell_corner_x(k) = rns*two*sin(pi4-half*lat)*cos(lon)
4737  cell_corner_y(k) = two*sin(pi4-half*lat)*sin(lon)
4738  enddo
4739  enddo
4740 
4741  !*** cell is so non-convex that no feasible center exists
4742  !*** we have to fall back on a different algorithm
4743 
4744  call ptinpolygen2(ptx, pty, k, cell_corner_x,
4745  & cell_corner_y, latlon, inpoly, lboundary, edgeid)
4746 
4747  if (lboundary) then
4748  edgeid = (edgeid-1)/npseg + 1 ! convert from index in
4749  ! multi-segmented cell to
4750  ! regular cell
4751  endif
4752  else
4753 
4754  latlon = .true.
4755 
4756  whichpole = 0
4757  if (cell_grid_id .eq. 1 .and.
4758  & cell_add .eq. grid1_spole_cell) then
4759 
4760  whichpole = -1 ! S pole
4761  call ptinpolarpoly(ptlat, ptlon, ncorners,
4762  & cell_corner_lat, cell_corner_lon,
4763  & latlon, whichpole, inpoly, lboundary, edgeid)
4764 
4765  else if (cell_grid_id .eq. 1 .and.
4766  & cell_add .eq. grid1_npole_cell) then
4767 
4768  whichpole = 1 ! N pole
4769  call ptinpolarpoly(ptlat, ptlon, ncorners,
4770  & cell_corner_lat, cell_corner_lon,
4771  & latlon, whichpole, inpoly, lboundary, edgeid)
4772 
4773  else if (cell_grid_id .eq. 2 .and.
4774  & cell_add .eq. grid2_spole_cell) then
4775 
4776  whichpole = -1 ! S pole
4777  call ptinpolarpoly(ptlat, ptlon, ncorners,
4778  & cell_corner_lat, cell_corner_lon,
4779  & latlon, whichpole, inpoly, lboundary, edgeid)
4780 
4781  else if (cell_grid_id .eq. 2 .and.
4782  & cell_add .eq. grid2_npole_cell) then
4783 
4784  whichpole = 1 ! N pole
4785  call ptinpolarpoly(ptlat, ptlon, ncorners,
4786  & cell_corner_lat, cell_corner_lon,
4787  & latlon, whichpole, inpoly, lboundary, edgeid)
4788 
4789  else
4790 
4791  !***
4792  !*** General cell
4793  !***
4794 
4795  call ptinpoly(ptlat, ptlon, ncorners,
4796  & cell_corner_lat, cell_corner_lon,
4797  & latlon, inpoly, lboundary, edgeid)
4798 
4799  endif
4800 
4801  endif
4802 
4803  return
4804 
4805 !----------------------------------------------------------------------
4806 

References scrip_grids::grid1_npole_cell, scrip_grids::grid1_spole_cell, scrip_grids::grid2_npole_cell, scrip_grids::grid2_spole_cell, scrip_constants::half, scrip_grids::north_thresh, scrip_grids::npseg, scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, scrip_constants::pih, ptinpolarpoly(), ptinpoly(), ptinpolygen2(), scrip_constants::quart, scrip_grids::south_thresh, scrip_constants::two, and scrip_constants::zero.

Referenced by cell_integrate(), converge_to_bdry(), and locate_point().

◆ ptinpolarpoly()

subroutine scrip_remap_conservative::ptinpolarpoly ( real (scrip_r8), intent(in)  ptx,
real (scrip_r8), intent(in)  pty,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_x,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_y,
logical (scrip_logical), intent(in)  latlon,
integer (scrip_i4), intent(in)  whichpole,
logical (scrip_logical), intent(out)  inpoly,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 5042 of file scrip_remap_conservative.f.

5042 
5043 !----------------------------------------------------------------------
5044 !
5045 ! Check if point is in polygonal cell overlapping the pole
5046 ! Cannot check the containment as is in latlon space - We have
5047 ! to check by connecting each edge of the polygon to the pole
5048 ! and check containment in the resulting quadrilateral in latlon
5049 ! space
5050 ! The cell can be non-convex as long as the pole is 'visible' to
5051 ! all the edges of the polygon, i.e., we can connect the pole to
5052 ! each edge of the polygon and form a triangle with positive area
5053 !
5054 !----------------------------------------------------------------------
5055 
5056  implicit none
5057 
5058 !----------------------------------------------------------------------
5059 !
5060 ! Input arguments
5061 !
5062 !----------------------------------------------------------------------
5063 
5064  real (SCRIP_r8), intent(in) ::
5065  & ptx, pty ! Point to check
5066 
5067  integer (SCRIP_i4), intent(in) ::
5068  & ncorners ! Number of polygon corners
5069 
5070  real (SCRIP_r8), dimension(ncorners), intent(in) ::
5071  & cell_corner_x, ! Coordinates of cell corners
5072  & cell_corner_y ! Could be x-y or lat-lon or ...
5073 
5074  logical (SCRIP_logical), intent(in) ::
5075  & latlon ! Are coordinates in latlon space?
5076 
5077  integer (SCRIP_i4), intent(in) ::
5078  & whichpole ! South or North pole
5079 
5080 !----------------------------------------------------------------------
5081 !
5082 ! Output arguments
5083 !
5084 !----------------------------------------------------------------------
5085 
5086  logical (SCRIP_logical), intent(out) ::
5087  & inpoly ! Is point in the polygon?
5088 
5089  logical (SCRIP_logical), intent(out) ::
5090  & lboundary ! Is point on the boundary of the polygon?
5091 
5092  integer (SCRIP_i4), intent(out) ::
5093  & edgeid ! if point is on boundary, which local
5094  ! edge is it on? (0 otherwise)
5095 
5096 !----------------------------------------------------------------------
5097 !
5098 ! Local variables
5099 !
5100 !----------------------------------------------------------------------
5101 
5102  integer (SCRIP_i4) :: n, next_n, ledgeid
5103 
5104  real (SCRIP_r8), dimension(4) ::
5105  & pquad_corner_x, ! Coordinates of polar quad
5106  & pquad_corner_y
5107 
5108  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
5109  & cross_product, pole_lat
5110 
5111  pole_lat = whichpole*pih
5112 
5113  !***
5114  !*** This is a polygon that overlaps the pole
5115  !*** A normal point in polygon check could fail
5116  !*** So, with each edge of the polygon form a quadrilateral
5117  !*** in latlon space using the polar latitude and the longitude
5118  !*** values of the endpoints of the edge. Then check containment
5119  !*** of the point in this quadrilateral
5120  !***
5121 
5122  inpoly = .false.
5123  lboundary = .false.
5124 
5125  do n = 1, ncorners
5126  next_n = mod(n,ncorners) + 1
5127 
5128  pquad_corner_x(1) = cell_corner_x(n)
5129  pquad_corner_y(1) = cell_corner_y(n)
5130  pquad_corner_x(2) = cell_corner_x(next_n)
5131  pquad_corner_y(2) = cell_corner_y(next_n)
5132  pquad_corner_x(3) = pole_lat
5133  pquad_corner_y(3) = cell_corner_y(next_n)
5134  pquad_corner_x(4) = pole_lat
5135  pquad_corner_y(4) = cell_corner_y(n)
5136 
5137 
5138  call ptinpoly(ptx,pty,4,pquad_corner_x,pquad_corner_y,
5139  & latlon,inpoly,lboundary, ledgeid)
5140 
5141  if (inpoly) then
5142 
5143  if (lboundary) then
5144 
5145  !***
5146  !*** Check to see if the lboundary flag is being
5147  !*** triggered by the outer edge of the polygon or
5148  !*** by one of the artificial internal edges
5149  !***
5150 
5151  vec1_x = pquad_corner_x(2) - pquad_corner_x(1)
5152  vec1_y = pquad_corner_y(2) - pquad_corner_y(1)
5153  vec2_x = ptx - pquad_corner_x(1)
5154  vec2_y = pty - pquad_corner_y(1)
5155 
5156 
5157  if (latlon) then
5158 
5159  !***
5160  !*** check for 0,2pi crossings
5161  !***
5162 
5163  if (vec1_y > pi) vec1_y = vec1_y - pi2
5164  if (vec1_y < -pi) vec1_y = vec1_y + pi2
5165  if (vec2_y > pi) vec2_y = vec2_y - pi2
5166  if (vec2_y < -pi) vec2_y = vec2_y + pi2
5167 
5168  endif
5169 
5170  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
5171 
5172  !***
5173  !*** if the cross product for a side is zero, the point
5174  !*** lies exactly on the side or the side is degenerate
5175  !*** (zero length). if degenerate, set the cross
5176  !*** product to a positive number.
5177  !***
5178 
5179  if (abs(cross_product) < tiny) then
5180  if (vec1_x .eq. zero .and. vec1_y .eq. zero) then
5181  cross_product = one
5182  lboundary = .false.
5183  else
5184  edgeid = n
5185  lboundary = .true.
5186  endif
5187  else
5188  lboundary = .false.
5189  endif
5190  endif ! if (lboundary)
5191 
5192  return ! pt in polygon
5193 
5194  endif ! if (inpoly)
5195 
5196  end do
5197 
5198  return ! pt outside polygon
5199 
5200 !----------------------------------------------------------------------
5201 

References scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, scrip_constants::pih, ptinpoly(), scrip_constants::tiny, and scrip_constants::zero.

Referenced by locate_segstart(), and ptincell().

◆ ptinpoly()

subroutine scrip_remap_conservative::ptinpoly ( real (scrip_r8), intent(in)  ptx,
real (scrip_r8), intent(in)  pty,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_x,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_y,
logical (scrip_logical), intent(in)  latlon,
logical (scrip_logical), intent(out)  inpoly,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 4815 of file scrip_remap_conservative.f.

4815 
4816 !----------------------------------------------------------------------
4817 !
4818 ! Check if point is in (convex) polygonal cell
4819 !
4820 !----------------------------------------------------------------------
4821 
4822  implicit none
4823 
4824 !----------------------------------------------------------------------
4825 !
4826 ! Input arguments
4827 !
4828 !----------------------------------------------------------------------
4829 
4830  real (SCRIP_r8), intent(in) ::
4831  & ptx, pty ! Point to check
4832 
4833  integer (SCRIP_i4), intent(in) ::
4834  & ncorners ! Number of polygon corners
4835 
4836  real (SCRIP_r8), dimension(ncorners), intent(in) ::
4837  & cell_corner_x, ! Coordinates of cell corners
4838  & cell_corner_y ! Could be x-y or lat-lon or ...
4839 
4840  logical (SCRIP_logical), intent(in) ::
4841  & latlon ! Are coordinates in latlon space?
4842 
4843 !----------------------------------------------------------------------
4844 !
4845 ! Output arguments
4846 !
4847 !----------------------------------------------------------------------
4848 
4849  logical (SCRIP_logical), intent(out) ::
4850  & inpoly ! Is point in the polygon?
4851 
4852  logical (SCRIP_logical), intent(out) ::
4853  & lboundary ! Is point on the boundary of the polygon?
4854 
4855  integer (SCRIP_i4), intent(out) ::
4856  & edgeid ! if point is on boundary, which local
4857  ! edge is it on? (0 otherwise)
4858 
4859 !----------------------------------------------------------------------
4860 !
4861 ! Local variables
4862 !
4863 !----------------------------------------------------------------------
4864 
4865  integer (SCRIP_i4) :: n, next_n
4866 
4867  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
4868  & cross_product, minlon, maxlon, ptx_loc, pty_loc
4869 
4870  real (SCRIP_r8), dimension(ncorners) ::
4871  & cell_corner_lat_loc, cell_corner_lon_loc
4872 
4873 
4874  !***********************************************************
4875  !*** We should just remove the latlon argument since that is
4876  !*** the only coordinate system we are using it for
4877  !***********************************************************
4878 
4879 
4880  !***
4881  !*** here we take the cross product of the vector making
4882  !*** up each cell side with the vector formed by the vertex
4883  !*** and search point. if all the cross products are
4884  !*** positive, the point is contained in the cell.
4885  !***
4886 
4887  inpoly = .false.
4888  lboundary = .false.
4889  edgeid = 0
4890 
4891  if (.not. latlon) then
4892 
4893  do n = 1, ncorners
4894  next_n = mod(n,ncorners) + 1
4895 
4896  x1 = cell_corner_x(n)
4897  y1 = cell_corner_y(n)
4898  x2 = cell_corner_x(next_n)
4899  y2 = cell_corner_y(next_n)
4900 
4901  vec1_x = x2 - x1
4902  vec1_y = y2 - y1
4903  vec2_x = ptx - x1
4904  vec2_y = pty - y1
4905 
4906  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
4907 
4908  !***
4909  !*** if the cross product for a side is zero, the point
4910  !*** lies exactly on the side or the side is degenerate
4911  !*** (zero length). if degenerate, set the cross
4912  !*** product to a positive number.
4913  !***
4914 
4915  if (abs(cross_product) < tiny) then
4916  if (vec1_x*vec1_x + vec1_y*vec1_y .le. tiny*tiny) then
4917  cross_product = one
4918  else
4919  lboundary = .true.
4920  edgeid = n
4921  endif
4922  else
4923 
4924  !***
4925  !*** if cross product is less than zero, this cell
4926  !*** doesn't work
4927  !***
4928  !*** Should we say "if (cp < zero .and. abs(cp) > tiny)" ?
4929 
4930  if (cross_product < zero) then
4931  inpoly = .false.
4932  lboundary = .false.
4933  return
4934  endif
4935  endif
4936 
4937  end do
4938 
4939  else
4940 
4941  !*** Checking in latlon space
4942  !*** If the grid cell coordinates spans more than pi radians
4943  !*** transform the coordinates so that they don't
4944 
4945  cell_corner_lat_loc = cell_corner_x
4946  cell_corner_lon_loc = cell_corner_y
4947 
4948  minlon = 9999.0
4949  maxlon = -9999.0
4950  do n = 1, ncorners
4951  if (cell_corner_lon_loc(n) < minlon) then
4952  minlon = cell_corner_lon_loc(n)
4953  endif
4954  if (cell_corner_lon_loc(n) > maxlon) then
4955  maxlon = cell_corner_lon_loc(n)
4956  endif
4957  enddo
4958 
4959  if (maxlon-minlon > pi) then
4960 
4961  do n = 1, ncorners
4962  if (cell_corner_lon_loc(n)-minlon > pi) then
4963  cell_corner_lon_loc(n) = cell_corner_lon_loc(n)-pi2
4964  endif
4965  enddo
4966 
4967  endif
4968 
4969  ptx_loc = ptx
4970  pty_loc = pty
4971  if (pty_loc - minlon > pi) then
4972  pty_loc = pty_loc - pi2
4973  else if (pty_loc - minlon < -pi) then
4974  pty_loc = pty_loc + pi2
4975  endif
4976 
4977 
4978  do n = 1, ncorners
4979  next_n = mod(n,ncorners) + 1
4980 
4981  x1 = cell_corner_lat_loc(n)
4982  y1 = cell_corner_lon_loc(n)
4983  x2 = cell_corner_lat_loc(next_n)
4984  y2 = cell_corner_lon_loc(next_n)
4985 
4986  vec1_x = x2 - x1
4987  vec1_y = y2 - y1
4988  vec2_x = ptx_loc - x1
4989  vec2_y = pty_loc - y1
4990 
4991  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
4992 
4993  !***
4994  !*** if the cross product for a side is zero, the point
4995  !*** lies exactly on the side or the side is degenerate
4996  !*** (zero length). if degenerate, set the cross
4997  !*** product to a positive number.
4998  !***
4999 
5000  if (abs(cross_product) < tiny) then
5001  if (vec1_x*vec1_x + vec1_y*vec1_y .le. tiny*tiny) then
5002  cross_product = one
5003  else
5004  lboundary = .true.
5005  edgeid = n
5006  endif
5007  else
5008 
5009  !***
5010  !*** if cross product is less than zero, this cell
5011  !*** doesn't work
5012  !***
5013  !*** Should we say "if (cp < zero .and. abs(cp) > tiny)" ?
5014 
5015  if (cross_product < zero) then
5016  inpoly = .false.
5017  lboundary = .false.
5018  return
5019  endif
5020  endif
5021 
5022  end do
5023 
5024  endif
5025  !***
5026  !*** if cross products all positive, we found the location
5027  !***
5028 
5029  inpoly = .true.
5030  return
5031 
5032 !----------------------------------------------------------------------
5033 

References scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, scrip_constants::tiny, and scrip_constants::zero.

Referenced by locate_segstart(), ptincell(), ptinpolarpoly(), and ptinpolygen().

◆ ptinpolygen()

subroutine scrip_remap_conservative::ptinpolygen ( real (scrip_r8), intent(in)  ptx,
real (scrip_r8), intent(in)  pty,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_x,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_y,
real (scrip_r8), intent(in)  cell_center_x,
real (scrip_r8), intent(in)  cell_center_y,
logical (scrip_logical), intent(in)  latlon,
logical (scrip_logical), intent(out)  inpoly,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 5211 of file scrip_remap_conservative.f.

5211 
5212 !----------------------------------------------------------------------
5213 !
5214 ! Check if point is in general (convex or mildly non-convex)
5215 ! polygonal cell by connecting each edge of the polygon to a
5216 ! a central point (average of vertices) and check containment in
5217 ! the resulting triangle
5218 !
5219 ! The cell can be non-convex as long as the 'center' is 'visible' to
5220 ! all the edges of the polygon, i.e., we can connect the 'center' to
5221 ! each edge of the polygon and form a triangle with positive area
5222 !
5223 !----------------------------------------------------------------------
5224 
5225  implicit none
5226 
5227 !----------------------------------------------------------------------
5228 !
5229 ! Input arguments
5230 !
5231 !----------------------------------------------------------------------
5232 
5233  real (SCRIP_r8), intent(in) ::
5234  & ptx, pty ! Point to check
5235 
5236  integer (SCRIP_i4), intent(in) ::
5237  & ncorners ! Number of polygon corners
5238 
5239  real (SCRIP_r8), dimension(ncorners), intent(in) ::
5240  & cell_corner_x, ! Coordinates of cell corners
5241  & cell_corner_y ! Could be x-y or lat-lon or ...
5242 
5243  real (SCRIP_r8), intent(in) ::
5244  & cell_center_x,
5245  & cell_center_y
5246 
5247  logical (SCRIP_logical), intent(in) ::
5248  & latlon ! Are coordinates in latlon space?
5249 
5250 !----------------------------------------------------------------------
5251 !
5252 ! Output arguments
5253 !
5254 !----------------------------------------------------------------------
5255 
5256  logical (SCRIP_logical), intent(out) ::
5257  & inpoly ! Is point in the polygon?
5258 
5259  logical (SCRIP_logical), intent(out) ::
5260  & lboundary ! Is point on the boundary of the polygon?
5261 
5262  integer (SCRIP_i4), intent(out) ::
5263  & edgeid ! if point is on boundary, which local
5264  ! edge is it on? (0 otherwise)
5265 
5266 !----------------------------------------------------------------------
5267 !
5268 ! Local variables
5269 !
5270 !----------------------------------------------------------------------
5271 
5272  integer (SCRIP_i4) :: n, next_n, ledgeid
5273 
5274  real (SCRIP_r8), dimension(3) ::
5275  & tri_corner_x, ! Coordinates of triangle
5276  & tri_corner_y
5277 
5278  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
5279  & cross_product
5280 
5281 
5282  !***
5283  !*** So, with each edge of the polygon form a triangle
5284  !*** by connecting a 'central' point to the endpoints of
5285  !*** the edge. Then check containment of the point in this tri
5286  !***
5287 
5288  inpoly = .false.
5289  lboundary = .false.
5290 
5291  do n = 1, ncorners
5292  next_n = mod(n,ncorners) + 1
5293 
5294  tri_corner_x(1) = cell_corner_x(n)
5295  tri_corner_y(1) = cell_corner_y(n)
5296  tri_corner_x(2) = cell_corner_x(next_n)
5297  tri_corner_y(2) = cell_corner_y(next_n)
5298  tri_corner_x(3) = cell_center_x
5299  tri_corner_y(3) = cell_center_y
5300 
5301  vec1_x = tri_corner_x(2) - tri_corner_x(1)
5302  vec1_y = tri_corner_y(2) - tri_corner_y(1)
5303 
5304  !*** Skip triangles arising from degenerate edges
5305 
5306  if (vec1_x*vec1_x+vec1_y*vec1_y .le. tiny*tiny) cycle
5307 
5308  call ptinpoly(ptx,pty,3,tri_corner_x,tri_corner_y,
5309  & latlon,inpoly,lboundary, ledgeid)
5310 
5311  if (inpoly) then
5312 
5313  if (lboundary) then
5314 
5315  !***
5316  !*** Check to see if the lboundary flag is being
5317  !*** triggered by the outer edge of the polygon or
5318  !*** by one of the artificial internal edges
5319  !***
5320 
5321  vec2_x = ptx - tri_corner_x(1)
5322  vec2_y = pty - tri_corner_y(1)
5323 
5324 
5325  if (latlon) then
5326 
5327  !***
5328  !*** check for 0,2pi crossings
5329  !***
5330 
5331  if (vec1_y > pi) vec1_y = vec1_y - pi2
5332  if (vec1_y < -pi) vec1_y = vec1_y + pi2
5333  if (vec2_y > pi) vec2_y = vec2_y - pi2
5334  if (vec2_y < -pi) vec2_y = vec2_y + pi2
5335 
5336  endif
5337 
5338  cross_product = vec1_y*vec2_x - vec2_y*vec1_x
5339 
5340  !***
5341  !*** if the cross product for a side is zero, the point
5342  !*** lies exactly on the side or the side is degenerate
5343  !*** (zero length). if degenerate, set the cross
5344  !*** product to a positive number.
5345  !***
5346 
5347  if (abs(cross_product) < tiny) then
5348  if (vec1_x*vec1_x+vec1_y*vec1_y .le. tiny*tiny) then
5349  cross_product = one
5350  lboundary = .false.
5351  else
5352  edgeid = n
5353  lboundary = .true.
5354  endif
5355  else
5356  lboundary = .false.
5357  endif
5358  endif ! if (lboundary)
5359 
5360  return ! pt in polygon
5361 
5362  endif ! if (inpoly)
5363 
5364  end do
5365 
5366  return ! pt outside polygon
5367 
5368 !----------------------------------------------------------------------
5369 

References scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, ptinpoly(), and scrip_constants::tiny.

◆ ptinpolygen2()

subroutine scrip_remap_conservative::ptinpolygen2 ( real (scrip_r8), intent(in)  ptx,
real (scrip_r8), intent(in)  pty,
integer (scrip_i4), intent(in)  ncorners,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_x,
real (scrip_r8), dimension(ncorners), intent(in)  cell_corner_y,
logical (scrip_logical), intent(in)  latlon,
logical (scrip_logical), intent(out)  inpoly,
logical (scrip_logical), intent(out)  lboundary,
integer (scrip_i4), intent(out)  edgeid 
)

Definition at line 5378 of file scrip_remap_conservative.f.

5378 
5379 !----------------------------------------------------------------------
5380 !
5381 ! Check if point is in general (convex or mildly non-convex)
5382 ! polygonal cell by connecting each edge of the polygon to a
5383 ! a central point (average of vertices) and check containment in
5384 ! the resulting triangle
5385 !
5386 ! The cell can be non-convex as long as the 'center' is 'visible' to
5387 ! all the edges of the polygon, i.e., we can connect the 'center' to
5388 ! each edge of the polygon and form a triangle with positive area
5389 !
5390 !----------------------------------------------------------------------
5391 
5392  implicit none
5393 
5394 !----------------------------------------------------------------------
5395 !
5396 ! Input arguments
5397 !
5398 !----------------------------------------------------------------------
5399 
5400  real (SCRIP_r8), intent(in) ::
5401  & ptx, pty ! Point to check
5402 
5403  integer (SCRIP_i4), intent(in) ::
5404  & ncorners ! Number of polygon corners
5405 
5406  real (SCRIP_r8), dimension(ncorners), intent(in) ::
5407  & cell_corner_x, ! Coordinates of cell corners
5408  & cell_corner_y ! Could be x-y or lat-lon or ...
5409 
5410  logical (SCRIP_logical), intent(in) ::
5411  & latlon ! Are coordinates in latlon space?
5412 
5413 !----------------------------------------------------------------------
5414 !
5415 ! Output arguments
5416 !
5417 !----------------------------------------------------------------------
5418 
5419  logical (SCRIP_logical), intent(out) ::
5420  & inpoly ! Is point in the polygon?
5421 
5422  logical (SCRIP_logical), intent(out) ::
5423  & lboundary ! Is point on the boundary of the polygon?
5424 
5425  integer (SCRIP_i4), intent(out) ::
5426  & edgeid ! if point is on boundary, which local
5427  ! edge is it on? (0 otherwise)
5428 
5429 !----------------------------------------------------------------------
5430 !
5431 ! Local variables
5432 !
5433 !----------------------------------------------------------------------
5434 
5435  integer (SCRIP_i4) :: c, n, next_n
5436 
5437  real (SCRIP_r8) :: x1, y1, x2, y2, vec1_x, vec1_y, vec2_x, vec2_y,
5438  & vec3_x, vec3_y, vec1_len, vec2_len, vec3_len,
5439  & cross_product, dot_product
5440 
5441 
5442  !***
5443  !*** So, with each edge of the polygon form a triangle
5444  !*** by connecting a 'central' point to the endpoints of
5445  !*** the edge. Then check containment of the point in this tri
5446  !***
5447 
5448  inpoly = .false.
5449  lboundary = .false.
5450 
5451  c = 0
5452  do n = 1, ncorners
5453  next_n = mod(n,ncorners) + 1
5454 
5455  x1 = cell_corner_x(n)
5456  y1 = cell_corner_y(n)
5457  x2 = cell_corner_x(next_n)
5458  y2 = cell_corner_y(next_n)
5459 
5460  if (((y1 > pty .and. y2 <= pty) .or.
5461  & (y2 > pty .and. y1 <= pty)) .and.
5462  & (ptx <= (x1 + (pty-y1)*(x2-x1)/(y2-y1)))) then
5463 
5464  c = 1 - c
5465 
5466  endif
5467  enddo
5468 
5469  if (c .eq. 1) inpoly = .true.
5470 
5471 
5472  !*** Check if the point is on the boundary of the polygon
5473 
5474  do n = 1, ncorners
5475 
5476  next_n = mod(n,ncorners) + 1
5477 
5478  x1 = cell_corner_x(n)
5479  y1 = cell_corner_y(n)
5480  x2 = cell_corner_x(next_n)
5481  y2 = cell_corner_y(next_n)
5482 
5483  vec1_x = x2 - x1
5484  vec1_y = y2 - y1
5485  vec1_len = sqrt(vec1_x*vec1_x + vec1_y*vec1_y)
5486  vec1_x = vec1_x/vec1_len
5487  vec1_y = vec1_y/vec1_len
5488 
5489  vec2_x = ptx - x1
5490  vec2_y = pty - y1
5491  vec2_len = sqrt(vec2_x*vec2_x + vec2_y*vec2_y)
5492 
5493  cross_product = vec1_x*vec2_y - vec2_x*vec1_y
5494  if (abs(cross_product) > tiny .and. vec2_len > tiny) then
5495  cross_product = cross_product/vec2_len
5496  endif
5497 
5498  if (abs(cross_product) < 1e5*tiny .and.
5499  & abs(cross_product) > 10*tiny) then
5500 
5501  !*** Sometimes when the point is too close to a vertex
5502  !*** then the cross product computation has errors due
5503  !*** to subtraction of two small numbers - So check w.r.t.
5504  !*** other vertex of the segment as well
5505 
5506  vec3_x = ptx - x2
5507  vec3_y = pty - y2
5508  vec3_len = sqrt(vec3_x*vec3_x + vec3_y*vec3_y)
5509 
5510  cross_product = -vec1_x*vec3_y + vec1_y*vec3_x
5511  if (abs(cross_product) > tiny .and. vec3_len > tiny) then
5512  !***
5513  !*** Normalize only if we won't be dividing two small
5514  !*** numbers
5515  cross_product = cross_product/vec3_len
5516  endif
5517  endif
5518 
5519  if (abs(cross_product) < 10*tiny) then
5520 
5521  if (vec1_x*vec1_x+vec1_y*vec1_y .le. tiny*tiny) then
5522  cross_product = one
5523  else
5524  dot_product = vec1_x*vec2_x + vec1_y*vec2_y
5525 
5526  if (dot_product >= 0 .and. dot_product <= vec1_len) then
5527  inpoly = .true.
5528  lboundary = .true.
5529  edgeid = n
5530  exit
5531  endif
5532  endif
5533  endif
5534 
5535  enddo
5536 
5537  return
5538 
5539 !----------------------------------------------------------------------
5540 

References scrip_constants::one, and scrip_constants::tiny.

Referenced by locate_segstart(), and ptincell().

◆ remap_conserv()

subroutine scrip_remap_conservative::remap_conserv ( logical(scrip_logical), intent(in)  l_master,
logical(scrip_logical), intent(in)  l_test 
)

Definition at line 246 of file scrip_remap_conservative.f.

246 
247 !-----------------------------------------------------------------------
248 !
249 ! this routine traces the perimeters of every grid cell on each
250 ! grid checking for intersections with the other grid and computing
251 ! line integrals for each subsegment.
252 !
253 !-----------------------------------------------------------------------
254 
255  logical(SCRIP_Logical), intent(in) :: l_master ! Am I the master
256  ! processor (do I/O)?
257  logical(SCRIP_Logical), intent(in) :: l_test ! Whether to
258  !include test output
259 
260 !-----------------------------------------------------------------------
261 !
262 ! local variables
263 !
264 !-----------------------------------------------------------------------
265 
266  integer (SCRIP_i4), parameter ::
267  & phi_or_theta = 2 ! integrate w.r.t. phi (1) or theta (2)
268 
269 
270  integer (SCRIP_i4) ::
271  & i, inext, !
272  & n, nwgt,
273  & grid1_add, ! Current linear address for grid1 cell
274  & grid2_add, ! Current linear address for grid2 cell
275  & grid_num, ! Index (1,2) of grid that we are
276  ! processing
277  & opp_grid_num, ! Index of opposite grid (2,1)
278  & maxrd_cell, ! cell with the max. relative difference
279  ! in area
280  & progint, ! Intervals at which progress is to be
281  ! printed
282  & icount ! for counting
283 
284  real (SCRIP_r8) ::
285  & norm_factor ! factor for normalizing wts
286 
287  real (SCRIP_r8), dimension(6) ::
288  & weights ! Weights array
289 
290  real (SCRIP_r8) ::
291  & beglat, beglon,
292  & endlat, endlon,
293  & ave_reldiff, ! Average rel. diff. in areas
294  & max_reldiff, ! Maximum rel. diff in areas
295  & maxrd_area, ! Computed area for cell with max rel
296  ! diff
297  & maxrd_true ! True area for cell with max rel diff
298 
299  real (SCRIP_r8), dimension(:), allocatable ::
300  & reldiff, ! Relative difference in computed
301  ! and true area
302  & ref_area ! Area of cell as computed by direct
303  ! integration around its boundaries
304 
305 ! call OMP_SET_DYNAMIC(.FALSE.)
306 
307 !-----------------------------------------------------------------------
308 !
309 ! integrate around each cell on grid1
310 !
311 !-----------------------------------------------------------------------
312 
313  is_master=l_master ! set module variable using subroutine input
314  ! argument variable.
315  ! Use the former subsequently.
316 
317  if(is_master)print *,'grid1 sweep'
318 
319 !NRL Progress is slow when the other grid (grid 2) is large, so we use
320 !NRL that. Really, it would be a better to do this with a timer...
321  if (grid2_size > 500000) then
322  progint = 1000
323  elseif (grid2_size > 250000) then
324  progint = 2000
325  elseif (grid2_size > 100000) then
326  progint = 5000
327  else
328  progint = 10000
329  endif
330 
331  grid_num = 1
332  opp_grid_num = 2
333 
334  call timer_start(1)
335 
336 C$OMP PARALLEL DEFAULT(SHARED) PRIVATE(grid1_add) NUM_THREADS(nthreads)
337 
338 C$OMP DO SCHEDULE(DYNAMIC)
339 
340  do grid1_add = 1,grid1_size
341 
342  if (mod(grid1_add,progint) .eq. 0 .and. is_master) then
343  print *, grid1_add,' of ',grid1_size,' cells processed ...'
344  endif
345 
346  call cell_integrate(grid1_add, grid_num, phi_or_theta)
347 
348  end do ! do grid1_add=...
349 
350 C$OMP END DO
351 
352 C$OMP END PARALLEL
353 
354 !-----------------------------------------------------------------------
355 !
356 ! integrate around each cell on grid2
357 !
358 !-----------------------------------------------------------------------
359 
360  if(is_master)print *,'grid2 sweep '
361 
362 !NRL Progress is slow when the other grid (grid 1) is large, so we use
363 !NRL that.
364  if (grid1_size > 500000) then
365  progint = 1000
366  elseif (grid1_size > 250000) then
367  progint = 2000
368  elseif (grid1_size > 100000) then
369  progint = 5000
370  else
371  progint = 10000
372  endif
373 
374  grid_num = 2
375  opp_grid_num = 1
376 
377  call timer_start(2)
378 
379 C$OMP PARALLEL DEFAULT(SHARED) PRIVATE(grid2_add) NUM_THREADS(nthreads)
380 
381 C$OMP DO SCHEDULE(DYNAMIC)
382 
383  do grid2_add = 1,grid2_size
384 
385  if (mod(grid2_add,progint) .eq. 0 .and. is_master) then
386  print *, grid2_add,' of ',grid2_size,' cells processed ...'
387  endif
388 
389  call cell_integrate(grid2_add, grid_num, phi_or_theta)
390 
391  end do ! do grid2_add=...
392 
393 C$OMP END DO
394 
395 C$OMP END PARALLEL
396 
397  call timer_stop(2)
398 
399 !-----------------------------------------------------------------------
400 !
401 ! correct for situations where N/S pole not explicitly included in
402 ! grid (i.e. as a grid corner point). if pole is missing from only
403 ! one grid, need to correct only the area and centroid of that
404 ! grid. if missing from both, do complete weight calculation.
405 ! This is necessary only when integrating w.r.t. phi (longitude)
406 !
407 !-----------------------------------------------------------------------
408 
409  if (phi_or_theta .eq. 1) then
410 
411  !*** North Pole
412  weights(1) = pi2
413  weights(2) = pi*pi
414  weights(3) = zero
415  weights(4) = pi2
416  weights(5) = pi*pi
417  weights(6) = zero
418 
419  if (grid1_npole_cell /=0) then
420  grid1_area(grid1_npole_cell) = grid1_area(grid1_npole_cell)
421  & + weights(1)
422  grid1_centroid_lat(grid1_npole_cell) =
423  & grid1_centroid_lat(grid1_npole_cell) + weights(2)
424  grid1_centroid_lon(grid1_npole_cell) =
425  & grid1_centroid_lon(grid1_npole_cell) + weights(3)
426  endif
427 
428  if (grid2_npole_cell /=0) then
429  grid2_area(grid2_npole_cell) = grid2_area(grid2_npole_cell)
430  & + weights(num_wts+1)
431  grid2_centroid_lat(grid2_npole_cell) =
432  & grid2_centroid_lat(grid2_npole_cell) +
433  & weights(num_wts+2)
434  grid2_centroid_lon(grid2_npole_cell) =
435  & grid2_centroid_lon(grid2_npole_cell) +
436  & weights(num_wts+3)
437  endif
438 
439  if (grid1_npole_cell /= 0 .and. grid2_npole_cell /=0) then
440  call store_link_cnsrv(grid1_npole_cell,
441  & grid2_npole_cell, weights)
442 
443  grid1_frac(grid1_npole_cell) = grid1_frac(grid1_npole_cell)
444  & + weights(1)
445  grid2_frac(grid2_npole_cell) = grid2_frac(grid2_npole_cell)
446  & + weights(num_wts+1)
447  endif
448 
449 
450  !*** South Pole
451  weights(1) = pi2
452  weights(2) = -pi*pi
453  weights(3) = zero
454  weights(4) = pi2
455  weights(5) = -pi*pi
456  weights(6) = zero
457 
458  if (grid1_spole_cell /=0) then
459  grid1_area(grid1_spole_cell) = grid1_area(grid1_spole_cell)
460  & + weights(1)
461  grid1_centroid_lat(grid1_spole_cell) =
462  & grid1_centroid_lat(grid1_spole_cell) + weights(2)
463  grid1_centroid_lon(grid1_spole_cell) =
464  & grid1_centroid_lon(grid1_spole_cell) + weights(3)
465  endif
466 
467  if (grid2_spole_cell /=0) then
468  grid2_area(grid2_spole_cell) = grid2_area(grid2_spole_cell)
469  & + weights(num_wts+1)
470  grid2_centroid_lat(grid2_spole_cell) =
471  & grid2_centroid_lat(grid2_spole_cell) +
472  & weights(num_wts+2)
473  grid2_centroid_lon(grid2_spole_cell) =
474  & grid2_centroid_lon(grid2_spole_cell) +
475  & weights(num_wts+3)
476  endif
477 
478  if (grid1_spole_cell /= 0 .and. grid2_spole_cell /=0) then
479  call store_link_cnsrv(grid1_spole_cell,
480  & grid2_spole_cell, weights)
481 
482  grid1_frac(grid1_spole_cell) = grid1_frac(grid1_spole_cell)
483  & + weights(1)
484  grid2_frac(grid2_spole_cell) = grid2_frac(grid2_spole_cell)
485  & + weights(num_wts+1)
486  endif
487  endif
488 
489 
490 
491  if(is_master)print *, 'Grid sweeps completed'
492 
493 
494 !-----------------------------------------------------------------------
495 !
496 ! finish centroid computation
497 !
498 !-----------------------------------------------------------------------
499 
500  call timer_start(3)
501 
502 C$OMP PARALLEL
503 C$OMP WORKSHARE
504  where (grid1_area /= zero)
505  grid1_centroid_lat = grid1_centroid_lat/grid1_area
506  grid1_centroid_lon = grid1_centroid_lon/grid1_area
507  end where
508 C$OMP END WORKSHARE
509 
510 C$OMP WORKSHARE
511  where (grid2_area /= zero)
512  grid2_centroid_lat = grid2_centroid_lat/grid2_area
513  grid2_centroid_lon = grid2_centroid_lon/grid2_area
514  end where
515 C$OMP END WORKSHARE
516 C$OMP END PARALLEL
517 
518 
519 !-----------------------------------------------------------------------
520 !
521 ! include centroids in weights and normalize using destination
522 ! area if requested
523 !
524 !-----------------------------------------------------------------------
525 
526 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
527 C$OMP& PRIVATE(n,grid1_add,grid2_add,nwgt,weights,norm_factor)
528 
529 C$OMP DO SCHEDULE(DYNAMIC)
530 
531  do n=1,num_links_map1
532  grid1_add = grid1_add_map1(n)
533  grid2_add = grid2_add_map1(n)
534  do nwgt=1,num_wts
535  weights( nwgt) = wts_map1(nwgt,n)
536  if (num_maps > 1) then
537  weights(num_wts+nwgt) = wts_map2(nwgt,n)
538  endif
539  end do
540 
541  select case(norm_opt)
542  case (norm_opt_dstarea)
543  if (grid2_area(grid2_add) /= zero) then
544  if (luse_grid2_area) then
545  norm_factor = one/grid2_area_in(grid2_add)
546  else
547  norm_factor = one/grid2_area(grid2_add)
548  endif
549  else
550  norm_factor = zero
551  endif
552  case (norm_opt_frcarea)
553  if (grid2_frac(grid2_add) /= zero) then
554  if (luse_grid2_area) then
555  norm_factor = grid2_area(grid2_add)/
556  & (grid2_frac(grid2_add)*
557  & grid2_area_in(grid2_add))
558  else
559  norm_factor = one/grid2_frac(grid2_add)
560  endif
561  else
562  norm_factor = zero
563  endif
564  case (norm_opt_none)
565  norm_factor = one
566  end select
567 
568  wts_map1(1,n) = weights(1)*norm_factor
569  wts_map1(2,n) = (weights(2) - weights(1)*
570  & grid1_centroid_lat(grid1_add))*
571  & norm_factor
572  wts_map1(3,n) = (weights(3) - weights(1)*
573  & grid1_centroid_lon(grid1_add))*
574  & norm_factor
575 
576  if (num_maps > 1) then
577  select case(norm_opt)
578  case (norm_opt_dstarea)
579  if (grid1_area(grid1_add) /= zero) then
580  if (luse_grid1_area) then
581  norm_factor = one/grid1_area_in(grid1_add)
582  else
583  norm_factor = one/grid1_area(grid1_add)
584  endif
585  else
586  norm_factor = zero
587  endif
588  case (norm_opt_frcarea)
589  if (grid1_frac(grid1_add) /= zero) then
590  if (luse_grid1_area) then
591  norm_factor = grid1_area(grid1_add)/
592  & (grid1_frac(grid1_add)*
593  & grid1_area_in(grid1_add))
594  else
595  norm_factor = one/grid1_frac(grid1_add)
596  endif
597  else
598  norm_factor = zero
599  endif
600  case (norm_opt_none)
601  norm_factor = one
602  end select
603 
604  wts_map2(1,n) = weights(num_wts+1)*norm_factor
605  wts_map2(2,n) = (weights(num_wts+2) - weights(num_wts+1)*
606  & grid2_centroid_lat(grid2_add))*
607  & norm_factor
608  wts_map2(3,n) = (weights(num_wts+3) - weights(num_wts+1)*
609  & grid2_centroid_lon(grid2_add))*
610  & norm_factor
611  endif
612 
613  end do
614 
615 C$OMP END DO
616 
617 C$OMP END PARALLEL
618 
619  if(is_master)print *, 'Total number of links = ',num_links_map1
620 
621 C$OMP PARALLEL
622 C$OMP WORKSHARE
623  where (grid1_area /= zero) grid1_frac = grid1_frac/grid1_area
624 C$OMP END WORKSHARE
625 C$OMP WORKSHARE
626  where (grid2_area /= zero) grid2_frac = grid2_frac/grid2_area
627 C$OMP END WORKSHARE
628 C$OMP END PARALLEL
629 
630  call timer_stop(3)
631 
632 !-----------------------------------------------------------------------
633 !
634 ! perform some error checking on final weights
635 !
636 !-----------------------------------------------------------------------
637 
638  allocate(ref_area(grid1_size))
639  allocate(reldiff(grid1_size))
640 
641 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
642 C$OMP& PRIVATE(n, i, inext, beglat, beglon, endlat, endlon, weights)
643 C$OMP DO SCHEDULE(DYNAMIC)
644 
645  do n=1,grid1_size
646  if (grid1_area(n) < -.01 .and. is_master) then
647  print *,'Grid 1 area error: ',n,grid1_area(n)
648  endif
649  if ((grid1_centroid_lat(n) < -pih-.01 .or.
650  & grid1_centroid_lat(n) > pih+.01) .and. is_master) then
651  print *,'Grid 1 centroid lat error: ',n,grid1_centroid_lat(n)
652  endif
653 
654  ref_area(n) = 0.0
655  do i = 1, grid1_corners
656  inext = 1 + mod(i,grid1_corners)
657 
658  beglat = grid1_corner_lat(i,n)
659  beglon = grid1_corner_lon(i,n)
660  endlat = grid1_corner_lat(inext,n)
661  endlon = grid1_corner_lon(inext,n)
662 
663  if ((phi_or_theta .eq. 1 .and. beglon .eq. endlon) .or.
664  & (phi_or_theta .eq. 2 .and. beglat .eq. endlat)) cycle
665 
666  call line_integral(phi_or_theta, weights, num_wts, beglon,
667  & endlon, beglat, endlat, grid1_center_lat(n),
668  & grid1_center_lon(n), grid1_center_lat(n),
669  & grid1_center_lon(n))
670 
671  ref_area(n) = ref_area(n) + weights(1)
672  enddo
673  enddo
674 C$OMP END DO
675 C$OMP END PARALLEL
676 
677 
678 ! Correct for polar cells
679 
680  if (phi_or_theta .eq. 1) then
681 
682  !*** North Pole
683  weights(1) = pi2
684 
685  if (grid1_npole_cell /=0) then
686  ref_area(grid1_npole_cell) = ref_area(grid1_npole_cell)
687  & + weights(1)
688  endif
689 
690  !*** South Pole
691  weights(1) = pi2
692 
693  if (grid1_spole_cell /=0) then
694  ref_area(grid1_spole_cell) = ref_area(grid1_spole_cell)
695  & + weights(1)
696  endif
697 
698  endif
699 
700 
701  ave_reldiff = 0.0
702  max_reldiff = -1.0
703 
704  do n = 1, grid1_size
705  if(ref_area(n).gt.0.0)then ! added May 21 2013
706  reldiff(n) = abs(ref_area(n)-grid1_area(n))/abs(ref_area(n))
707  endif
708  ave_reldiff = ave_reldiff + reldiff(n)
709  if (reldiff(n) > max_reldiff) then
710  max_reldiff = reldiff(n)
711  maxrd_cell = n
712  maxrd_area = grid1_area(n)
713  maxrd_true = ref_area(n)
714  endif
715  end do
716 
717  ave_reldiff = ave_reldiff/grid1_size
718 
719  if(is_master.and.l_test)then
720  print *
721  print *
722  print *,'Grid 1: Ave. rel. diff. in areas: ',
723  & ave_reldiff
724  print *,' rel. diff. = abs(area-refarea)/refarea'
725  print *
726  print *,'Grid 1: Max. rel. diff. in areas: ',
727  & max_reldiff
728  print *, 'Max rel. diff. is in cell ',maxrd_cell
729  print *, 'Computed Area: ', maxrd_area
730  print *, 'Reference Area: ',maxrd_true
731  print *
732  endif
733 
734  deallocate(ref_area, reldiff)
735 
736 
737 
738  allocate(ref_area(grid2_size))
739  allocate(reldiff(grid2_size))
740 
741 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
742 C$OMP& PRIVATE(n, i, inext, beglat, beglon, endlat, endlon, weights)
743 C$OMP DO SCHEDULE(DYNAMIC)
744 
745  do n=1,grid2_size
746  if (grid2_area(n) < -.01 .and. is_master) then
747  print *,'Grid 2 area error: ',n,grid2_area(n)
748  endif
749  if ((grid2_centroid_lat(n) < -pih-.01 .or.
750  & grid2_centroid_lat(n) > pih+.01) .and. is_master) then
751  print *,'Grid 2 centroid lat error: ',n,grid2_centroid_lat(n)
752  endif
753 
754  ref_area(n) = 0.0
755  do i = 1, grid2_corners
756  inext = 1 + mod(i,grid2_corners)
757 
758  beglat = grid2_corner_lat(i,n)
759  beglon = grid2_corner_lon(i,n)
760  endlat = grid2_corner_lat(inext,n)
761  endlon = grid2_corner_lon(inext,n)
762 
763  if ((phi_or_theta .eq. 1 .and. beglon .eq. endlon) .or.
764  & (phi_or_theta .eq. 2 .and. beglat .eq. endlat)) cycle
765 
766  call line_integral(phi_or_theta, weights, num_wts, beglon,
767  & endlon, beglat, endlat, grid2_center_lat(n),
768  & grid2_center_lon(n), grid2_center_lat(n),
769  & grid2_center_lon(n))
770 
771  ref_area(n) = ref_area(n) + weights(1)
772  enddo
773  enddo
774 C$OMP END DO
775 C$OMP END PARALLEL
776 
777 
778 ! Correct for polar cells
779 
780  if (phi_or_theta .eq. 1) then
781 
782  !*** North Pole
783  weights(1) = pi2
784 
785  if (grid2_npole_cell /=0) then
786  ref_area(grid2_npole_cell) = ref_area(grid2_npole_cell)
787  & + weights(1)
788  endif
789 
790  !*** South Pole
791  weights(1) = pi2
792 
793  if (grid2_spole_cell /=0) then
794  ref_area(grid2_spole_cell) = ref_area(grid2_spole_cell)
795  & + weights(1)
796  endif
797 
798  endif
799 
800 
801  ave_reldiff = 0.0
802  max_reldiff = -1.0
803 
804  do n = 1, grid2_size
805  reldiff(n) = abs(ref_area(n)-grid2_area(n))/abs(ref_area(n))
806  ave_reldiff = ave_reldiff + reldiff(n)
807  if (reldiff(n) > max_reldiff) then
808  max_reldiff = reldiff(n)
809  maxrd_cell = n
810  maxrd_area = grid2_area(n)
811  maxrd_true = ref_area(n)
812  endif
813  end do
814 
815  ave_reldiff = ave_reldiff/grid2_size
816 
817  if(is_master.and.l_test)then
818  print *
819  print *,'Grid 2: Ave. rel. diff. in areas: ',
820  & ave_reldiff
821  print *,' rel. diff. = abs(area-refarea)/refarea'
822  print *
823  print *,'Grid 2: Max. rel. diff. in areas: ',
824  & max_reldiff
825  print *, 'Max rel. diff. is in cell ',maxrd_cell
826  print *, 'Computed Area: ', maxrd_area
827  print *, 'Reference Area: ',maxrd_true
828  print *
829  endif
830 
831  deallocate(ref_area,reldiff)
832 
833  if(is_master.and.l_test)then
834  print *, 'Computed area = Area of cell computed by adding areas'
835  print *, ' of intersection with other cells'
836  print *, 'Reference area = Area of cell by direct integration'
837  print *
838  endif
839 
840  !***
841  !*** In the following code, gridN_centroid_lat is being used to
842  !*** store running tallies of the cell areas - so it is a
843  !*** misnomer used to avoid allocation of a new variable
844  !***
845 
846  grid1_centroid_lat = zero
847  grid2_centroid_lat = zero
848  icount=0
849 
850 C$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nthreads)
851 C$OMP& PRIVATE(n,grid1_add,grid2_add,nwgt,weights)
852 C$OMP DO SCHEDULE(DYNAMIC)
853 
854  do n=1,num_links_map1
855  grid1_add = grid1_add_map1(n)
856  grid2_add = grid2_add_map1(n)
857 
858  do nwgt=1,num_wts
859  weights( nwgt) = wts_map1(nwgt,n)
860  if (num_maps > 1) then
861  weights(num_wts+nwgt) = wts_map2(nwgt,n)
862  endif
863  end do
864 
865 ! count warnings about weights that will be excluded
866  if (grid2_frac(grid2_add).gt.frac_lowest .and.
867  & grid2_frac(grid2_add).lt.frac_highest .and. is_master) then
868  if ( (wts_map1(1,n) < wt_lowest) )then
869  icount=icount+1
870 ! print statements that were here have been moved to another routine...
871  endif
872  if (norm_opt /= norm_opt_none .and. wts_map1(1,n) >
873  & wt_highest)then
874  icount=icount+1
875 ! print statements that were here have been moved to another routine...
876  endif
877  endif
878 C$OMP CRITICAL
879  grid2_centroid_lat(grid2_add) =
880  & grid2_centroid_lat(grid2_add) + wts_map1(1,n)
881 C$OMP END CRITICAL
882 
883  if (num_maps > 1) then
884  if (wts_map2(1,n) < -.01 .and. is_master) then
885  print *,'Map 2 weight < 0 ',grid1_add,grid2_add,
886  & wts_map2(1,n)
887  endif
888  if (norm_opt /= norm_opt_none .and. wts_map2(1,n) > 1.01
889  & .and. is_master) then
890  print *,'Map 2 weight > 1 ',grid1_add,grid2_add,
891  & wts_map2(1,n)
892  endif
893 C$OMP CRITICAL
894  grid1_centroid_lat(grid1_add) =
895  & grid1_centroid_lat(grid1_add) + wts_map2(1,n)
896 C$OMP END CRITICAL
897  endif
898  end do
899 
900 C$OMP END DO
901 C$OMP END PARALLEL
902 
903  if(icount.gt.0.and.is_master)then
904  print *,'We had problems in ',icount,' points.'
905  endif
906 ! stop condition was here...has been moved to another routine...
907 
908  !***
909  !*** If grid1 has masks, links between some cells of grid1 and
910  !*** grid2 do not exist even though they overlap. In such a case,
911  !*** the following code will generate errors even though nothing
912  !*** is wrong (grid1_centroid_lat or grid2_centroid_lat are never
913  !*** updated in the above loop)
914  !***
915 
916  do n=1,grid2_size
917  select case(norm_opt)
918  case (norm_opt_dstarea)
919  norm_factor = grid2_frac(n)
920  case (norm_opt_frcarea)
921  norm_factor = one
922  case (norm_opt_none)
923  if (luse_grid2_area) then
924  norm_factor = grid2_area_in(n)
925  else
926  norm_factor = grid2_area(n)
927  endif
928  end select
929 ! if (abs(grid2_centroid_lat(n)-norm_factor) > .01
930 ! & .and. is_master) then
931 ! print *,'Warning: sum of wts for map1 ',n,
932 ! & grid2_centroid_lat(n),norm_factor
933 ! endif
934 ! write(501,*)n,grid2_centroid_lat(n)
935  end do
936 
937 
938  if (num_maps > 1) then
939  do n=1,grid1_size
940  select case(norm_opt)
941  case (norm_opt_dstarea)
942  norm_factor = grid1_frac(n)
943  case (norm_opt_frcarea)
944  norm_factor = one
945  case (norm_opt_none)
946  if (luse_grid1_area) then
947  norm_factor = grid1_area_in(n)
948  else
949  norm_factor = grid1_area(n)
950  endif
951  end select
952  if (abs(grid1_centroid_lat(n)-norm_factor) > .01
953  & .and. is_master) then
954  print *,'Error: sum of wts for map2 ',n,
955  & grid1_centroid_lat(n),norm_factor
956  endif
957  end do
958  endif
959 !-----------------------------------------------------------------------
960 
961  call timer_stop(4)
962 
963  if(is_master)print *, 'Finished Conservative Remapping'
964 
965  if(l_test)then
966  call timer_print(1)
967  call timer_print(2)
968  call timer_print(3)
969  call timer_print(4)
970  endif
971 

References cell_integrate(), scrip_remap_vars::frac_highest, scrip_remap_vars::frac_lowest, scrip_remap_vars::grid1_add_map1, scrip_grids::grid1_area, scrip_grids::grid1_area_in, scrip_grids::grid1_center_lat, scrip_grids::grid1_center_lon, scrip_grids::grid1_centroid_lat, scrip_grids::grid1_centroid_lon, scrip_grids::grid1_corner_lat, scrip_grids::grid1_corner_lon, scrip_grids::grid1_corners, scrip_grids::grid1_frac, scrip_grids::grid1_npole_cell, scrip_grids::grid1_size, scrip_grids::grid1_spole_cell, scrip_remap_vars::grid2_add_map1, scrip_grids::grid2_area, scrip_grids::grid2_area_in, scrip_grids::grid2_center_lat, scrip_grids::grid2_center_lon, scrip_grids::grid2_centroid_lat, scrip_grids::grid2_centroid_lon, scrip_grids::grid2_corner_lat, scrip_grids::grid2_corner_lon, scrip_grids::grid2_corners, scrip_grids::grid2_frac, scrip_grids::grid2_npole_cell, scrip_grids::grid2_size, scrip_grids::grid2_spole_cell, line_integral(), scrip_grids::luse_grid1_area, scrip_grids::luse_grid2_area, scrip_remap_vars::norm_opt, scrip_remap_vars::norm_opt_dstarea, scrip_remap_vars::norm_opt_frcarea, scrip_remap_vars::norm_opt_none, scrip_remap_vars::num_links_map1, scrip_remap_vars::num_maps, scrip_remap_vars::num_wts, scrip_constants::one, scrip_constants::pi, scrip_constants::pi2, scrip_constants::pih, store_link_cnsrv(), scrip_timers::timer_print(), scrip_timers::timer_start(), scrip_timers::timer_stop(), scrip_remap_vars::wt_highest, scrip_remap_vars::wt_lowest, scrip_remap_vars::wts_map1, scrip_remap_vars::wts_map2, and scrip_constants::zero.

Referenced by scrip_interface::scrip().

◆ store_link_cnsrv()

subroutine scrip_remap_conservative::store_link_cnsrv ( integer (scrip_i4), intent(in)  add1,
integer (scrip_i4), intent(in)  add2,
real (scrip_r8), dimension(:), intent(in)  weights 
)

Definition at line 3951 of file scrip_remap_conservative.f.

3951 
3952 !-----------------------------------------------------------------------
3953 !
3954 ! this routine stores the address and weight for this link in
3955 ! the appropriate address and weight arrays and resizes those
3956 ! arrays if necessary.
3957 !
3958 !-----------------------------------------------------------------------
3959 
3960 !-----------------------------------------------------------------------
3961 !
3962 ! input variables
3963 !
3964 !-----------------------------------------------------------------------
3965 
3966  integer (SCRIP_i4), intent(in) ::
3967  & add1, ! address on grid1
3968  & add2 ! address on grid2
3969 
3970  real (SCRIP_r8), dimension(:), intent(in) ::
3971  & weights ! array of remapping weights for this link
3972 
3973 !-----------------------------------------------------------------------
3974 !
3975 ! local variables
3976 !
3977 !-----------------------------------------------------------------------
3978 
3979  integer (SCRIP_i4) :: nlink, min_link, max_link ! link index
3980 
3981  logical (SCRIP_logical) :: found
3982 
3983 
3984 !-----------------------------------------------------------------------
3985 !
3986 ! if all weights are zero, do not bother storing the link
3987 !
3988 !-----------------------------------------------------------------------
3989 
3990  if (all(weights == zero)) return
3991 
3992 !-----------------------------------------------------------------------
3993 !
3994 ! restrict the range of links to search for existing links
3995 !
3996 !-----------------------------------------------------------------------
3997 
3998 C$OMP CRITICAL(block5)
3999 ! first_call should be within critical block or else multiple
4000 ! threads will see it as true the first time around
4001 
4002  if (first_call_store_link_cnsrv) then
4003  allocate(link_add1(2,grid1_size), link_add2(2,grid2_size))
4004  link_add1 = 0
4005  link_add2 = 0
4006  first_call_store_link_cnsrv = .false.
4007  min_link = 1
4008  max_link = 0
4009  else
4010  min_link = min(link_add1(1,add1),link_add2(1,add2))
4011  max_link = max(link_add1(2,add1),link_add2(2,add2))
4012  if (min_link == 0) then
4013  min_link = 1
4014  max_link = 0
4015  endif
4016  endif
4017 C$OMP END CRITICAL(block5)
4018 
4019 !-----------------------------------------------------------------------
4020 !
4021 ! if the link already exists, add the weight to the current weight
4022 ! arrays
4023 !
4024 !-----------------------------------------------------------------------
4025 
4026  found = .false.
4027 
4028  do nlink=min_link,max_link
4029  if (add1 == grid1_add_map1(nlink)) then
4030  if (add2 == grid2_add_map1(nlink)) then
4031 
4032 C$OMP CRITICAL(block3a)
4033  wts_map1(:,nlink) = wts_map1(:,nlink) + weights(1:num_wts)
4034  if (num_maps == 2) then
4035  wts_map2(:,nlink) = wts_map2(:,nlink) +
4036  & weights(num_wts+1:2*num_wts)
4037  endif
4038 C$OMP END CRITICAL(block3a)
4039  found = .true.
4040  exit
4041 
4042  endif
4043  endif
4044  end do
4045 
4046 
4047  if (found) return
4048 
4049 !-----------------------------------------------------------------------
4050 !
4051 ! if the link does not yet exist, increment number of links and
4052 ! check to see if remap arrays need to be increased to accomodate
4053 ! the new link. then store the link.
4054 !
4055 !-----------------------------------------------------------------------
4056 
4057 C$OMP CRITICAL(block6)
4058 
4059  num_links_map1 = num_links_map1 + 1
4060  if (num_links_map1 > max_links_map1)
4061  & call resize_remap_vars(1,resize_increment)
4062 
4063  grid1_add_map1(num_links_map1) = add1
4064  grid2_add_map1(num_links_map1) = add2
4065  wts_map1(:,num_links_map1) = weights(1:num_wts)
4066 
4067  if (num_maps > 1) then
4068  num_links_map2 = num_links_map2 + 1
4069  if (num_links_map2 > max_links_map2)
4070  & call resize_remap_vars(2,resize_increment)
4071 
4072  grid1_add_map2(num_links_map2) = add1
4073  grid2_add_map2(num_links_map2) = add2
4074  wts_map2(:,num_links_map2) = weights(num_wts+1:2*num_wts)
4075  endif
4076 
4077  if (link_add1(1,add1) == 0) link_add1(1,add1) = num_links_map1
4078  if (link_add2(1,add2) == 0) link_add2(1,add2) = num_links_map1
4079  link_add1(2,add1) = num_links_map1
4080  link_add2(2,add2) = num_links_map1
4081 
4082 C$OMP END CRITICAL(block6)
4083 
4084 !-----------------------------------------------------------------------
4085 

References first_call_store_link_cnsrv, scrip_remap_vars::grid1_add_map1, scrip_remap_vars::grid1_add_map2, scrip_grids::grid1_size, scrip_remap_vars::grid2_add_map1, scrip_remap_vars::grid2_add_map2, scrip_grids::grid2_size, link_add1, link_add2, scrip_remap_vars::max_links_map1, scrip_remap_vars::max_links_map2, scrip_remap_vars::num_links_map1, scrip_remap_vars::num_links_map2, scrip_remap_vars::num_maps, scrip_remap_vars::num_wts, scrip_remap_vars::resize_increment, scrip_remap_vars::resize_remap_vars(), scrip_remap_vars::wts_map1, scrip_remap_vars::wts_map2, and scrip_constants::zero.

Referenced by cell_integrate(), and remap_conserv().

Variable Documentation

◆ avoid_pole_count

integer (scrip_i4), save scrip_remap_conservative::avoid_pole_count = 0

Definition at line 80 of file scrip_remap_conservative.f.

80  integer (SCRIP_i4), save ::
81  & avoid_pole_count = 0 ! count attempts to avoid pole

Referenced by pole_intersection(), and scrip_interface::scrip_clear().

◆ avoid_pole_offset

real (scrip_r8), save scrip_remap_conservative::avoid_pole_offset = tiny

Definition at line 83 of file scrip_remap_conservative.f.

83  real (SCRIP_r8), save ::
84  & avoid_pole_offset = tiny ! endpoint offset to avoid pole

Referenced by pole_intersection(), and scrip_interface::scrip_clear().

◆ first_call_find_adj_cell

logical (scrip_logical), save scrip_remap_conservative::first_call_find_adj_cell =.true.

Definition at line 174 of file scrip_remap_conservative.f.

174  logical (SCRIP_logical), save ::
175  & first_call_find_adj_cell=.true.

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ first_call_get_srch_cells

logical (scrip_logical), save scrip_remap_conservative::first_call_get_srch_cells =.true.

Definition at line 171 of file scrip_remap_conservative.f.

171  logical (SCRIP_logical), save ::
172  & first_call_get_srch_cells=.true.

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ first_call_locate_point

logical (scrip_logical), save scrip_remap_conservative::first_call_locate_point = .true.

Definition at line 121 of file scrip_remap_conservative.f.

121  logical (SCRIP_logical), save ::
122  & first_call_locate_point= .true.

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ first_call_locate_segstart

logical (scrip_logical), save scrip_remap_conservative::first_call_locate_segstart = .true.

Definition at line 93 of file scrip_remap_conservative.f.

93  logical (SCRIP_logical), save ::
94  & first_call_locate_segstart= .true.

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ first_call_store_link_cnsrv

logical (scrip_logical), save scrip_remap_conservative::first_call_store_link_cnsrv = .true.

Definition at line 90 of file scrip_remap_conservative.f.

90  logical (SCRIP_logical), save ::
91  & first_call_store_link_cnsrv = .true.

Referenced by scrip_interface::scrip_clear(), and store_link_cnsrv().

◆ last_cell_add_get_srch_cells

integer (scrip_i4), save scrip_remap_conservative::last_cell_add_get_srch_cells

Definition at line 166 of file scrip_remap_conservative.f.

166  integer (SCRIP_i4), save ::
167  & last_cell_add_get_srch_cells,
168  & last_cell_grid_num_get_srch_cells,
169  & last_srch_grid_num_get_srch_cells

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ last_cell_find_adj_cell

integer (scrip_i4), save scrip_remap_conservative::last_cell_find_adj_cell

Definition at line 180 of file scrip_remap_conservative.f.

180  integer (SCRIP_i4), save ::
181  & last_cell_find_adj_cell,
182  & last_cell_grid_num_find_adj_cell,
183  & num_srch_cells_find_adj_cell,
184  & srch_corners_find_adj_cell

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ last_cell_grid_num_find_adj_cell

integer (scrip_i4), save scrip_remap_conservative::last_cell_grid_num_find_adj_cell

Definition at line 180 of file scrip_remap_conservative.f.

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ last_cell_grid_num_get_srch_cells

integer (scrip_i4), save scrip_remap_conservative::last_cell_grid_num_get_srch_cells

Definition at line 166 of file scrip_remap_conservative.f.

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ last_cell_grid_num_locate_point

integer (scrip_i4), save scrip_remap_conservative::last_cell_grid_num_locate_point =0

Definition at line 124 of file scrip_remap_conservative.f.

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ last_cell_grid_num_locate_segstart

integer (scrip_i4), save scrip_remap_conservative::last_cell_grid_num_locate_segstart =0

Definition at line 96 of file scrip_remap_conservative.f.

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ last_cell_locate_point

integer (scrip_i4), save scrip_remap_conservative::last_cell_locate_point =0

Definition at line 124 of file scrip_remap_conservative.f.

124  integer (SCRIP_i4), save ::
125  & last_cell_locate_point=0, ! save the search parameters
126  & last_cell_grid_num_locate_point=0, ! if unchanged, reuse
127  ! search lists
128  & last_srch_grid_num_locate_point=0

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ last_cell_locate_segstart

integer (scrip_i4), save scrip_remap_conservative::last_cell_locate_segstart =0

Definition at line 96 of file scrip_remap_conservative.f.

96  integer (SCRIP_i4), save ::
97  & last_cell_locate_segstart=0, ! save the search parameters
98  & last_cell_grid_num_locate_segstart=0, ! if unchanged, reuse
99  ! search lists
100  & last_srch_grid_num_locate_segstart=0

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ last_srch_grid_num_get_srch_cells

integer (scrip_i4), save scrip_remap_conservative::last_srch_grid_num_get_srch_cells

Definition at line 166 of file scrip_remap_conservative.f.

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ last_srch_grid_num_locate_point

integer (scrip_i4), save scrip_remap_conservative::last_srch_grid_num_locate_point =0

Definition at line 124 of file scrip_remap_conservative.f.

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ last_srch_grid_num_locate_segstart

integer (scrip_i4), save scrip_remap_conservative::last_srch_grid_num_locate_segstart =0

Definition at line 96 of file scrip_remap_conservative.f.

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ link_add1

integer (scrip_i4), dimension(:,:), allocatable, save scrip_remap_conservative::link_add1

Definition at line 86 of file scrip_remap_conservative.f.

86  integer (SCRIP_i4), dimension(:,:), allocatable, save ::
87  & link_add1, ! min,max link add to restrict search
88  & link_add2 ! min,max link add to restrict search

Referenced by scrip_interface::scrip_clear(), and store_link_cnsrv().

◆ link_add2

integer (scrip_i4), dimension(:,:), allocatable, save scrip_remap_conservative::link_add2

Definition at line 86 of file scrip_remap_conservative.f.

Referenced by scrip_interface::scrip_clear(), and store_link_cnsrv().

◆ nthreads

integer (scrip_i4) scrip_remap_conservative::nthreads =2

Definition at line 73 of file scrip_remap_conservative.f.

73  integer (SCRIP_i4) :: nthreads=2 ! Number of parallel threads

Referenced by scrip_interface::scrip().

◆ num_srch_cell_locate_points

integer (scrip_i4), save scrip_remap_conservative::num_srch_cell_locate_points =0

Definition at line 130 of file scrip_remap_conservative.f.

130  integer (SCRIP_i4), save ::
131  & num_srch_cell_locate_points=0,
132  & srch_corners_locate_point ! number of corners for each cell

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ num_srch_cells_find_adj_cell

integer (scrip_i4), save scrip_remap_conservative::num_srch_cells_find_adj_cell

Definition at line 180 of file scrip_remap_conservative.f.

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ num_srch_cells_loc_get_srch_cells

integer (scrip_i4), save scrip_remap_conservative::num_srch_cells_loc_get_srch_cells

Definition at line 148 of file scrip_remap_conservative.f.

148  integer (SCRIP_i4), save ::
149  & num_srch_cells_loc_get_srch_cells, ! Number of srch cells
150  ! found
151  & srch_corners_loc_get_srch_cells ! Number of corners for

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ num_srch_cells_locate_segstart

integer (scrip_i4), save scrip_remap_conservative::num_srch_cells_locate_segstart =0

Definition at line 102 of file scrip_remap_conservative.f.

102  integer (SCRIP_i4), save ::
103  & num_srch_cells_locate_segstart=0,
104  & srch_corners_locate_segstart ! number of corners for

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ srch_add_find_adj_cell

integer (scrip_i4), dimension(:), allocatable, save scrip_remap_conservative::srch_add_find_adj_cell

Definition at line 186 of file scrip_remap_conservative.f.

186  integer (SCRIP_i4), dimension(:), allocatable, save ::
187  & srch_add_find_adj_cell

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ srch_add_loc_get_srch_cells

integer (scrip_i4), dimension(:), allocatable, save scrip_remap_conservative::srch_add_loc_get_srch_cells

Definition at line 154 of file scrip_remap_conservative.f.

154  integer (SCRIP_i4), dimension(:), allocatable, save ::
155  & srch_add_loc_get_srch_cells ! Global addresses of

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ srch_add_locate_point

integer (scrip_i4), dimension(:), allocatable, save scrip_remap_conservative::srch_add_locate_point

Definition at line 134 of file scrip_remap_conservative.f.

134  integer (SCRIP_i4), dimension(:), allocatable, save ::
135  & srch_add_locate_point ! global address of cells in

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ srch_add_locate_segstart

integer (scrip_i4), dimension(:), allocatable, save scrip_remap_conservative::srch_add_locate_segstart

Definition at line 107 of file scrip_remap_conservative.f.

107  integer (SCRIP_i4), dimension(:), allocatable, save ::
108  & srch_add_locate_segstart ! global address of cells

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ srch_center_lat_find_adj_cell

real (scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lat_find_adj_cell

Definition at line 191 of file scrip_remap_conservative.f.

191  real (SCRIP_r8), dimension(:), allocatable, save ::
192  & srch_center_lat_find_adj_cell, srch_center_lon_find_adj_cell

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ srch_center_lat_loc_get_srch_cells

real (scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lat_loc_get_srch_cells

Definition at line 162 of file scrip_remap_conservative.f.

162  real (SCRIP_r8), dimension(:), allocatable, save ::
163  & srch_center_lat_loc_get_srch_cells,
164  & srch_center_lon_loc_get_srch_cells

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ srch_center_lat_locate_point

real (scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lat_locate_point

Definition at line 144 of file scrip_remap_conservative.f.

144  real (SCRIP_r8), dimension(:), allocatable, save ::
145  & srch_center_lat_locate_point, ! lat of center of srch cells
146  & srch_center_lon_locate_point ! lon of center of srch cells

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ srch_center_lat_locate_segstart

real(scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lat_locate_segstart

Definition at line 117 of file scrip_remap_conservative.f.

117  real(SCRIP_r8), dimension(:), allocatable, save ::
118  & srch_center_lat_locate_segstart,! lat of center of srch cells
119  & srch_center_lon_locate_segstart ! lon of center of srch cells

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ srch_center_lon_find_adj_cell

real (scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lon_find_adj_cell

Definition at line 191 of file scrip_remap_conservative.f.

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ srch_center_lon_loc_get_srch_cells

real (scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lon_loc_get_srch_cells

Definition at line 162 of file scrip_remap_conservative.f.

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ srch_center_lon_locate_point

real (scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lon_locate_point

Definition at line 144 of file scrip_remap_conservative.f.

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ srch_center_lon_locate_segstart

real(scrip_r8), dimension(:), allocatable, save scrip_remap_conservative::srch_center_lon_locate_segstart

Definition at line 117 of file scrip_remap_conservative.f.

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ srch_corner_lat_find_adj_cell

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lat_find_adj_cell

Definition at line 188 of file scrip_remap_conservative.f.

188  real (SCRIP_r8), dimension(:,:), allocatable, save ::
189  & srch_corner_lat_find_adj_cell, srch_corner_lon_find_adj_cell

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ srch_corner_lat_loc_get_srch_cells

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lat_loc_get_srch_cells

Definition at line 158 of file scrip_remap_conservative.f.

158  real (SCRIP_r8), dimension(:,:), allocatable, save ::
159  & srch_corner_lat_loc_get_srch_cells,
160  & srch_corner_lon_loc_get_srch_cells

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ srch_corner_lat_locate_point

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lat_locate_point

Definition at line 138 of file scrip_remap_conservative.f.

138  real (SCRIP_r8), dimension(:,:), allocatable, save ::
139  & srch_corner_lat_locate_point, ! lat of each corner of srch
140  ! cells
141  & srch_corner_lon_locate_point ! lon of each corner of srch

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ srch_corner_lat_locate_segstart

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lat_locate_segstart

Definition at line 111 of file scrip_remap_conservative.f.

111  real (SCRIP_r8), dimension(:,:), allocatable, save ::
112  & srch_corner_lat_locate_segstart, ! lat of each corner of
113  ! srch cells
114  & srch_corner_lon_locate_segstart ! lon of each corner of

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ srch_corner_lon_find_adj_cell

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lon_find_adj_cell

Definition at line 188 of file scrip_remap_conservative.f.

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ srch_corner_lon_loc_get_srch_cells

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lon_loc_get_srch_cells

Definition at line 158 of file scrip_remap_conservative.f.

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ srch_corner_lon_locate_point

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lon_locate_point

Definition at line 138 of file scrip_remap_conservative.f.

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ srch_corner_lon_locate_segstart

real (scrip_r8), dimension(:,:), allocatable, save scrip_remap_conservative::srch_corner_lon_locate_segstart

Definition at line 111 of file scrip_remap_conservative.f.

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

◆ srch_corners_find_adj_cell

integer (scrip_i4), save scrip_remap_conservative::srch_corners_find_adj_cell

Definition at line 180 of file scrip_remap_conservative.f.

Referenced by find_adj_cell(), and scrip_interface::scrip_clear().

◆ srch_corners_loc_get_srch_cells

integer (scrip_i4), save scrip_remap_conservative::srch_corners_loc_get_srch_cells

Definition at line 148 of file scrip_remap_conservative.f.

Referenced by get_srch_cells(), and scrip_interface::scrip_clear().

◆ srch_corners_locate_point

integer (scrip_i4), save scrip_remap_conservative::srch_corners_locate_point

Definition at line 130 of file scrip_remap_conservative.f.

Referenced by locate_point(), and scrip_interface::scrip_clear().

◆ srch_corners_locate_segstart

integer (scrip_i4), save scrip_remap_conservative::srch_corners_locate_segstart

Definition at line 102 of file scrip_remap_conservative.f.

Referenced by locate_segstart(), and scrip_interface::scrip_clear().

scrip_remap_conservative::srch_corners_locate_point
integer(scrip_i4), save srch_corners_locate_point
Definition: scrip_remap_conservative.f:130
scrip_remap_conservative::avoid_pole_offset
real(scrip_r8), save avoid_pole_offset
Definition: scrip_remap_conservative.f:83
scrip_grids::grid2_centroid_lat
real(scrip_r8), dimension(:), allocatable, target, save grid2_centroid_lat
Definition: scrip_grids.f:103
scrip_remap_vars::norm_opt_none
integer(scrip_i4), parameter norm_opt_none
Definition: scrip_remap_vars.f:54
scrip_grids::grid1_corners
integer(scrip_i4), save grid1_corners
Definition: scrip_grids.f:68
scrip_remap_vars::grid1_add_map1
integer(scrip_i4), dimension(:), allocatable, save grid1_add_map1
Definition: scrip_remap_vars.f:77
scrip_constants::half
real(kind=scrip_r8), parameter half
Definition: scrip_constants.f:50
constants::pi
real, parameter pi
PI Value of Pi.
Definition: constants.F90:71
scrip_grids::grid1_mask
logical(scrip_logical), dimension(:), allocatable, target, save grid1_mask
Definition: scrip_grids.f:93
scrip_remap_conservative::srch_corners_locate_segstart
integer(scrip_i4), save srch_corners_locate_segstart
Definition: scrip_remap_conservative.f:102
scrip_remap_conservative::srch_center_lat_locate_point
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_locate_point
Definition: scrip_remap_conservative.f:144
scrip_remap_conservative::srch_corner_lon_locate_segstart
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_locate_segstart
Definition: scrip_remap_conservative.f:111
scrip_remap_conservative::srch_center_lon_locate_segstart
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_locate_segstart
Definition: scrip_remap_conservative.f:117
scrip_remap_conservative::last_cell_grid_num_get_srch_cells
integer(scrip_i4), save last_cell_grid_num_get_srch_cells
Definition: scrip_remap_conservative.f:166
scrip_grids::grid1_centroid_lat
real(scrip_r8), dimension(:), allocatable, target, save grid1_centroid_lat
Definition: scrip_grids.f:103
scrip_remap_conservative::find_adj_cell
subroutine find_adj_cell(cell_add, edge_id, cell_grid_num, adj_add)
Definition: scrip_remap_conservative.f:5987
scrip_grids::grid1_spole_cell
integer(scrip_i4), save grid1_spole_cell
Definition: scrip_grids.f:135
scrip_remap_vars::wt_lowest
real(kind=scrip_r8) wt_lowest
Definition: scrip_remap_vars.f:88
scrip_remap_vars::num_maps
integer(scrip_i4), save num_maps
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::last_cell_locate_segstart
integer(scrip_i4), save last_cell_locate_segstart
Definition: scrip_remap_conservative.f:96
scrip_remap_vars::wt_highest
real(kind=scrip_r8) wt_highest
Definition: scrip_remap_vars.f:88
scrip_remap_conservative::converge_to_bdry
subroutine converge_to_bdry(cell_add, cell_grid_num, ncorners, cell_corner_lat, cell_corner_lon, cell_center_lat, cell_center_lon, inpt_lat, inpt_lon, outpt_lat, outpt_lon, bpt_lat, bpt_lon, bedgeid)
Definition: scrip_remap_conservative.f:6217
scrip_grids::luse_grid2_area
logical(scrip_logical), save luse_grid2_area
Definition: scrip_grids.f:126
scrip_remap_conservative::first_call_locate_point
logical(scrip_logical), save first_call_locate_point
Definition: scrip_remap_conservative.f:121
scrip_remap_conservative::link_add2
integer(scrip_i4), dimension(:,:), allocatable, save link_add2
Definition: scrip_remap_conservative.f:86
scrip_timers::timer_start
subroutine timer_start(timer)
Definition: scrip_timers.f:216
scrip_grids::grid2_centroid_lon
real(scrip_r8), dimension(:), allocatable, target, save grid2_centroid_lon
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_corner_lon_loc_get_srch_cells
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_loc_get_srch_cells
Definition: scrip_remap_conservative.f:158
scrip_grids::num_srch_bins
integer(scrip_i4), save num_srch_bins
Definition: scrip_grids.f:152
scrip_remap_conservative::ptincell
subroutine ptincell(ptlat, ptlon, cell_add, ncorners, cell_corner_lat, cell_corner_lon, cell_center_lat, cell_center_lon, cell_grid_id, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4608
scrip_remap_conservative::last_srch_grid_num_locate_segstart
integer(scrip_i4), save last_srch_grid_num_locate_segstart
Definition: scrip_remap_conservative.f:96
scrip_remap_conservative::store_link_cnsrv
subroutine store_link_cnsrv(add1, add2, weights)
Definition: scrip_remap_conservative.f:3951
scrip_remap_conservative::last_cell_grid_num_locate_point
integer(scrip_i4), save last_cell_grid_num_locate_point
Definition: scrip_remap_conservative.f:124
scrip_remap_vars::num_links_map1
integer(scrip_i4), save num_links_map1
Definition: scrip_remap_vars.f:66
scrip_timers::timer_print
subroutine timer_print(timer)
Definition: scrip_timers.f:176
scrip_grids::grid2_mask
logical(scrip_logical), dimension(:), allocatable, target, save grid2_mask
Definition: scrip_grids.f:93
scrip_remap_vars::norm_opt
integer(scrip_i4), save norm_opt
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::cell_integrate
subroutine cell_integrate(cell_add, grid_num, phi_or_theta)
Definition: scrip_remap_conservative.f:1005
scrip_grids::grid2_center_lat
real(scrip_r8), dimension(:), allocatable, target, save grid2_center_lat
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_add_find_adj_cell
integer(scrip_i4), dimension(:), allocatable, save srch_add_find_adj_cell
Definition: scrip_remap_conservative.f:186
scrip_remap_vars::max_links_map2
integer(scrip_i4), save max_links_map2
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::last_cell_locate_point
integer(scrip_i4), save last_cell_locate_point
Definition: scrip_remap_conservative.f:124
scrip_remap_conservative::first_call_locate_segstart
logical(scrip_logical), save first_call_locate_segstart
Definition: scrip_remap_conservative.f:93
scrip_remap_conservative::srch_center_lon_locate_point
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_locate_point
Definition: scrip_remap_conservative.f:144
m_constants::pih
real pih
pi/2
Definition: mod_constants.f90:31
scrip_constants::one
real(kind=scrip_r8), parameter one
Definition: scrip_constants.f:50
scrip_remap_vars::num_links_map2
integer(scrip_i4), save num_links_map2
Definition: scrip_remap_vars.f:66
scrip_grids::grid1_area
real(scrip_r8), dimension(:), allocatable, target, save grid1_area
Definition: scrip_grids.f:103
scrip_grids::north_thresh
real(scrip_r8), save north_thresh
Definition: scrip_grids.f:169
scrip_grids::grid1_npole_cell
integer(scrip_i4), save grid1_npole_cell
Definition: scrip_grids.f:135
scrip_remap_conservative::srch_corner_lat_locate_point
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_locate_point
Definition: scrip_remap_conservative.f:138
scrip_remap_conservative::get_srch_cells
subroutine get_srch_cells(cell_add, cell_grid_num, srch_grid_num, num_srch_cells, srch_add, srch_corners, srch_corner_lat, srch_corner_lon, srch_center_lat, srch_center_lon)
Definition: scrip_remap_conservative.f:5557
scrip_constants::tiny
real(kind=scrip_r8), parameter tiny
Definition: scrip_constants.f:50
scrip_remap_conservative::srch_corner_lon_find_adj_cell
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_find_adj_cell
Definition: scrip_remap_conservative.f:188
scrip_grids::special_polar_cell2
logical(scrip_logical), dimension(:), allocatable, target, save special_polar_cell2
Definition: scrip_grids.f:93
scrip_remap_conservative::last_srch_grid_num_get_srch_cells
integer(scrip_i4), save last_srch_grid_num_get_srch_cells
Definition: scrip_remap_conservative.f:166
scrip_grids::grid2_corners
integer(scrip_i4), save grid2_corners
Definition: scrip_grids.f:68
scrip_grids::grid1_area_in
real(scrip_r8), dimension(:), allocatable, target, save grid1_area_in
Definition: scrip_grids.f:103
scrip_remap_conservative::num_srch_cell_locate_points
integer(scrip_i4), save num_srch_cell_locate_points
Definition: scrip_remap_conservative.f:130
scrip_constants::zero
real(kind=scrip_r8), parameter zero
Definition: scrip_constants.f:50
scrip_kindsmod::scrip_r8
integer, parameter, public scrip_r8
Definition: scrip_kindsmod.f90:38
scrip_remap_conservative::last_cell_find_adj_cell
integer(scrip_i4), save last_cell_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_constants::quart
real(kind=scrip_r8), parameter quart
Definition: scrip_constants.f:50
scrip_grids::grid1_centroid_lon
real(scrip_r8), dimension(:), allocatable, target, save grid1_centroid_lon
Definition: scrip_grids.f:103
scrip_grids::grid1_bound_box
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_bound_box
Definition: scrip_grids.f:131
scrip_remap_vars::grid2_add_map2
integer(scrip_i4), dimension(:), allocatable, save grid2_add_map2
Definition: scrip_remap_vars.f:77
scrip_grids::grid2_corner_lat
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_corner_lat
Definition: scrip_grids.f:120
scrip_remap_vars::norm_opt_frcarea
integer(scrip_i4), parameter norm_opt_frcarea
Definition: scrip_remap_vars.f:54
scrip_remap_vars::max_links_map1
integer(scrip_i4), save max_links_map1
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::num_srch_cells_locate_segstart
integer(scrip_i4), save num_srch_cells_locate_segstart
Definition: scrip_remap_conservative.f:102
scrip_remap_conservative::srch_center_lat_find_adj_cell
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_find_adj_cell
Definition: scrip_remap_conservative.f:191
scrip_grids::luse_grid1_area
logical(scrip_logical), save luse_grid1_area
Definition: scrip_grids.f:126
scrip_grids::grid1_center_lon
real(scrip_r8), dimension(:), allocatable, target, save grid1_center_lon
Definition: scrip_grids.f:103
scrip_constants::two
real(kind=scrip_r8), parameter two
Definition: scrip_constants.f:50
scrip_grids::grid2_npole_cell
integer(scrip_i4), save grid2_npole_cell
Definition: scrip_grids.f:135
scrip_grids::grid2_area
real(scrip_r8), dimension(:), allocatable, target, save grid2_area
Definition: scrip_grids.f:103
scrip_remap_conservative::intersection
subroutine intersection(seg_cell_id, seg_grid_id, beglat, beglon, endlat, endlon, begseg, begedge, cell_id, ncorners, cell_corner_lat, cell_corner_lon, cell_grid_id, intrsct_lat, intrsct_lon, intedge, sinang2, lcoinc, lthresh)
Definition: scrip_remap_conservative.f:2196
scrip_grids::grid1_corner_lat
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_corner_lat
Definition: scrip_grids.f:120
scrip_remap_conservative::srch_center_lon_loc_get_srch_cells
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_loc_get_srch_cells
Definition: scrip_remap_conservative.f:162
scrip_remap_vars::frac_highest
real(kind=scrip_r8) frac_highest
Definition: scrip_remap_vars.f:87
scrip_remap_conservative::srch_corner_lat_find_adj_cell
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_find_adj_cell
Definition: scrip_remap_conservative.f:188
scrip_grids::grid2_size
integer(scrip_i4), save grid2_size
Definition: scrip_grids.f:68
scrip_grids::bin_addr2
integer(scrip_i4), dimension(:,:), allocatable, save bin_addr2
Definition: scrip_grids.f:155
scrip_remap_vars::frac_lowest
real(kind=scrip_r8) frac_lowest
Definition: scrip_remap_vars.f:87
scrip_remap_conservative::last_cell_grid_num_locate_segstart
integer(scrip_i4), save last_cell_grid_num_locate_segstart
Definition: scrip_remap_conservative.f:96
scrip_remap_conservative::ptinpoly
subroutine ptinpoly(ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4815
scrip_remap_conservative::last_cell_add_get_srch_cells
integer(scrip_i4), save last_cell_add_get_srch_cells
Definition: scrip_remap_conservative.f:166
scrip_remap_vars::wts_map1
real(scrip_r8), dimension(:,:), allocatable, save wts_map1
Definition: scrip_remap_vars.f:83
scrip_remap_conservative::link_add1
integer(scrip_i4), dimension(:,:), allocatable, save link_add1
Definition: scrip_remap_conservative.f:86
scrip_remap_conservative::srch_add_locate_point
integer(scrip_i4), dimension(:), allocatable, save srch_add_locate_point
Definition: scrip_remap_conservative.f:134
scrip_remap_vars::resize_remap_vars
subroutine resize_remap_vars(nmap, increment)
Definition: scrip_remap_vars.f:173
scrip_remap_conservative::line_integral
subroutine line_integral(phi_or_theta, weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
Definition: scrip_remap_conservative.f:3460
scrip_remap_conservative::last_cell_grid_num_find_adj_cell
integer(scrip_i4), save last_cell_grid_num_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_grids::grid2_area_in
real(scrip_r8), dimension(:), allocatable, target, save grid2_area_in
Definition: scrip_grids.f:103
scrip_remap_conservative::num_srch_cells_loc_get_srch_cells
integer(scrip_i4), save num_srch_cells_loc_get_srch_cells
Definition: scrip_remap_conservative.f:148
scrip_grids::grid1_center_lat
real(scrip_r8), dimension(:), allocatable, target, save grid1_center_lat
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_corner_lat_loc_get_srch_cells
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_loc_get_srch_cells
Definition: scrip_remap_conservative.f:158
scrip_remap_conservative::num_srch_cells_find_adj_cell
integer(scrip_i4), save num_srch_cells_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_remap_conservative::modify_polar_cell
subroutine modify_polar_cell(ncorners, nalloc, cell_corner_lat, cell_corner_lon)
Definition: scrip_remap_conservative.f:2082
scrip_remap_conservative::ptinpolarpoly
subroutine ptinpolarpoly(ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, whichpole, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:5042
scrip_grids::grid1_corner_lon
real(scrip_r8), dimension(:,:), allocatable, target, save grid1_corner_lon
Definition: scrip_grids.f:120
scrip_remap_conservative::srch_corner_lon_locate_point
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lon_locate_point
Definition: scrip_remap_conservative.f:138
scrip_grids::npseg
integer(scrip_i4), save npseg
Definition: scrip_grids.f:178
scrip_remap_vars::grid1_add_map2
integer(scrip_i4), dimension(:), allocatable, save grid1_add_map2
Definition: scrip_remap_vars.f:77
scrip_remap_vars::grid2_add_map1
integer(scrip_i4), dimension(:), allocatable, save grid2_add_map1
Definition: scrip_remap_vars.f:77
scrip_remap_conservative::pole_intersection
subroutine pole_intersection(location, ncorners, cell_corners_lat, cell_corners_lon, cell_grid_id, beglat, beglon, endlat, endlon, begseg, begedge, intedge, intrsct_lat, intrsct_lon, sinang2, lcoinc, lthresh)
Definition: scrip_remap_conservative.f:2728
scrip_grids::south_thresh
real(scrip_r8), save south_thresh
Definition: scrip_grids.f:169
scrip_remap_conservative::srch_corner_lat_locate_segstart
real(scrip_r8), dimension(:,:), allocatable, save srch_corner_lat_locate_segstart
Definition: scrip_remap_conservative.f:111
scrip_grids::grid1_size
integer(scrip_i4), save grid1_size
Definition: scrip_grids.f:68
scrip_remap_conservative::first_call_find_adj_cell
logical(scrip_logical), save first_call_find_adj_cell
Definition: scrip_remap_conservative.f:174
scrip_remap_conservative::srch_center_lon_find_adj_cell
real(scrip_r8), dimension(:), allocatable, save srch_center_lon_find_adj_cell
Definition: scrip_remap_conservative.f:191
scrip_remap_vars::wts_map2
real(scrip_r8), dimension(:,:), allocatable, save wts_map2
Definition: scrip_remap_vars.f:83
scrip_grids::grid2_bound_box
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_bound_box
Definition: scrip_grids.f:131
scrip_grids::grid2_corner_lon
real(scrip_r8), dimension(:,:), allocatable, target, save grid2_corner_lon
Definition: scrip_grids.f:120
scrip_grids::grid1_frac
real(scrip_r8), dimension(:), allocatable, target, save grid1_frac
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_corners_find_adj_cell
integer(scrip_i4), save srch_corners_find_adj_cell
Definition: scrip_remap_conservative.f:180
scrip_remap_conservative::srch_center_lat_loc_get_srch_cells
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_loc_get_srch_cells
Definition: scrip_remap_conservative.f:162
scrip_remap_conservative::srch_corners_loc_get_srch_cells
integer(scrip_i4), save srch_corners_loc_get_srch_cells
Definition: scrip_remap_conservative.f:148
scrip_remap_conservative::srch_add_loc_get_srch_cells
integer(scrip_i4), dimension(:), allocatable, save srch_add_loc_get_srch_cells
Definition: scrip_remap_conservative.f:154
scrip_remap_conservative::line_integral_theta
subroutine line_integral_theta(weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
Definition: scrip_remap_conservative.f:3682
scrip_remap_conservative::line_integral_phi
subroutine line_integral_phi(weights, num_wts, in_phi1, in_phi2, theta1, theta2, grid1_lat, grid1_lon, grid2_lat, grid2_lon)
Definition: scrip_remap_conservative.f:3522
scrip_remap_vars::norm_opt_dstarea
integer(scrip_i4), parameter norm_opt_dstarea
Definition: scrip_remap_vars.f:54
scrip_remap_conservative::ptinpolygen2
subroutine ptinpolygen2(ptx, pty, ncorners, cell_corner_x, cell_corner_y, latlon, inpoly, lboundary, edgeid)
Definition: scrip_remap_conservative.f:5378
scrip_remap_conservative::srch_add_locate_segstart
integer(scrip_i4), dimension(:), allocatable, save srch_add_locate_segstart
Definition: scrip_remap_conservative.f:107
scrip_grids::grid2_spole_cell
integer(scrip_i4), save grid2_spole_cell
Definition: scrip_grids.f:135
scrip_remap_conservative::avoid_pole_count
integer(scrip_i4), save avoid_pole_count
Definition: scrip_remap_conservative.f:80
scrip_remap_conservative::last_srch_grid_num_locate_point
integer(scrip_i4), save last_srch_grid_num_locate_point
Definition: scrip_remap_conservative.f:124
w3sic3md::delta
real function, dimension(:), allocatable delta(X)
Definition: w3sic3md.F90:1733
scrip_remap_conservative::first_call_store_link_cnsrv
logical(scrip_logical), save first_call_store_link_cnsrv
Definition: scrip_remap_conservative.f:90
scrip_grids::special_polar_cell1
logical(scrip_logical), dimension(:), allocatable, target, save special_polar_cell1
Definition: scrip_grids.f:93
scrip_remap_vars::num_wts
integer(scrip_i4), save num_wts
Definition: scrip_remap_vars.f:66
scrip_remap_conservative::first_call_get_srch_cells
logical(scrip_logical), save first_call_get_srch_cells
Definition: scrip_remap_conservative.f:171
scrip_timers::timer_stop
subroutine timer_stop(timer)
Definition: scrip_timers.f:250
scrip_grids::grid2_frac
real(scrip_r8), dimension(:), allocatable, target, save grid2_frac
Definition: scrip_grids.f:103
m_constants::pi2
real pi2
2*pi
Definition: mod_constants.f90:30
scrip_remap_conservative::locate_point
subroutine locate_point(ptlat, ptlon, cell, cell_grid_num, srch_grid_num, cont_cell, lboundary, edgeid)
Definition: scrip_remap_conservative.f:4423
scrip_grids::bin_addr1
integer(scrip_i4), dimension(:,:), allocatable, save bin_addr1
Definition: scrip_grids.f:155
scrip_grids::grid2_center_lon
real(scrip_r8), dimension(:), allocatable, target, save grid2_center_lon
Definition: scrip_grids.f:103
scrip_remap_conservative::srch_center_lat_locate_segstart
real(scrip_r8), dimension(:), allocatable, save srch_center_lat_locate_segstart
Definition: scrip_remap_conservative.f:117
scrip_remap_vars::resize_increment
integer(scrip_i4), save resize_increment
Definition: scrip_remap_vars.f:66