|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 C**** 00002 C ************************ 00003 C * OASIS MODULE * 00004 C * ------------ * 00005 C ************************ 00006 C**** 00007 C*********************************************************************** 00008 C This module belongs to the SCRIP library. It is modified to run 00009 C within OASIS. 00010 C Modifications: 00011 C - routine does not read SCRIP grid description files, but gets 00012 C arrays from the calling routine 00013 C - unit conversion : only from degrees (OASIS unit) to radians 00014 C - some allocated array will be freed in the end to allow multiple 00015 C calls of SCRIP 00016 C - map-methods and restriction-types are written in capital letters 00017 C - added bin definition for reduced grid 00018 C 00019 C Modified by V. Gayler, M&D 20.09.2001 00020 C Modified by D. Declat, CERFACS 27.06.2002 00021 C*********************************************************************** 00022 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00023 ! 00024 ! This module reads in and initializes two grids for remapping. 00025 ! NOTE: grid1 must be the master grid -- the grid that determines 00026 ! which cells participate (e.g. land mask) and the fractional 00027 ! area of grid2 cells that participate in the remapping. 00028 ! 00029 !----------------------------------------------------------------------- 00030 ! 00031 ! CVS:$Id: grids.f 2826 2010-12-10 11:14:21Z valcke $ 00032 ! 00033 ! Copyright (c) 1997, 1998 the Regents of the University of 00034 ! California. 00035 ! 00036 ! This software and ancillary information (herein called software) 00037 ! called SCRIP is made available under the terms described here. 00038 ! The software has been approved for release with associated 00039 ! LA-CC Number 98-45. 00040 ! 00041 ! Unless otherwise indicated, this software has been authored 00042 ! by an employee or employees of the University of California, 00043 ! operator of the Los Alamos National Laboratory under Contract 00044 ! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S. 00045 ! Government has rights to use, reproduce, and distribute this 00046 ! software. The public may copy and use this software without 00047 ! charge, provided that this Notice and any statement of authorship 00048 ! are reproduced on all copies. Neither the Government nor the 00049 ! University makes any warranty, express or implied, or assumes 00050 ! any liability or responsibility for the use of this software. 00051 ! 00052 ! If software is modified to produce derivative works, such modified 00053 ! software should be clearly marked, so as not to confuse it with 00054 ! the version available from Los Alamos National Laboratory. 00055 ! 00056 !*********************************************************************** 00057 00058 module grids 00059 00060 !----------------------------------------------------------------------- 00061 00062 use kinds_mod ! defines data types 00063 use constants ! common constants 00064 use iounits ! I/O unit manager 00065 USE mod_unit 00066 USE mod_printing 00067 00068 implicit none 00069 00070 !----------------------------------------------------------------------- 00071 ! 00072 ! variables that describe each grid 00073 ! 00074 !----------------------------------------------------------------------- 00075 00076 integer (kind=int_kind), save :: 00077 & grid1_size, grid2_size, ! total points on each grid 00078 & grid1_rank, grid2_rank, ! rank of each grid 00079 & grid1_corners, grid2_corners ! number of corners 00080 ! for each grid cell 00081 00082 integer (kind=int_kind), dimension(:), allocatable, save :: 00083 & grid1_dims, grid2_dims ! size of each grid dimension 00084 00085 character(char_len), save :: 00086 & grid1_name, grid2_name ! name for each grid 00087 00088 character (char_len), save :: 00089 & grid1_units, ! units for grid coords (degs/radians) 00090 & grid2_units ! units for grid coords 00091 00092 real (kind=dbl_kind), parameter :: 00093 & deg2rad = pi/180. ! conversion for deg to rads 00094 00095 !----------------------------------------------------------------------- 00096 ! 00097 ! grid coordinates and masks 00098 ! 00099 !----------------------------------------------------------------------- 00100 00101 logical (kind=log_kind), dimension(:), allocatable, save :: 00102 & grid1_mask, ! flag which cells participate 00103 & grid2_mask ! flag which cells participate 00104 00105 real (kind=dbl_kind), dimension(:), allocatable, save :: 00106 & grid1_center_lat, ! lat/lon coordinates for 00107 & grid1_center_lon, ! each grid center in radians 00108 & grid2_center_lat, 00109 & grid2_center_lon, 00110 & grid1_area, ! tot area of each grid1 cell 00111 & grid2_area, ! tot area of each grid2 cell 00112 & grid1_area_in, ! area of grid1 cell from file 00113 & grid2_area_in, ! area of grid2 cell from file 00114 & grid1_frac, ! fractional area of grid cells 00115 & grid2_frac ! participating in remapping 00116 00117 real (kind=dbl_kind), dimension(:,:), allocatable, save :: 00118 & grid1_corner_lat, ! lat/lon coordinates for 00119 & grid1_corner_lon, ! each grid corner in radians 00120 & grid2_corner_lat, 00121 & grid2_corner_lon 00122 00123 logical (kind=log_kind), save :: 00124 & luse_grid_centers ! use centers for bounding boxes 00125 &, luse_grid1_area ! use area from grid file 00126 &, luse_grid2_area ! use area from grid file 00127 00128 real (kind=dbl_kind), dimension(:,:), allocatable, save :: 00129 & grid1_bound_box, ! lat/lon bounding box for use 00130 & grid2_bound_box ! in restricting grid searches 00131 00132 !----------------------------------------------------------------------- 00133 ! 00134 ! bins for restricting searches 00135 ! 00136 !----------------------------------------------------------------------- 00137 00138 character (char_len), save :: 00139 & restrict_type ! type of bins to use 00140 00141 integer (kind=int_kind), save :: 00142 & num_srch_bins, ! num of bins for restricted srch 00143 & num_srch_red ! num of bins for reduced case 00144 00145 integer (kind=int_kind), dimension(:,:), allocatable, save :: 00146 & bin_addr1, ! min,max adds for grid1 cells in this lat bin 00147 & bin_addr2 ! min,max adds for grid2 cells in this lat bin 00148 integer (kind=int_kind), dimension(:,:), allocatable, save :: 00149 & bin_addr1_r ! min,max adds for red grid1 cells 00150 00151 real(kind=dbl_kind), dimension(:,:), allocatable, save :: 00152 & bin_lats ! min,max latitude for each search bin 00153 &, bin_lons ! min,max longitude for each search bin 00154 real(kind=dbl_kind), dimension(:,:), allocatable, save :: 00155 & bin_lats_r ! min,max lat for each search bin for red grid 00156 &, bin_lons_r ! min,max lon for each search bin for red grid 00157 00158 !*********************************************************************** 00159 00160 contains 00161 00162 !*********************************************************************** 00163 00164 subroutine grid_init(m_method, rst_type, n_srch_bins, 00165 $ src_size, dst_size, src_dims, dst_dims, 00166 & src_rank, dst_rank, ncrn_src, ncrn_dst, 00167 & src_mask, dst_mask, src_name, dst_name, 00168 & src_lat, src_lon, dst_lat, dst_lon, 00169 & src_corner_lat, src_corner_lon, 00170 & dst_corner_lat, dst_corner_lon) 00171 00172 !----------------------------------------------------------------------- 00173 ! 00174 ! this routine gets grid info from routine scriprmp and makes any 00175 ! necessary changes (e.g. for 0,2pi longitude range) 00176 ! 00177 !----------------------------------------------------------------------- 00178 00179 !----------------------------------------------------------------------- 00180 ! 00181 ! input variables 00182 ! 00183 !----------------------------------------------------------------------- 00184 00185 integer (kind=int_kind), intent(in) :: 00186 & n_srch_bins, ! num of bins for restricted srch 00187 & src_size, ! source grid size 00188 & dst_size, ! target grid size 00189 & src_rank, ! source grid rank 00190 & dst_rank, ! target grid rank 00191 & src_dims(src_rank), ! source grid dimensions 00192 & dst_dims(dst_rank), ! target grid dimensions 00193 & ncrn_src, ! number of corners of a source grid cell 00194 & ncrn_dst, ! number of corners of a target grid cell 00195 & src_mask(src_size), ! source grid mask (master mask) 00196 & dst_mask(dst_size) ! target grid mask 00197 00198 character*8, intent(in) :: 00199 & m_method, ! remapping method 00200 & rst_type, ! restriction type 00201 & src_name, ! source grid name 00202 & dst_name ! target grid name 00203 00204 real (kind=real_kind), intent (in) :: 00205 & src_lat(src_size), ! source grid latitudes 00206 & src_lon(src_size), ! sourde grid longitudes 00207 & dst_lat(dst_size), ! target grid latitudes 00208 & dst_lon(dst_size), ! target grid longitudes 00209 & src_corner_lat(ncrn_src,src_size), 00210 & src_corner_lon(ncrn_src,src_size), 00211 & dst_corner_lat(ncrn_dst,dst_size), 00212 & dst_corner_lon(ncrn_dst,dst_size) 00213 00214 !----------------------------------------------------------------------- 00215 ! 00216 ! local variables 00217 ! 00218 !----------------------------------------------------------------------- 00219 00220 integer (kind=int_kind) :: 00221 & n ! loop counter 00222 &, nele ! element loop counter 00223 &, i,j ! logical 2d addresses 00224 &, ip1,jp1 00225 &, n_add, e_add, ne_add 00226 &, nx, ny 00227 00228 real (kind=dbl_kind) :: 00229 & dlat, dlon ! lat/lon intervals for search bins 00230 00231 real (kind=dbl_kind), dimension(4) :: 00232 & tmp_lats, tmp_lons ! temps for computing bounding boxes 00233 ! 00234 !----------------------------------------------------------------------- 00235 ! 00236 IF (nlogprt .GE. 2) THEN 00237 WRITE (UNIT = nulou,FMT = *)' ' 00238 WRITE (UNIT = nulou,FMT = *)'Entering routine grid_init' 00239 WRITE (UNIT = nulou,FMT = *)' ' 00240 CALL FLUSH(nulou) 00241 ENDIF 00242 ! 00243 !----------------------------------------------------------------------- 00244 ! 00245 ! allocate grid coordinates/masks and read data 00246 ! 00247 !----------------------------------------------------------------------- 00248 00249 select case(m_method) 00250 case ('CONSERV') 00251 luse_grid_centers = .false. 00252 case ('BILINEAR') 00253 luse_grid_centers = .true. 00254 case ('BICUBIC') 00255 luse_grid_centers = .true. 00256 case ('DISTWGT') 00257 luse_grid_centers = .true. 00258 case ('GAUSWGT') 00259 luse_grid_centers = .true. 00260 case default 00261 stop 'unknown mapping method' 00262 end select 00263 00264 allocate( grid1_mask (src_size), 00265 & grid2_mask (dst_size), 00266 & grid1_center_lat(src_size), 00267 & grid1_center_lon(src_size), 00268 & grid2_center_lat(dst_size), 00269 & grid2_center_lon(dst_size), 00270 & grid1_area (src_size), 00271 & grid2_area (dst_size), 00272 & grid1_frac (src_size), 00273 & grid2_frac (dst_size), 00274 & grid1_dims (src_rank), 00275 & grid2_dims (dst_rank), 00276 & grid1_bound_box (4 , src_size), 00277 & grid2_bound_box (4 , dst_size)) 00278 00279 if (.not. luse_grid_centers) then 00280 allocate( grid1_corner_lat(ncrn_src, src_size), 00281 & grid1_corner_lon(ncrn_src, src_size), 00282 & grid2_corner_lat(ncrn_dst, dst_size), 00283 & grid2_corner_lon(ncrn_dst, dst_size)) 00284 endif 00285 00286 !----------------------------------------------------------------------- 00287 ! 00288 ! copy input data to module data 00289 ! 00290 !----------------------------------------------------------------------- 00291 00292 restrict_type = rst_type 00293 num_srch_bins = n_srch_bins 00294 grid1_size = src_size 00295 grid2_size = dst_size 00296 grid1_dims = src_dims 00297 grid2_dims = dst_dims 00298 grid1_rank = src_rank 00299 grid2_rank = dst_rank 00300 grid1_corners = ncrn_src 00301 grid2_corners = ncrn_dst 00302 grid1_name = src_name 00303 grid2_name = dst_name 00304 grid1_center_lat = src_lat 00305 grid1_center_lon = src_lon 00306 grid2_center_lat = dst_lat 00307 grid2_center_lon = dst_lon 00308 00309 if (.not. luse_grid_centers) then 00310 grid1_corner_lat = src_corner_lat 00311 grid1_corner_lon = src_corner_lon 00312 grid2_corner_lat = dst_corner_lat 00313 grid2_corner_lon = dst_corner_lon 00314 endif 00315 00316 c if (luse_grid1_area) then 00317 c grid1_area_in 00318 c endif 00319 c if (luse_grid2_area) then 00320 c grid2_area_in 00321 c endif 00322 00323 grid1_area = zero 00324 grid1_frac = zero 00325 grid2_area = zero 00326 grid2_frac = zero 00327 00328 00329 !----------------------------------------------------------------------- 00330 ! 00331 ! initialize logical mask and convert lat/lon units if required 00332 ! 00333 !----------------------------------------------------------------------- 00334 where (src_mask == 1) 00335 grid1_mask = .true. 00336 elsewhere 00337 grid1_mask = .false. 00338 endwhere 00339 00340 where (dst_mask == 1) 00341 grid2_mask = .true. 00342 elsewhere 00343 grid2_mask = .false. 00344 endwhere 00345 00346 C 00347 C* -- convert unit from degrees (OASIS unit) to radians 00348 C 00349 grid1_center_lat = grid1_center_lat*deg2rad 00350 grid1_center_lon = grid1_center_lon*deg2rad 00351 grid2_center_lat = grid2_center_lat*deg2rad 00352 grid2_center_lon = grid2_center_lon*deg2rad 00353 00354 if (.not. luse_grid_centers) then 00355 grid1_corner_lat = grid1_corner_lat*deg2rad 00356 grid1_corner_lon = grid1_corner_lon*deg2rad 00357 grid2_corner_lat = grid2_corner_lat*deg2rad 00358 grid2_corner_lon = grid2_corner_lon*deg2rad 00359 endif 00360 00361 grid1_units='radians' 00362 grid2_units='radians' 00363 !----------------------------------------------------------------------- 00364 ! 00365 ! convert longitudes to 0,2pi interval 00366 ! 00367 !----------------------------------------------------------------------- 00368 00369 where (grid1_center_lon .gt. pi2) grid1_center_lon = 00370 & grid1_center_lon - pi2 00371 where (grid1_center_lon .lt. zero) grid1_center_lon = 00372 & grid1_center_lon + pi2 00373 where (grid2_center_lon .gt. pi2) grid2_center_lon = 00374 & grid2_center_lon - pi2 00375 where (grid2_center_lon .lt. zero) grid2_center_lon = 00376 & grid2_center_lon + pi2 00377 00378 if (.not. luse_grid_centers) then 00379 where (grid1_corner_lon .gt. pi2) grid1_corner_lon = 00380 & grid1_corner_lon - pi2 00381 where (grid1_corner_lon .lt. zero) grid1_corner_lon = 00382 & grid1_corner_lon + pi2 00383 where (grid2_corner_lon .gt. pi2) grid2_corner_lon = 00384 & grid2_corner_lon - pi2 00385 where (grid2_corner_lon .lt. zero) grid2_corner_lon = 00386 & grid2_corner_lon + pi2 00387 endif 00388 !----------------------------------------------------------------------- 00389 ! 00390 ! make sure input latitude range is within the machine values 00391 ! for +/- pi/2 00392 ! 00393 !----------------------------------------------------------------------- 00394 00395 where (grid1_center_lat > pih) grid1_center_lat = pih 00396 where (grid1_center_lat < -pih) grid1_center_lat = -pih 00397 where (grid2_center_lat > pih) grid2_center_lat = pih 00398 where (grid2_center_lat < -pih) grid2_center_lat = -pih 00399 00400 if (.not. luse_grid_centers) then 00401 where (grid1_corner_lat > pih) grid1_corner_lat = pih 00402 where (grid1_corner_lat < -pih) grid1_corner_lat = -pih 00403 where (grid2_corner_lat > pih) grid2_corner_lat = pih 00404 where (grid2_corner_lat < -pih) grid2_corner_lat = -pih 00405 endif 00406 00407 !----------------------------------------------------------------------- 00408 ! 00409 ! compute bounding boxes for restricting future grid searches 00410 ! 00411 !----------------------------------------------------------------------- 00412 00413 if (.not. luse_grid_centers) then 00414 00415 grid1_bound_box(1,:) = minval(grid1_corner_lat, DIM=1) 00416 grid1_bound_box(2,:) = maxval(grid1_corner_lat, DIM=1) 00417 grid1_bound_box(3,:) = minval(grid1_corner_lon, DIM=1) 00418 grid1_bound_box(4,:) = maxval(grid1_corner_lon, DIM=1) 00419 00420 grid2_bound_box(1,:) = minval(grid2_corner_lat, DIM=1) 00421 grid2_bound_box(2,:) = maxval(grid2_corner_lat, DIM=1) 00422 grid2_bound_box(3,:) = minval(grid2_corner_lon, DIM=1) 00423 grid2_bound_box(4,:) = maxval(grid2_corner_lon, DIM=1) 00424 else 00425 00426 nx = grid1_dims(1) 00427 ny = grid1_dims(2) 00428 00429 do n=1,grid1_size 00430 00431 !*** find N,S and NE points to this grid point 00432 00433 j = (n - 1)/nx +1 00434 i = n - (j-1)*nx 00435 00436 if (i < nx) then 00437 ip1 = i + 1 00438 else 00439 !*** assume cyclic 00440 ip1 = 1 00441 !*** but if it is not, correct 00442 e_add = (j - 1)*nx + ip1 00443 if (abs(grid1_center_lat(e_add) - 00444 & grid1_center_lat(n )) > pih) then 00445 ip1 = i 00446 endif 00447 endif 00448 00449 if (j < ny) then 00450 jp1 = j+1 00451 else 00452 !*** assume cyclic 00453 jp1 = 1 00454 !*** but if it is not, correct 00455 n_add = (jp1 - 1)*nx + i 00456 if (abs(grid1_center_lat(n_add) - 00457 & grid1_center_lat(n )) > pih) then 00458 jp1 = j 00459 endif 00460 endif 00461 00462 n_add = (jp1 - 1)*nx + i 00463 e_add = (j - 1)*nx + ip1 00464 ne_add = (jp1 - 1)*nx + ip1 00465 00466 !*** find N,S and NE lat/lon coords and check bounding box 00467 00468 tmp_lats(1) = grid1_center_lat(n) 00469 tmp_lats(2) = grid1_center_lat(e_add) 00470 tmp_lats(3) = grid1_center_lat(ne_add) 00471 tmp_lats(4) = grid1_center_lat(n_add) 00472 00473 tmp_lons(1) = grid1_center_lon(n) 00474 tmp_lons(2) = grid1_center_lon(e_add) 00475 tmp_lons(3) = grid1_center_lon(ne_add) 00476 tmp_lons(4) = grid1_center_lon(n_add) 00477 00478 grid1_bound_box(1,n) = minval(tmp_lats) 00479 grid1_bound_box(2,n) = maxval(tmp_lats) 00480 grid1_bound_box(3,n) = minval(tmp_lons) 00481 grid1_bound_box(4,n) = maxval(tmp_lons) 00482 end do 00483 00484 nx = grid2_dims(1) 00485 ny = grid2_dims(2) 00486 00487 do n=1,grid2_size 00488 00489 !*** find N,S and NE points to this grid point 00490 00491 j = (n - 1)/nx +1 00492 i = n - (j-1)*nx 00493 00494 if (i < nx) then 00495 ip1 = i + 1 00496 else 00497 !*** assume cyclic 00498 ip1 = 1 00499 !*** but if it is not, correct 00500 e_add = (j - 1)*nx + ip1 00501 if (abs(grid2_center_lat(e_add) - 00502 & grid2_center_lat(n )) > pih) then 00503 ip1 = i 00504 endif 00505 endif 00506 00507 if (j < ny) then 00508 jp1 = j+1 00509 else 00510 !*** assume cyclic 00511 jp1 = 1 00512 !*** but if it is not, correct 00513 n_add = (jp1 - 1)*nx + i 00514 if (abs(grid2_center_lat(n_add) - 00515 & grid2_center_lat(n )) > pih) then 00516 jp1 = j 00517 endif 00518 endif 00519 00520 n_add = (jp1 - 1)*nx + i 00521 e_add = (j - 1)*nx + ip1 00522 ne_add = (jp1 - 1)*nx + ip1 00523 00524 !*** find N,S and NE lat/lon coords and check bounding box 00525 00526 tmp_lats(1) = grid2_center_lat(n) 00527 tmp_lats(2) = grid2_center_lat(e_add) 00528 tmp_lats(3) = grid2_center_lat(ne_add) 00529 tmp_lats(4) = grid2_center_lat(n_add) 00530 00531 tmp_lons(1) = grid2_center_lon(n) 00532 tmp_lons(2) = grid2_center_lon(e_add) 00533 tmp_lons(3) = grid2_center_lon(ne_add) 00534 tmp_lons(4) = grid2_center_lon(n_add) 00535 00536 grid2_bound_box(1,n) = minval(tmp_lats) 00537 grid2_bound_box(2,n) = maxval(tmp_lats) 00538 grid2_bound_box(3,n) = minval(tmp_lons) 00539 grid2_bound_box(4,n) = maxval(tmp_lons) 00540 end do 00541 00542 endif 00543 00544 where (abs(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi) 00545 grid1_bound_box(3,:) = zero 00546 grid1_bound_box(4,:) = pi2 00547 end where 00548 00549 where (abs(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi) 00550 grid2_bound_box(3,:) = zero 00551 grid2_bound_box(4,:) = pi2 00552 end where 00553 00554 !*** 00555 !*** try to check for cells that overlap poles 00556 !*** 00557 00558 where (grid1_center_lat > grid1_bound_box(2,:)) 00559 & grid1_bound_box(2,:) = pih 00560 00561 where (grid1_center_lat < grid1_bound_box(1,:)) 00562 & grid1_bound_box(1,:) = -pih 00563 00564 where (grid2_center_lat > grid2_bound_box(2,:)) 00565 & grid2_bound_box(2,:) = pih 00566 00567 where (grid2_center_lat < grid2_bound_box(1,:)) 00568 & grid2_bound_box(1,:) = -pih 00569 00570 !----------------------------------------------------------------------- 00571 ! 00572 ! set up and assign address ranges to search bins in order to 00573 ! further restrict later searches 00574 ! 00575 !----------------------------------------------------------------------- 00576 00577 select case (restrict_type) 00578 00579 case ('LATITUDE') 00580 write(stdout,*) 'Using latitude bins to restrict search.' 00581 00582 allocate(bin_addr1(2,num_srch_bins)) 00583 allocate(bin_addr2(2,num_srch_bins)) 00584 allocate(bin_lats (2,num_srch_bins)) 00585 allocate(bin_lons (2,num_srch_bins)) 00586 00587 dlat = pi/num_srch_bins 00588 00589 do n=1,num_srch_bins 00590 bin_lats(1,n) = (n-1)*dlat - pih 00591 bin_lats(2,n) = n*dlat - pih 00592 bin_lons(1,n) = zero 00593 bin_lons(2,n) = pi2 00594 bin_addr1(1,n) = grid1_size + 1 00595 bin_addr1(2,n) = 0 00596 bin_addr2(1,n) = grid2_size + 1 00597 bin_addr2(2,n) = 0 00598 end do 00599 00600 do nele=1,grid1_size 00601 do n=1,num_srch_bins 00602 if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and. 00603 & grid1_bound_box(2,nele) >= bin_lats(1,n)) then 00604 bin_addr1(1,n) = min(nele,bin_addr1(1,n)) 00605 bin_addr1(2,n) = max(nele,bin_addr1(2,n)) 00606 endif 00607 end do 00608 end do 00609 00610 do nele=1,grid2_size 00611 do n=1,num_srch_bins 00612 if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and. 00613 & grid2_bound_box(2,nele) >= bin_lats(1,n)) then 00614 bin_addr2(1,n) = min(nele,bin_addr2(1,n)) 00615 bin_addr2(2,n) = max(nele,bin_addr2(2,n)) 00616 endif 00617 end do 00618 end do 00619 00620 case ('LATLON') 00621 write(stdout,*) 'Using lat/lon boxes to restrict search.' 00622 00623 dlat = pi /num_srch_bins 00624 dlon = pi2/num_srch_bins 00625 00626 allocate(bin_addr1(2,num_srch_bins*num_srch_bins)) 00627 allocate(bin_addr2(2,num_srch_bins*num_srch_bins)) 00628 allocate(bin_lats (2,num_srch_bins*num_srch_bins)) 00629 allocate(bin_lons (2,num_srch_bins*num_srch_bins)) 00630 00631 n = 0 00632 do j=1,num_srch_bins 00633 do i=1,num_srch_bins 00634 n = n + 1 00635 00636 bin_lats(1,n) = (j-1)*dlat - pih 00637 bin_lats(2,n) = j*dlat - pih 00638 bin_lons(1,n) = (i-1)*dlon 00639 bin_lons(2,n) = i*dlon 00640 bin_addr1(1,n) = grid1_size + 1 00641 bin_addr1(2,n) = 0 00642 bin_addr2(1,n) = grid2_size + 1 00643 bin_addr2(2,n) = 0 00644 end do 00645 end do 00646 00647 num_srch_bins = num_srch_bins**2 00648 00649 do nele=1,grid1_size 00650 do n=1,num_srch_bins 00651 if (grid1_bound_box(1,nele) <= bin_lats(2,n) .and. 00652 & grid1_bound_box(2,nele) >= bin_lats(1,n) .and. 00653 & grid1_bound_box(3,nele) <= bin_lons(2,n) .and. 00654 & grid1_bound_box(4,nele) >= bin_lons(1,n)) then 00655 bin_addr1(1,n) = min(nele,bin_addr1(1,n)) 00656 bin_addr1(2,n) = max(nele,bin_addr1(2,n)) 00657 endif 00658 end do 00659 end do 00660 00661 do nele=1,grid2_size 00662 do n=1,num_srch_bins 00663 if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and. 00664 & grid2_bound_box(2,nele) >= bin_lats(1,n) .and. 00665 & grid2_bound_box(3,nele) <= bin_lons(2,n) .and. 00666 & grid2_bound_box(4,nele) >= bin_lons(1,n)) then 00667 bin_addr2(1,n) = min(nele,bin_addr2(1,n)) 00668 bin_addr2(2,n) = max(nele,bin_addr2(2,n)) 00669 endif 00670 end do 00671 end do 00672 00673 case ('REDUCED') 00674 write(stdout,*) 00675 & '| Using reduced bins to restrict search. Reduced grids with' 00676 write(stdout,*) 00677 & '| a maximum of 500*NBRBINS latitude circles can be treated' 00678 00679 allocate(bin_addr2(2,num_srch_bins)) 00680 allocate(bin_lats (2,num_srch_bins)) 00681 allocate(bin_lons (2,num_srch_bins)) 00682 allocate(bin_addr1_r(4,500*num_srch_bins)) 00683 allocate(bin_lats_r (2,500*num_srch_bins)) 00684 allocate(bin_lons_r (2,500*num_srch_bins)) 00685 00686 dlat = pi/num_srch_bins 00687 00688 do n=1,num_srch_bins 00689 bin_lats(1,n) = (n-1)*dlat - pih 00690 bin_lats(2,n) = n*dlat - pih 00691 bin_lons(1,n) = zero 00692 bin_lons(2,n) = pi2 00693 bin_addr2(1,n) = grid2_size + 1 00694 bin_addr2(2,n) = 0 00695 end do 00696 00697 do nele=1,grid2_size 00698 do n=1,num_srch_bins 00699 if (grid2_bound_box(1,nele) <= bin_lats(2,n) .and. 00700 & grid2_bound_box(2,nele) >= bin_lats(1,n)) then 00701 bin_addr2(1,n) = min(nele,bin_addr2(1,n)) 00702 bin_addr2(2,n) = max(nele,bin_addr2(2,n)) 00703 endif 00704 end do 00705 end DO 00706 00707 00708 bin_addr1_r(1,1) = 0 00709 bin_lats_r(1,1) = grid1_center_lat(1) 00710 num_srch_red = 1 00711 do nele=1,grid1_size-1 00712 if (grid1_center_lat(nele+1) == 00713 & grid1_center_lat(nele)) THEN 00714 bin_addr1_r(2,num_srch_red) = nele+1 00715 else IF (grid1_center_lat(nele+1) /= 00716 & grid1_center_lat(nele)) THEN 00717 bin_addr1_r(1,num_srch_red+1) = nele+1 00718 bin_lats_r(2,num_srch_red) = grid1_center_lat(nele+1) 00719 bin_lats_r(1,num_srch_red+1) = 00720 & bin_lats_r(2,num_srch_red) 00721 num_srch_red = num_srch_red+1 00722 ENDIF 00723 end DO 00724 00725 DO nele = 1,num_srch_red-1 00726 bin_addr1_r(3,nele)=bin_addr1_r(1,nele+1) 00727 bin_addr1_r(4,nele)=bin_addr1_r(2,nele+1) 00728 enddo 00729 bin_addr1_r(3,num_srch_red) = bin_addr1_r(1,num_srch_red-1) 00730 bin_addr1_r(4,num_srch_red) = bin_addr1_r(2,num_srch_red-1) 00731 00732 case default 00733 stop 'unknown search restriction method' 00734 end select 00735 ! 00736 IF (nlogprt .GE. 2) THEN 00737 WRITE (UNIT = nulou,FMT = *)' ' 00738 WRITE (UNIT = nulou,FMT = *)'Leaving routine grid_init' 00739 WRITE (UNIT = nulou,FMT = *)' ' 00740 CALL FLUSH(nulou) 00741 ENDIF 00742 ! 00743 !----------------------------------------------------------------------- 00744 00745 end subroutine grid_init 00746 00747 !*********************************************************************** 00748 00749 subroutine free_grids 00750 00751 !----------------------------------------------------------------------- 00752 ! 00753 ! subroutine to deallocate allocated arrays 00754 ! 00755 !----------------------------------------------------------------------- 00756 ! 00757 IF (nlogprt .GE. 2) THEN 00758 WRITE (UNIT = nulou,FMT = *)' ' 00759 WRITE (UNIT = nulou,FMT = *)'Entering routine free_grid' 00760 WRITE (UNIT = nulou,FMT = *)' ' 00761 CALL FLUSH(nulou) 00762 ENDIF 00763 ! 00764 deallocate(grid1_mask, grid2_mask, 00765 & grid1_center_lat, grid1_center_lon, 00766 & grid2_center_lat, grid2_center_lon, 00767 & grid1_area, grid2_area, 00768 & grid1_frac, grid2_frac, 00769 & grid1_dims, grid2_dims) 00770 00771 IF (restrict_TYPE == 'REDUCED') then 00772 deallocate( grid1_bound_box, grid2_bound_box, 00773 & bin_addr1_r, bin_addr2, 00774 & bin_lons, bin_lats, 00775 & bin_lats_r, bin_lons_r) 00776 else 00777 deallocate( grid1_bound_box, grid2_bound_box, 00778 & bin_addr1, bin_addr2, 00779 & bin_lats, bin_lons) 00780 endif 00781 00782 if (.not. luse_grid_centers) then 00783 deallocate(grid1_corner_lat, grid1_corner_lon, 00784 & grid2_corner_lat, grid2_corner_lon) 00785 ENDIF 00786 ! 00787 IF (nlogprt .GE. 2) THEN 00788 WRITE (UNIT = nulou,FMT = *)' ' 00789 WRITE (UNIT = nulou,FMT = *)'Leaving routine free_grid' 00790 WRITE (UNIT = nulou,FMT = *)' ' 00791 CALL FLUSH(nulou) 00792 ENDIF 00793 !----------------------------------------------------------------------- 00794 00795 end subroutine free_grids 00796 00797 !*********************************************************************** 00798 00799 end module grids 00800 00801 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00802
1.7.4 
