|
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 a modified version 00009 C of SCRIP 1.4 remap_ditswgt.f in order to weight a neighbour by 00010 C exp[-1/2 * d^^2/sigma^^2] where d is the distance between the source 00011 C and the target points and sigma^^2=VAR*dm^^2, where dm is the 00012 C average distance between the source grid points and VAR is a value 00013 C given by the user. 00014 C 00015 C Additional Modifications: 00016 C - restrict types are written in capital letters 00017 C - bug line 429: bin_lons(3,n) instead of bin_lons(2,n) 00018 C 00019 C Modified by D. Declat, CERFACS 27.06.2002 00020 C*********************************************************************** 00021 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00022 ! 00023 ! this module contains necessary routines for performing an 00024 ! interpolation using a distance-weighted average of n nearest 00025 ! neighbors. 00026 ! 00027 !----------------------------------------------------------------------- 00028 ! 00029 ! CVS:$Id: remap_gauswgt.f 2826 2010-12-10 11:14:21Z valcke $ 00030 ! 00031 ! Copyright (c) 1997, 1998 the Regents of the University of 00032 ! California. 00033 ! 00034 ! This software and ancillary information (herein called software) 00035 ! called SCRIP is made available under the terms described here. 00036 ! The software has been approved for release with associated 00037 ! LA-CC Number 98-45. 00038 ! 00039 ! Unless otherwise indicated, this software has been authored 00040 ! by an employee or employees of the University of California, 00041 ! operator of the Los Alamos National Laboratory under Contract 00042 ! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S. 00043 ! Government has rights to use, reproduce, and distribute this 00044 ! software. The public may copy and use this software without 00045 ! charge, provided that this Notice and any statement of authorship 00046 ! are reproduced on all copies. Neither the Government nor the 00047 ! University makes any warranty, express or implied, or assumes 00048 ! any liability or responsibility for the use of this software. 00049 ! 00050 ! If software is modified to produce derivative works, such modified 00051 ! software should be clearly marked, so as not to confuse it with 00052 ! the version available from Los Alamos National Laboratory. 00053 ! 00054 !*********************************************************************** 00055 00056 module remap_gaussian_weight 00057 00058 !----------------------------------------------------------------------- 00059 00060 use kinds_mod ! defines common data types 00061 use constants ! defines common constants 00062 use grids ! module containing grid info 00063 use remap_vars ! module containing remap info 00064 00065 implicit none 00066 00067 !----------------------------------------------------------------------- 00068 ! 00069 ! module variables 00070 ! 00071 !----------------------------------------------------------------------- 00072 00073 real (kind=dbl_kind), dimension(:), allocatable, save :: 00074 & coslat, sinlat, ! cosine, sine of grid lats (for distance) 00075 & coslon, sinlon, ! cosine, sine of grid lons (for distance) 00076 & wgtstmp ! an array to hold the link weight 00077 00078 !*********************************************************************** 00079 00080 contains 00081 00082 !*********************************************************************** 00083 00084 subroutine remap_gauswgt (lextrapdone, num_neighbors, rl_varmul) 00085 00086 !----------------------------------------------------------------------- 00087 ! 00088 ! this routine computes the inverse-distance weights for a 00089 ! nearest-neighbor interpolation. 00090 ! 00091 !----------------------------------------------------------------------- 00092 00093 LOGICAL :: 00094 & lextrapdone ! logical, true if EXTRAP done on field 00095 00096 REAL (kind=real_kind) :: 00097 $ rl_varmul ! Gaussian variance 00098 00099 INTEGER (kind=int_kind) :: 00100 & num_neighbors ! number of neighbours 00101 00102 !----------------------------------------------------------------------- 00103 ! 00104 ! local variables 00105 ! 00106 !----------------------------------------------------------------------- 00107 00108 logical (kind=log_kind), dimension(num_neighbors) :: 00109 & nbr_mask ! mask at nearest neighbors 00110 00111 integer (kind=int_kind) :: n, 00112 & dst_add, ! destination address 00113 & nmap ! index of current map being computed 00114 00115 integer (kind=int_kind), dimension(num_neighbors) :: 00116 & nbr_add ! source address at nearest neighbors 00117 00118 real (kind=dbl_kind), dimension(num_neighbors) :: 00119 & nbr_dist ! angular distance four nearest neighbors 00120 00121 real (kind=dbl_kind) :: 00122 & coslat_dst, ! cos(lat) of destination grid point 00123 & coslon_dst, ! cos(lon) of destination grid point 00124 & sinlat_dst, ! sin(lat) of destination grid point 00125 & sinlon_dst, ! sin(lon) of destination grid point 00126 & dist_tot, ! sum of neighbor distances (for normalizing) 00127 & dist_average ! average distance between the source points 00128 logical (kind=log_kind) :: ll_allmask, ll_nnei 00129 real (kind=dbl_kind) :: 00130 & distance ,plat,plon,src_latsnn, arg ! angular distance 00131 real (kind=dbl_kind), dimension (1) :: wgts_new 00132 integer (kind=int_kind) :: min_add_out, max_add_out, 00133 & srch_add, src_addnn 00134 !----------------------------------------------------------------------- 00135 ! 00136 IF (nlogprt .GE. 2) THEN 00137 WRITE (UNIT = nulou,FMT = *)' ' 00138 WRITE (UNIT = nulou,FMT = *)'Entering routine remap_gauswgt' 00139 CALL FLUSH(nulou) 00140 ENDIF 00141 !----------------------------------------------------------------------- 00142 ! 00143 ! compute mappings from grid1 to grid2 00144 ! 00145 !----------------------------------------------------------------------- 00146 00147 nmap = 1 00148 00149 !*** 00150 !*** allocate wgtstmp to be consistent with store_link interface 00151 !*** 00152 00153 allocate (wgtstmp(num_wts)) 00154 00155 !*** 00156 !*** compute cos, sin of lat/lon on source grid for distance 00157 !*** calculations 00158 !*** 00159 00160 allocate (coslat(grid1_size), coslon(grid1_size), 00161 & sinlat(grid1_size), sinlon(grid1_size)) 00162 00163 coslat = cos(grid1_center_lat) 00164 coslon = cos(grid1_center_lon) 00165 sinlat = sin(grid1_center_lat) 00166 sinlon = sin(grid1_center_lon) 00167 00168 !*** 00169 !*** compute the average of the distances between the source points 00170 !*** 00171 00172 call grid_dist_average(grid1_size, 00173 & coslat, coslon, 00174 & sinlat, sinlon, 00175 & dist_average) 00176 !*** 00177 !*** loop over destination grid 00178 !*** 00179 00180 grid_loop1: do dst_add = 1, grid2_size 00181 00182 if (.not. grid2_mask(dst_add)) cycle grid_loop1 00183 00184 coslat_dst = cos(grid2_center_lat(dst_add)) 00185 coslon_dst = cos(grid2_center_lon(dst_add)) 00186 sinlat_dst = sin(grid2_center_lat(dst_add)) 00187 sinlon_dst = sin(grid2_center_lon(dst_add)) 00188 00189 !*** 00190 !*** find nearest grid points on source grid and 00191 !*** distances to each point 00192 !*** 00193 00194 call grid_search_nbrg(nbr_add, nbr_dist, 00195 & min_add_out, max_add_out, 00196 & grid2_center_lat(dst_add), 00197 & grid2_center_lon(dst_add), 00198 & coslat_dst, coslon_dst, 00199 & sinlat_dst, sinlon_dst, 00200 & bin_addr1, 00201 & dist_average, num_neighbors, rl_varmul) 00202 00203 !*** 00204 !*** compute weights based on inverse distance 00205 !*** if mask is false, eliminate those points 00206 !*** 00207 00208 dist_tot = zero 00209 do n=1,num_neighbors 00210 if ((grid1_mask(nbr_add(n))) .or. 00211 & (.not. grid1_mask(nbr_add(n)) .and. lextrapdone)) THEN 00212 nbr_dist(n) = one/nbr_dist(n) 00213 dist_tot = dist_tot + nbr_dist(n) 00214 nbr_mask(n) = .true. 00215 else 00216 nbr_mask(n) = .false. 00217 endif 00218 end do 00219 00220 !*** 00221 !*** normalize weights and store the link 00222 !*** 00223 00224 do n=1,num_neighbors 00225 if (nbr_mask(n)) then 00226 wgtstmp(1) = nbr_dist(n)/dist_tot 00227 call store_link_nbr(nbr_add(n), dst_add, wgtstmp, nmap) 00228 grid2_frac(dst_add) = one 00229 endif 00230 end do 00231 ll_nnei = .false. 00232 IF (ll_nnei) THEN 00233 ll_allmask= .true. 00234 do n=1,num_neighbors 00235 if (nbr_mask(n)) then 00236 ll_allmask=.false. 00237 endif 00238 end do 00239 if ( ll_allmask) THEN 00240 IF (nlogprt .GE. 2) THEN 00241 WRITE(nulou,*)'ll_allmask true',src_addnn 00242 CALL FLUSH(nulou) 00243 ENDIF 00244 src_latsnn = bignum 00245 do srch_add = min_add_out,max_add_out 00246 if (grid1_mask(srch_add) .or. 00247 & (.not. grid1_mask(srch_add) 00248 & .and. lextrapdone)) THEN 00249 arg = coslat_dst*cos(grid1_center_lat(srch_add))* 00250 & (coslon_dst*cos(grid1_center_lon(srch_add)) + 00251 & sinlon_dst*sin(grid1_center_lon(srch_add)))+ 00252 & sinlat_dst*sin(grid1_center_lat(srch_add)) 00253 IF (arg < -1.0d0) THEN 00254 arg = -1.0d0 00255 ELSE IF (arg > 1.0d0) THEN 00256 arg = 1.0d0 00257 END IF 00258 distance=acos(arg) 00259 00260 if (distance < src_latsnn) then 00261 src_addnn = srch_add 00262 src_latsnn = distance 00263 endif 00264 endif 00265 end do 00266 wgts_new(1) = 1. 00267 grid2_frac(dst_add) = one 00268 call store_link_nbr(src_addnn,dst_add ,wgts_new, nmap) 00269 endif 00270 endif 00271 00272 end do grid_loop1 00273 00274 deallocate (coslat, coslon, sinlat, sinlon) 00275 00276 !----------------------------------------------------------------------- 00277 ! 00278 ! compute mappings from grid2 to grid1 if necessary 00279 ! 00280 !----------------------------------------------------------------------- 00281 00282 if (num_maps > 1) then 00283 00284 nmap = 2 00285 00286 !*** 00287 !*** compute cos, sin of lat/lon on source grid for distance 00288 !*** calculations 00289 !*** 00290 00291 allocate (coslat(grid2_size), coslon(grid2_size), 00292 & sinlat(grid2_size), sinlon(grid2_size)) 00293 00294 coslat = cos(grid2_center_lat) 00295 coslon = cos(grid2_center_lon) 00296 sinlat = sin(grid2_center_lat) 00297 sinlon = sin(grid2_center_lon) 00298 00299 !*** 00300 !*** compute the average of the distances between the source points 00301 !*** 00302 00303 call grid_dist_average(grid2_size, 00304 & coslat, coslon, 00305 & sinlat, sinlon, 00306 & dist_average) 00307 00308 !*** 00309 !*** loop over destination grid 00310 !*** 00311 00312 grid_loop2: do dst_add = 1, grid1_size 00313 00314 if (.not. grid1_mask(dst_add)) cycle grid_loop2 00315 00316 coslat_dst = cos(grid1_center_lat(dst_add)) 00317 coslon_dst = cos(grid1_center_lon(dst_add)) 00318 sinlat_dst = sin(grid1_center_lat(dst_add)) 00319 sinlon_dst = sin(grid1_center_lon(dst_add)) 00320 00321 !*** 00322 !*** find four nearest grid points on source grid and 00323 !*** distances to each point 00324 !*** 00325 00326 call grid_search_nbrg(nbr_add, nbr_dist, 00327 & min_add_out, max_add_out, 00328 & grid1_center_lat(dst_add), 00329 & grid1_center_lon(dst_add), 00330 & coslat_dst, coslon_dst, 00331 & sinlat_dst, sinlon_dst, 00332 & bin_addr2, 00333 & dist_average, num_neighbors, rl_varmul) 00334 00335 !*** 00336 !*** compute weights based on inverse distance 00337 !*** if mask is false, eliminate those points 00338 !*** 00339 00340 dist_tot = zero 00341 do n=1,num_neighbors 00342 if (grid2_mask(nbr_add(n))) then 00343 nbr_dist(n) = one/nbr_dist(n) 00344 dist_tot = dist_tot + nbr_dist(n) 00345 nbr_mask(n) = .true. 00346 else 00347 nbr_mask(n) = .false. 00348 endif 00349 end do 00350 00351 !*** 00352 !*** normalize weights and store the link 00353 !*** 00354 00355 do n=1,num_neighbors 00356 if (nbr_mask(n)) then 00357 wgtstmp(1) = nbr_dist(n)/dist_tot 00358 call store_link_nbr(dst_add, nbr_add(n), wgtstmp, nmap) 00359 grid1_frac(dst_add) = one 00360 endif 00361 end do 00362 00363 IF (ll_nnei) THEN 00364 ll_allmask= .true. 00365 do n=1,num_neighbors 00366 if (nbr_mask(n)) then 00367 ll_allmask=.false. 00368 endif 00369 end do 00370 if ( ll_allmask) then 00371 PRINT*, 'll_allmask true',src_addnn 00372 src_latsnn = bignum 00373 do srch_add = min_add_out,max_add_out 00374 if (grid2_mask(srch_add) .or. 00375 & (.not. grid2_mask(srch_add) 00376 & .and. lextrapdone)) THEN 00377 arg = coslat_dst*cos(grid2_center_lat(srch_add))* 00378 & (coslon_dst*cos(grid2_center_lon(srch_add)) + 00379 & sinlon_dst*sin(grid2_center_lon(srch_add)))+ 00380 & sinlat_dst*sin(grid2_center_lat(srch_add)) 00381 IF (arg < -1.0d0) THEN 00382 arg = -1.0d0 00383 ELSE IF (arg > 1.0d0) THEN 00384 arg = 1.0d0 00385 END IF 00386 distance=acos(arg) 00387 00388 if (distance < src_latsnn) then 00389 src_addnn = srch_add 00390 src_latsnn = distance 00391 endif 00392 endif 00393 end do 00394 wgts_new = 1. 00395 grid1_frac(dst_add) = one 00396 call store_link_nbr(dst_add, src_addnn, wgts_new, nmap) 00397 endif 00398 endif 00399 end do grid_loop2 00400 00401 deallocate (coslat, coslon, sinlat, sinlon) 00402 00403 endif 00404 00405 deallocate(wgtstmp) 00406 ! 00407 IF (nlogprt .GE. 2) THEN 00408 WRITE (UNIT = nulou,FMT = *)' ' 00409 WRITE (UNIT = nulou,FMT = *)'Leaving routine remap_gauswgt' 00410 CALL FLUSH(nulou) 00411 ENDIF 00412 ! 00413 !----------------------------------------------------------------------- 00414 00415 end subroutine remap_gauswgt 00416 00417 !*********************************************************************** 00418 00419 subroutine grid_search_nbrg(nbr_add, nbr_dist, min_add, max_add, 00420 & plat, plon, 00421 & coslat_dst, coslon_dst, sinlat_dst, sinlon_dst, 00422 & src_bin_add, dst_average, 00423 $ num_neighbors, rl_varmul) 00424 00425 !----------------------------------------------------------------------- 00426 ! 00427 ! this routine finds the closest num_neighbor points to a search 00428 ! point and computes a distance to each of the neighbors. 00429 ! 00430 !----------------------------------------------------------------------- 00431 00432 !----------------------------------------------------------------------- 00433 ! 00434 ! input variables 00435 ! 00436 !----------------------------------------------------------------------- 00437 00438 integer (kind=int_kind), dimension(:,:), intent(in) :: 00439 & src_bin_add ! search bins for restricting search 00440 00441 real (kind=dbl_kind), intent(in) :: 00442 & plat, ! latitude of the search point 00443 & plon, ! longitude of the search point 00444 & coslat_dst, ! cos(lat) of the search point 00445 & coslon_dst, ! cos(lon) of the search point 00446 & sinlat_dst, ! sin(lat) of the search point 00447 & sinlon_dst, ! sin(lon) of the search point 00448 & dst_average ! average distance between the source points 00449 00450 REAL (kind=real_kind), intent(in) :: 00451 & rl_varmul ! Gaussian variance 00452 00453 INTEGER (kind=int_kind) :: 00454 & num_neighbors ! number of neighbours 00455 00456 !----------------------------------------------------------------------- 00457 ! 00458 ! output variables 00459 ! 00460 !----------------------------------------------------------------------- 00461 00462 integer (kind=int_kind), dimension(num_neighbors), intent(out) :: 00463 & nbr_add ! address of each of the closest points 00464 00465 real (kind=dbl_kind), dimension(num_neighbors), intent(out) :: 00466 & nbr_dist ! distance to each of the closest points 00467 00468 integer (kind=int_kind),intent(out) :: min_add, max_add 00469 00470 !----------------------------------------------------------------------- 00471 ! 00472 ! local variables 00473 ! 00474 !----------------------------------------------------------------------- 00475 00476 integer (kind=int_kind) :: n, nmax, nadd, nchk, ! dummy indices 00477 & nm1, np1, i, j, ip1, im1, jp1, jm1 00478 00479 real (kind=dbl_kind) :: 00480 & distance, arg ! angular distance 00481 00482 real (kind=dbl_kind) :: 00483 & variance ! variance for the gaussian FUNCTION 00484 00485 00486 !----------------------------------------------------------------------- 00487 ! 00488 ! loop over source grid and find nearest neighbors 00489 ! 00490 !----------------------------------------------------------------------- 00491 00492 !*** 00493 !*** restrict the search using search bins 00494 !*** expand the bins to catch neighbors 00495 !*** 00496 00497 select case (restrict_type) 00498 case('LATITUDE') 00499 00500 do n=1,num_srch_bins 00501 if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n)) then 00502 min_add = src_bin_add(1,n) 00503 max_add = src_bin_add(2,n) 00504 00505 nm1 = max(n-1,1) 00506 np1 = min(n+1,num_srch_bins) 00507 00508 min_add = min(min_add,src_bin_add(1,nm1)) 00509 max_add = max(max_add,src_bin_add(2,nm1)) 00510 min_add = min(min_add,src_bin_add(1,np1)) 00511 max_add = max(max_add,src_bin_add(2,np1)) 00512 endif 00513 end do 00514 00515 case('LATLON') 00516 00517 n = 0 00518 nmax = nint(sqrt(real(num_srch_bins))) 00519 do j=1,nmax 00520 jp1 = min(j+1,nmax) 00521 jm1 = max(j-1,1) 00522 do i=1,nmax 00523 ip1 = min(i+1,nmax) 00524 im1 = max(i-1,1) 00525 00526 n = n+1 00527 if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. 00528 & plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then 00529 min_add = src_bin_add(1,n) 00530 max_add = src_bin_add(2,n) 00531 00532 nm1 = (jm1-1)*nmax + im1 00533 np1 = (jp1-1)*nmax + ip1 00534 nm1 = max(nm1,1) 00535 np1 = min(np1,num_srch_bins) 00536 00537 min_add = min(min_add,src_bin_add(1,nm1)) 00538 max_add = max(max_add,src_bin_add(2,nm1)) 00539 min_add = min(min_add,src_bin_add(1,np1)) 00540 max_add = max(max_add,src_bin_add(2,np1)) 00541 endif 00542 end do 00543 end do 00544 00545 end select 00546 00547 !*** 00548 !*** initialize distance and address arrays 00549 !*** 00550 00551 nbr_add = 0 00552 nbr_dist = bignum 00553 variance = rl_varmul*dst_average*dst_average 00554 do nadd=min_add,max_add 00555 00556 !*** 00557 !*** find distance to this point 00558 !*** 00559 arg = sinlat_dst*sinlat(nadd) + 00560 & coslat_dst*coslat(nadd)* 00561 & (coslon_dst*coslon(nadd) + 00562 & sinlon_dst*sinlon(nadd)) 00563 IF (arg < -1.0d0) THEN 00564 arg = -1.0d0 00565 ELSE IF (arg > 1.0d0) THEN 00566 arg = 1.0d0 00567 END IF 00568 distance = acos(arg) 00569 !distance = min(500., distance) !ts-sv 00570 !distance = exp(.5*distance*distance/variance) 00571 !*** 00572 !*** store the address and distance if this is one of the 00573 !*** smallest four so far 00574 !*** 00575 00576 check_loop: do nchk=1,num_neighbors 00577 if (distance .lt. nbr_dist(nchk)) THEN 00578 do n=num_neighbors,nchk+1,-1 00579 nbr_add(n) = nbr_add(n-1) 00580 nbr_dist(n) = nbr_dist(n-1) 00581 end do 00582 nbr_add(nchk) = nadd 00583 nbr_dist(nchk) = distance 00584 exit check_loop 00585 endif 00586 end do check_loop 00587 00588 end do 00589 00590 exp_loop: do nchk=1,num_neighbors 00591 nbr_dist(nchk) = 00592 & exp(.5*nbr_dist(nchk)*nbr_dist(nchk)/variance) 00593 end do exp_loop 00594 00595 !----------------------------------------------------------------------- 00596 00597 end subroutine grid_search_nbrg 00598 00599 !*********************************************************************** 00600 00601 subroutine store_link_nbr(add1, add2, weights, nmap) 00602 00603 !----------------------------------------------------------------------- 00604 ! 00605 ! this routine stores the address and weight for this link in 00606 ! the appropriate address and weight arrays and resizes those 00607 ! arrays if necessary. 00608 ! 00609 !----------------------------------------------------------------------- 00610 00611 !----------------------------------------------------------------------- 00612 ! 00613 ! input variables 00614 ! 00615 !----------------------------------------------------------------------- 00616 00617 integer (kind=int_kind), intent(in) :: 00618 & add1, ! address on grid1 00619 & add2, ! address on grid2 00620 & nmap ! identifies which direction for mapping 00621 00622 real (kind=dbl_kind), dimension(:), intent(in) :: 00623 & weights ! array of remapping weights for this link 00624 00625 !----------------------------------------------------------------------- 00626 ! 00627 ! increment number of links and check to see if remap arrays need 00628 ! to be increased to accomodate the new link. then store the 00629 ! link. 00630 ! 00631 !----------------------------------------------------------------------- 00632 00633 select case (nmap) 00634 case(1) 00635 00636 num_links_map1 = num_links_map1 + 1 00637 00638 if (num_links_map1 > max_links_map1) 00639 & call resize_remap_vars(1,resize_increment) 00640 00641 grid1_add_map1(num_links_map1) = add1 00642 grid2_add_map1(num_links_map1) = add2 00643 wts_map1 (:,num_links_map1) = weights 00644 00645 case(2) 00646 00647 num_links_map2 = num_links_map2 + 1 00648 00649 if (num_links_map2 > max_links_map2) 00650 & call resize_remap_vars(2,resize_increment) 00651 00652 grid1_add_map2(num_links_map2) = add1 00653 grid2_add_map2(num_links_map2) = add2 00654 wts_map2 (:,num_links_map2) = weights 00655 00656 end select 00657 00658 !----------------------------------------------------------------------- 00659 00660 end subroutine store_link_nbr 00661 00662 !*********************************************************************** 00663 00664 subroutine grid_dist_average(grid_size, 00665 & coslat_grid, coslon_grid, 00666 & sinlat_grid, sinlon_grid, 00667 & dist_average) 00668 00669 !----------------------------------------------------------------------- 00670 ! 00671 ! this routine computes the average distance between the points of a 00672 ! grid. 00673 ! 00674 !----------------------------------------------------------------------- 00675 00676 !----------------------------------------------------------------------- 00677 ! 00678 ! output variables 00679 ! 00680 !----------------------------------------------------------------------- 00681 00682 real (kind=dbl_kind), intent(out) :: 00683 & dist_average ! distance to each of the closest points 00684 00685 !----------------------------------------------------------------------- 00686 ! 00687 ! input variables 00688 ! 00689 !----------------------------------------------------------------------- 00690 00691 integer (kind=int_kind), intent(in) :: 00692 & grid_size 00693 00694 real (kind=dbl_kind), dimension(:), intent(in) :: 00695 & coslat_grid, ! cos(lat) of the grid points 00696 & coslon_grid, ! cos(lon) of the grid points 00697 & sinlat_grid, ! sin(lat) of the grid points 00698 & sinlon_grid ! sin(lon) of the grid points 00699 REAL (kind=dbl_kind) :: arg 00700 00701 !----------------------------------------------------------------------- 00702 ! 00703 ! local variables 00704 ! 00705 !----------------------------------------------------------------------- 00706 00707 integer (kind=int_kind) :: i 00708 ! 00709 IF (nlogprt .GE. 2) THEN 00710 WRITE (UNIT = nulou,FMT = *)' ' 00711 WRITE (UNIT = nulou,FMT = *) 00712 & 'Entering routine remap_dist_average' 00713 CALL FLUSH(nulou) 00714 ENDIF 00715 ! 00716 !----------------------------------------------------------------------- 00717 ! 00718 ! compute the distance over the grid and average 00719 ! 00720 !----------------------------------------------------------------------- 00721 00722 dist_average = 0.0 00723 DO i = 1, grid_size 00724 IF (i .eq. 1) THEN 00725 arg = sinlat_grid(grid_size)*sinlat_grid(i) + 00726 & coslat_grid(grid_size)*coslat_grid(i)* 00727 & (coslon_grid(grid_size)*coslon_grid(i) + 00728 & sinlon_grid(grid_size)*sinlon_grid(i)) 00729 IF (arg < -1.0d0) THEN 00730 arg = -1.0d0 00731 ELSE IF (arg > 1.0d0) THEN 00732 arg = 1.0d0 00733 END IF 00734 dist_average = dist_average + acos(arg) 00735 00736 arg = sinlat_grid(i)*sinlat_grid(i+1) + 00737 & coslat_grid(i)*coslat_grid(i+1)* 00738 & (coslon_grid(i)*coslon_grid(i+1) + 00739 & sinlon_grid(i)*sinlon_grid(i+1)) 00740 IF (arg < -1.0d0) THEN 00741 arg = -1.0d0 00742 ELSE IF (arg > 1.0d0) THEN 00743 arg = 1.0d0 00744 END IF 00745 dist_average = dist_average + acos(arg) 00746 00747 ELSE IF (i .eq. grid_size) THEN 00748 arg = sinlat_grid(i-1)*sinlat_grid(i) + 00749 & coslat_grid(i-1)*coslat_grid(i)* 00750 & (coslon_grid(i-1)*coslon_grid(i) + 00751 & sinlon_grid(i-1)*sinlon_grid(i)) 00752 IF (arg < -1.0d0) THEN 00753 arg = -1.0d0 00754 ELSE IF (arg > 1.0d0) THEN 00755 arg = 1.0d0 00756 END IF 00757 dist_average = dist_average + acos(arg) 00758 00759 arg = sinlat_grid(i)*sinlat_grid(1) + 00760 & coslat_grid(i)*coslat_grid(1)* 00761 & (coslon_grid(i)*coslon_grid(1) + 00762 & sinlon_grid(i)*sinlon_grid(1)) 00763 IF (arg < -1.0d0) THEN 00764 arg = -1.0d0 00765 ELSE IF (arg > 1.0d0) THEN 00766 arg = 1.0d0 00767 END IF 00768 dist_average = dist_average + acos(arg) 00769 00770 ELSE 00771 arg = sinlat_grid(i-1)*sinlat_grid(i) + 00772 & coslat_grid(i-1)*coslat_grid(i)* 00773 & (coslon_grid(i-1)*coslon_grid(i) + 00774 & sinlon_grid(i-1)*sinlon_grid(i)) 00775 IF (arg < -1.0d0) THEN 00776 arg = -1.0d0 00777 ELSE IF (arg > 1.0d0) THEN 00778 arg = 1.0d0 00779 END IF 00780 dist_average = dist_average + acos(arg) 00781 00782 arg = sinlat_grid(i)*sinlat_grid(i+1) + 00783 & coslat_grid(i)*coslat_grid(i+1)* 00784 & (coslon_grid(i)*coslon_grid(i+1) + 00785 & sinlon_grid(i)*sinlon_grid(i+1)) 00786 IF (arg < -1.0d0) THEN 00787 arg = -1.0d0 00788 ELSE IF (arg > 1.0d0) THEN 00789 arg = 1.0d0 00790 END IF 00791 dist_average = dist_average + acos(arg) 00792 END IF 00793 END DO 00794 dist_average = dist_average / (2*grid_size) 00795 ! 00796 IF (nlogprt .GE. 2) THEN 00797 WRITE (UNIT = nulou,FMT = *)' ' 00798 WRITE (UNIT = nulou,FMT = *) 00799 & 'Leaving routine remap_dist_average' 00800 CALL FLUSH(nulou) 00801 ENDIF 00802 ! 00803 END subroutine grid_dist_average 00804 00805 !*********************************************************************** 00806 00807 end module remap_gaussian_weight 00808 00809 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00810
1.7.4 
