- Timestamp:
- 2016-04-07T16:32:24+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/TOOLS/SIREN/src/grid.f90
r5037 r6440 80 80 !> point:<br/> 81 81 !> @code 82 !> il_index(:)=grid_get_closest(dd_lon0(:,:), dd_lat0(:,:), dd_lon1, dd_lat1) 82 !> il_index(:)=grid_get_closest(dd_lon0(:,:), dd_lat0(:,:), dd_lon1, dd_lat1 83 !> [,dd_fill] [,cd_pos]) 83 84 !> @endcode 84 85 !> - il_index(:) is coarse grid indices (/ i0, j0 /) … … 87 88 !> - dd_lon1 is fine grid longitude value (real(8)) 88 89 !> - dd_lat1 is fine grid latitude value (real(8)) 90 !> - dd_fill 91 !> - cd_pos 89 92 !> 90 93 !> to compute distance between a point A and grid points:<br/> … … 149 152 !> CALL grid_check_coincidence(td_coord0, td_coord1, 150 153 !> id_imin0, id_imax0, id_jmin0, id_jmax0 151 !> [,id_rho])154 !> ,id_rho) 152 155 !> @endcode 153 156 !> - td_coord0 is coarse grid coordinate mpp structure … … 161 164 !> - id_jmax0 is coarse grid upper right corner j-indice of fine grid 162 165 !> domain 163 !> - id_rho is array of refinement factor (default 1)166 !> - id_rho is array of refinement factor 164 167 !> 165 168 !> to add ghost cell at boundaries:<br/> … … 213 216 !> @date October, 2014 214 217 !> - use mpp file structure instead of file 218 !> @date February, 2015 219 !> - add function grid_fill_small_msk to fill small domain inside bigger one 220 !> @February, 2016 221 !> - improve way to check coincidence (bug fix) 222 !> - manage grid cases for T,U,V or F point, with even or odd refinment (bug fix) 215 223 ! 216 224 !> @note Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 255 263 PUBLIC :: grid_split_domain !< compute closed sea domain 256 264 PUBLIC :: grid_fill_small_dom !< fill small closed sea with fill value 265 PUBLIC :: grid_fill_small_msk !< fill small domain inside bigger one 257 266 258 267 ! get closest coarse grid indices of fine grid domain … … 352 361 !> @note need all processor files to be there 353 362 !> @author J.Paul 354 !> - October, 2014- Initial Version363 !> @date October, 2014 - Initial Version 355 364 !> 356 365 !> @param[inout] td_file file structure … … 466 475 !> - compute East West overlap 467 476 !> 468 !> @note need all processor files to be there477 !> @note need all processor files 469 478 !> @author J.Paul 470 !> - October, 2014- Initial Version479 !> @date October, 2014 - Initial Version 471 480 !> 472 481 !> @param[in] td_mpp mpp structure … … 496 505 il_ew =-1 497 506 507 CALL logger_info("GRID GET INFO: look for "//TRIM(td_mpp%c_name)) 498 508 ! copy structure 499 509 tl_mpp=mpp_copy(td_mpp) … … 523 533 ENDIF 524 534 535 CALL logger_info("GRID GET INFO: perio "//TRIM(fct_str(il_perio))) 536 525 537 SELECT CASE(il_perio) 526 538 CASE(3,4) 539 il_pivot=1 540 CASE(5,6) 527 541 il_pivot=0 528 CASE(5,6)529 il_pivot=1530 542 CASE(0,1,2) 531 543 il_pivot=1 … … 534 546 IF( il_pivot < 0 .OR. il_pivot > 1 )THEN 535 547 ! get pivot 548 CALL logger_info("GRID GET INFO: look for pivot ") 536 549 il_pivot=grid_get_pivot(tl_mpp) 537 550 ENDIF … … 539 552 IF( il_perio < 0 .OR. il_perio > 6 )THEN 540 553 ! get periodicity 554 CALL logger_info("GRID GET INFO: look for perio ") 541 555 il_perio=grid_get_perio(tl_mpp, il_pivot) 542 556 ENDIF … … 544 558 IF( il_ew < 0 )THEN 545 559 ! get periodicity 560 CALL logger_info("GRID GET INFO: look for overlap ") 546 561 il_ew=grid_get_ew_overlap(tl_mpp) 547 562 ENDIF … … 595 610 !> 596 611 !> @author J.Paul 597 !> - November, 2013- Subroutine written612 !> @date November, 2013 - Initial version 598 613 !> @date September, 2014 599 614 !> - add dummy loop in case variable not over right point. … … 655 670 656 671 ! no pivot point found 657 CALL logger_ error("GRID GET PIVOT: something wrong "//&672 CALL logger_warn("GRID GET PIVOT: something wrong "//& 658 673 & "when computing pivot point with variable "//& 659 674 & TRIM(td_var%c_name)) … … 676 691 677 692 IF( grid__get_pivot_var /= -1 )THEN 678 CALL logger_ warn("GRID GET PIVOT: variable "//&693 CALL logger_info("GRID GET PIVOT: variable "//& 679 694 & TRIM(td_var%c_name)//" seems to be on grid point "//& 680 695 & TRIM(cp_grid_point(jj)) ) … … 708 723 !> 709 724 !> @author J.Paul 710 !> -October, 2014 - Initial version725 !> @date October, 2014 - Initial version 711 726 ! 712 727 !> @param[in] dd_value array of value … … 783 798 784 799 IF( ll_check )THEN 785 CALL logger_info("GRID GET PIVOT: T-pivot")800 CALL logger_info("GRID GET PIVOT: F-pivot") 786 801 grid__get_pivot_varT=0 787 802 ENDIF … … 801 816 !> 802 817 !> @author J.Paul 803 !> -October, 2014 - Initial version818 !> @date October, 2014 - Initial version 804 819 ! 805 820 !> @param[in] dd_value array of value … … 876 891 877 892 IF( ll_check )THEN 878 CALL logger_info("GRID GET PIVOT: T-pivot")893 CALL logger_info("GRID GET PIVOT: F-pivot") 879 894 grid__get_pivot_varU=0 880 895 ENDIF … … 894 909 !> 895 910 !> @author J.Paul 896 !> -October, 2014 - Initial version911 !> @date October, 2014 - Initial version 897 912 ! 898 913 !> @param[in] dd_value array of value … … 969 984 970 985 IF( ll_check )THEN 971 CALL logger_info("GRID GET PIVOT: T-pivot")986 CALL logger_info("GRID GET PIVOT: F-pivot") 972 987 grid__get_pivot_varV=0 973 988 ENDIF … … 987 1002 !> 988 1003 !> @author J.Paul 989 !> -October, 2014 - Initial version1004 !> @date October, 2014 - Initial version 990 1005 ! 991 1006 !> @param[in] dd_value array of value … … 1062 1077 1063 1078 IF( ll_check )THEN 1064 CALL logger_info("GRID GET PIVOT: T-pivot")1079 CALL logger_info("GRID GET PIVOT: F-pivot") 1065 1080 grid__get_pivot_varF=0 1066 1081 ENDIF … … 1083 1098 !> 1084 1099 !> @author J.Paul 1085 !> - Ocotber, 2014- Initial version1100 !> @date Ocotber, 2014 - Initial version 1086 1101 ! 1087 1102 !> @param[in] td_file file structure … … 1172 1187 !> 1173 1188 !> @author J.Paul 1174 !> -October, 2014 - Initial version1189 !> @date October, 2014 - Initial version 1175 1190 ! 1176 1191 !> @param[in] td_mpp mpp file structure … … 1277 1292 !> 1: cyclic east-west boundary 1278 1293 !> 2: symmetric boundary condition across the equator 1279 !> 3: North fold boundary (with a F-point pivot)1280 !> 4: North fold boundary (with a F-point pivot) and cyclic east-west boundary1281 !> 5: North fold boundary (with a T-point pivot)1282 !> 6: North fold boundary (with a T-point pivot) and cyclic east-west boundary1294 !> 3: North fold boundary (with a T-point pivot) 1295 !> 4: North fold boundary (with a T-point pivot) and cyclic east-west boundary 1296 !> 5: North fold boundary (with a F-point pivot) 1297 !> 6: North fold boundary (with a F-point pivot) and cyclic east-west boundary 1283 1298 !> 1284 1299 !> @warning pivot point should have been computed before run this script. see grid_get_pivot. 1285 1300 !> 1286 1301 !> @author J.Paul 1287 !> - November, 2013- Subroutine written1302 !> @date November, 2013 - Initial version 1288 1303 !> @date October, 2014 1289 1304 !> - work on variable structure instead of file structure … … 1326 1341 il_dim(:)=td_var%t_dim(:)%i_len 1327 1342 1328 CALL logger_ info("GRID GET PERIO: use varibale "//TRIM(td_var%c_name))1329 CALL logger_ info("GRID GET PERIO: fill value "//TRIM(fct_str(td_var%d_fill)))1330 CALL logger_ info("GRID GET PERIO: fillvalue "//TRIM(fct_str(td_var%d_value(1,1,1,1))))1343 CALL logger_debug("GRID GET PERIO: use varibale "//TRIM(td_var%c_name)) 1344 CALL logger_debug("GRID GET PERIO: fill value "//TRIM(fct_str(td_var%d_fill))) 1345 CALL logger_debug("GRID GET PERIO: first value "//TRIM(fct_str(td_var%d_value(1,1,1,1)))) 1331 1346 1332 1347 IF(ALL(td_var%d_value( 1 , : ,1,1)/=td_var%d_fill).AND.& … … 1335 1350 & ALL(td_var%d_value( : ,il_dim(2),1,1)/=td_var%d_fill))THEN 1336 1351 ! no boundary closed 1337 CALL logger_ warn("GRID GET PERIO: can't determined periodicity. "//&1352 CALL logger_error("GRID GET PERIO: can't determined periodicity. "//& 1338 1353 & "there is no boundary closed for variable "//& 1339 1354 & TRIM(td_var%c_name) ) 1355 ! check pivot 1356 SELECT CASE(id_pivot) 1357 CASE(0) 1358 ! F pivot 1359 CALL logger_warn("GRID GET PERIO: assume domain is global") 1360 grid__get_perio_var=6 1361 CASE(1) 1362 ! T pivot 1363 CALL logger_warn("GRID GET PERIO: assume domain is global") 1364 grid__get_perio_var=4 1365 END SELECT 1340 1366 ELSE 1341 1367 ! check periodicity … … 1452 1478 !> 1453 1479 !> @author J.Paul 1454 !> -October, 2014 - Initial version1480 !> @date October, 2014 - Initial version 1455 1481 !> 1456 1482 !> @param[in] td_file file structure … … 1537 1563 !> 1: cyclic east-west boundary 1538 1564 !> 2: symmetric boundary condition across the equator 1539 !> 3: North fold boundary (with a F-point pivot)1540 !> 4: North fold boundary (with a F-point pivot) and cyclic east-west boundary1541 !> 5: North fold boundary (with a T-point pivot)1542 !> 6: North fold boundary (with a T-point pivot) and cyclic east-west boundary1565 !> 3: North fold boundary (with a T-point pivot) 1566 !> 4: North fold boundary (with a T-point pivot) and cyclic east-west boundary 1567 !> 5: North fold boundary (with a F-point pivot) 1568 !> 6: North fold boundary (with a F-point pivot) and cyclic east-west boundary 1543 1569 !> 1544 1570 !> @warning pivot point should have been computed before run this script. see grid_get_pivot. 1545 1571 !> 1546 1572 !> @author J.Paul 1547 !> -October, 2014 - Initial version1573 !> @date October, 2014 - Initial version 1548 1574 ! 1549 1575 !> @param[in] td_mpp mpp file structure … … 1634 1660 ! 1635 1661 !> @author J.Paul 1636 !> - November, 2013- Initial Version1662 !> @date November, 2013 - Initial Version 1637 1663 !> @date October, 2014 1638 1664 !> - work on mpp file structure instead of file structure … … 1746 1772 !> 1747 1773 !> @author J.Paul 1748 !> - October, 2014- Initial Version1774 !> @date October, 2014 - Initial Version 1749 1775 !> 1750 1776 !> @param[in] td_file file structure … … 1797 1823 ! 1798 1824 !> @author J.Paul 1799 !> - November, 2013- Initial Version1825 !> @date November, 2013 - Initial Version 1800 1826 !> @date October, 2014 1801 1827 !> - work on mpp file structure instead of file structure … … 1853 1879 !> 1854 1880 !> @author J.Paul 1855 !> - November, 2013- Initial Version1881 !> @date November, 2013 - Initial Version 1856 1882 !> 1857 1883 !> @param[in] td_lat latitude variable structure … … 1890 1916 ! 1891 1917 !> @author J.Paul 1892 !> - November, 2013- Initial Version1918 !> @date November, 2013 - Initial Version 1893 1919 !> @date October, 2014 1894 1920 !> - work on mpp file structure instead of file structure … … 1978 2004 !> 1979 2005 !> @author J.Paul 1980 !> - November, 2013- Initial Version2006 !> @date November, 2013 - Initial Version 1981 2007 !> @date September, 2014 1982 2008 !> - use grid point to read coordinates variable. 1983 2009 !> @date October, 2014 1984 2010 !> - work on mpp file structure instead of file structure 2011 !> @date February, 2015 2012 !> - use longitude or latitude as standard name, if can not find 2013 !> longitude_T, latitude_T... 1985 2014 !> 1986 2015 !> @param[in] td_coord0 coarse grid coordinate mpp structure … … 2004 2033 2005 2034 ! local variable 2006 TYPE(TMPP) :: tl_coord0 2007 TYPE(TMPP) :: tl_coord1 2008 2009 TYPE(TVAR) :: tl_lon0 2010 TYPE(TVAR) :: tl_lat0 2011 TYPE(TVAR) :: tl_lon1 2012 TYPE(TVAR) :: tl_lat1 2013 2014 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 2015 2016 INTEGER(i4), DIMENSION(2,2) :: il_xghost0 2017 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 2018 2019 INTEGER(i4) :: il_imin0 2020 INTEGER(i4) :: il_imax0 2021 INTEGER(i4) :: il_jmin0 2022 INTEGER(i4) :: il_jmax0 2023 2024 CHARACTER(LEN= 1) :: cl_point 2025 CHARACTER(LEN=lc) :: cl_name 2035 CHARACTER(LEN= 1) :: cl_point 2036 CHARACTER(LEN=lc) :: cl_name 2037 2038 INTEGER(i4) :: il_imin0 2039 INTEGER(i4) :: il_imax0 2040 INTEGER(i4) :: il_jmin0 2041 INTEGER(i4) :: il_jmax0 2042 INTEGER(i4) :: il_ind 2043 2044 INTEGER(i4), DIMENSION(2,2) :: il_xghost0 2045 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 2046 2047 INTEGER(i4), DIMENSION(:) , ALLOCATABLE :: il_rho 2048 2049 TYPE(TVAR) :: tl_lon0 2050 TYPE(TVAR) :: tl_lat0 2051 TYPE(TVAR) :: tl_lon1 2052 TYPE(TVAR) :: tl_lat1 2053 2054 TYPE(TMPP) :: tl_coord0 2055 TYPE(TMPP) :: tl_coord1 2026 2056 2027 2057 ! loop indices … … 2057 2087 ! read coarse longitue and latitude 2058 2088 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 2089 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 2090 IF( il_ind == 0 )THEN 2091 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2092 & TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". & 2093 & try to use longitude.") 2094 WRITE(cl_name,*) 'longitude' 2095 ENDIF 2059 2096 tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 2097 2060 2098 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 2099 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 2100 IF( il_ind == 0 )THEN 2101 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2102 & TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". & 2103 & try to use latitude.") 2104 WRITE(cl_name,*) 'latitude' 2105 ENDIF 2061 2106 tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 2062 2107 … … 2077 2122 ! read fine longitue and latitude 2078 2123 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 2124 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 2125 IF( il_ind == 0 )THEN 2126 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2127 & TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". & 2128 & try to use longitude.") 2129 WRITE(cl_name,*) 'longitude' 2130 ENDIF 2079 2131 tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 2132 2080 2133 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 2134 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 2135 IF( il_ind == 0 )THEN 2136 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2137 & TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". & 2138 & try to use latitude.") 2139 WRITE(cl_name,*) 'latitude' 2140 ENDIF 2081 2141 tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 2082 2142 … … 2127 2187 !> 2128 2188 !> @author J.Paul 2129 !> - November, 2013- Initial Version2189 !> @date November, 2013 - Initial Version 2130 2190 !> @date September, 2014 2131 2191 !> - use grid point to read coordinates variable. 2132 2192 !> @date October, 2014 2133 2193 !> - work on mpp file structure instead of file structure 2194 !> @date February, 2015 2195 !> - use longitude or latitude as standard name, if can not find 2196 !> longitude_T, latitude_T... 2134 2197 !> 2135 2198 !> @param[in] td_longitude0 coarse grid longitude … … 2154 2217 2155 2218 ! local variable 2156 TYPE(TMPP) :: tl_coord1 2157 2158 TYPE(TVAR) :: tl_lon1 2159 TYPE(TVAR) :: tl_lat1 2160 2161 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 2162 2163 INTEGER(i4), DIMENSION(2,2) :: il_xghost 2164 2165 CHARACTER(LEN= 1) :: cl_point 2166 CHARACTER(LEN=lc) :: cl_name 2219 CHARACTER(LEN= 1) :: cl_point 2220 CHARACTER(LEN=lc) :: cl_name 2221 2222 INTEGER(i4) :: il_ind 2223 2224 INTEGER(i4), DIMENSION(:) , ALLOCATABLE :: il_rho 2225 2226 INTEGER(i4), DIMENSION(2,2) :: il_xghost 2227 2228 TYPE(TVAR) :: tl_lon1 2229 TYPE(TVAR) :: tl_lat1 2230 2231 TYPE(TMPP) :: tl_coord1 2167 2232 2168 2233 ! loop indices … … 2209 2274 ! read fine longitue and latitude 2210 2275 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 2276 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 2277 IF( il_ind == 0 )THEN 2278 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2279 & TRIM(cl_name)//"in file "//TRIM(tl_coord1%c_name)//". & 2280 & try to use longitude.") 2281 WRITE(cl_name,*) 'longitude' 2282 ENDIF 2211 2283 tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 2284 2212 2285 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 2286 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 2287 IF( il_ind == 0 )THEN 2288 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2289 & TRIM(cl_name)//"in file "//TRIM(tl_coord1%c_name)//". & 2290 & try to use longitude.") 2291 WRITE(cl_name,*) 'latitude' 2292 ENDIF 2213 2293 tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 2214 2294 … … 2224 2304 & il_rho(:), cl_point ) 2225 2305 2226 2227 2306 CALL var_clean(tl_lon1) 2228 2307 CALL var_clean(tl_lat1) … … 2244 2323 !> 2245 2324 !> @author J.Paul 2246 !> - November, 2013- Initial Version2325 !> @date November, 2013 - Initial Version 2247 2326 !> @date September, 2014 2248 2327 !> - use grid point to read coordinates variable. 2249 2328 !> @date October, 2014 2250 2329 !> - work on mpp file structure instead of file structure 2330 !> @date February, 2015 2331 !> - use longitude or latitude as standard name, if can not find 2332 !> longitude_T, latitude_T... 2251 2333 !> 2252 2334 !> @param[in] td_coord0 coarse grid coordinate mpp structure … … 2271 2353 2272 2354 ! local variable 2273 TYPE(TMPP) :: tl_coord0 2274 2275 TYPE(TVAR) :: tl_lon0 2276 TYPE(TVAR) :: tl_lat0 2277 2278 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 2279 2280 INTEGER(i4), DIMENSION(2,2) :: il_xghost 2281 2282 INTEGER(i4) :: il_imin0 2283 INTEGER(i4) :: il_imax0 2284 INTEGER(i4) :: il_jmin0 2285 INTEGER(i4) :: il_jmax0 2286 2287 CHARACTER(LEN= 1) :: cl_point 2288 CHARACTER(LEN=lc) :: cl_name 2355 CHARACTER(LEN= 1) :: cl_point 2356 CHARACTER(LEN=lc) :: cl_name 2357 2358 INTEGER(i4) :: il_imin0 2359 INTEGER(i4) :: il_imax0 2360 INTEGER(i4) :: il_jmin0 2361 INTEGER(i4) :: il_jmax0 2362 INTEGER(i4) :: il_ind 2363 2364 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 2365 2366 INTEGER(i4), DIMENSION(2,2) :: il_xghost 2367 2368 TYPE(TVAR) :: tl_lon0 2369 TYPE(TVAR) :: tl_lat0 2370 2371 TYPE(TMPP) :: tl_coord0 2289 2372 2290 2373 ! loop indices … … 2330 2413 ! read coarse longitue and latitude 2331 2414 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 2415 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 2416 IF( il_ind == 0 )THEN 2417 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2418 & TRIM(cl_name)//"in file "//TRIM(tl_coord0%c_name)//". & 2419 & try to use longitude.") 2420 WRITE(cl_name,*) 'longitude' 2421 ENDIF 2332 2422 tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 2423 2333 2424 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 2425 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 2426 IF( il_ind == 0 )THEN 2427 CALL logger_warn("GRID GET COARSE INDEX: no variable "//& 2428 & TRIM(cl_name)//"in file "//TRIM(tl_coord0%c_name)//". & 2429 & try to use latitude.") 2430 WRITE(cl_name,*) 'latitude' 2431 ENDIF 2334 2432 tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 2335 2433 … … 2377 2475 !> 2378 2476 !> @author J.Paul 2379 !> - November, 2013- Initial Version2477 !> @date November, 2013 - Initial Version 2380 2478 !> @date September, 2014 2381 2479 !> - check grid point 2382 2480 !> - take into account EW overlap 2481 !> @date February, 2016 2482 !> - use delta (lon or lat) 2483 !> - manage cases for T,U,V or F point, with even or odd refinment 2383 2484 !> 2384 2485 !> @param[in] td_lon0 coarse grid longitude … … 2408 2509 2409 2510 ! local variable 2410 REAL(dp) :: dl_lon1_ll 2411 REAL(dp) :: dl_lon1_ul 2412 REAL(dp) :: dl_lon1_lr 2413 REAL(dp) :: dl_lon1_ur 2414 2415 REAL(dp) :: dl_lat1_ll 2416 REAL(dp) :: dl_lat1_ul 2417 REAL(dp) :: dl_lat1_lr 2418 REAL(dp) :: dl_lat1_ur 2511 CHARACTER(LEN= 1) :: cl_point0 2512 CHARACTER(LEN= 1) :: cl_point1 2513 2514 LOGICAL , DIMENSION(2) :: ll_even 2515 2516 REAL(dp) :: dl_lon1 2517 REAL(dp) :: dl_dlon 2518 REAL(dp) :: dl_lat1 2519 REAL(dp) :: dl_dlat 2520 2521 INTEGER(i4) :: il_ew0 2522 INTEGER(i4) :: il_imin0 2523 INTEGER(i4) :: il_imax0 2524 INTEGER(i4) :: il_jmin0 2525 INTEGER(i4) :: il_jmax0 2526 2527 INTEGER(i4) :: il_ew1 2528 INTEGER(i4) :: il_imin1 2529 INTEGER(i4) :: il_imax1 2530 INTEGER(i4) :: il_jmin1 2531 INTEGER(i4) :: il_jmax1 2532 2533 INTEGER(i4) :: il_imin 2534 INTEGER(i4) :: il_imax 2535 INTEGER(i4) :: il_jmin 2536 INTEGER(i4) :: il_jmax 2419 2537 2420 2538 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 2421 2539 2422 INTEGER(i4), DIMENSION(2) :: il_ill 2423 INTEGER(i4), DIMENSION(2) :: il_ilr 2424 INTEGER(i4), DIMENSION(2) :: il_iul 2425 INTEGER(i4), DIMENSION(2) :: il_iur 2426 2427 INTEGER(i4) :: il_ew0 2428 INTEGER(i4) :: il_imin0 2429 INTEGER(i4) :: il_imax0 2430 INTEGER(i4) :: il_jmin0 2431 INTEGER(i4) :: il_jmax0 2432 2433 INTEGER(i4) :: il_ew1 2434 INTEGER(i4) :: il_imin1 2435 INTEGER(i4) :: il_imax1 2436 INTEGER(i4) :: il_jmin1 2437 INTEGER(i4) :: il_jmax1 2438 2439 INTEGER(i4) :: il_imin 2440 INTEGER(i4) :: il_imax 2441 INTEGER(i4) :: il_jmin 2442 INTEGER(i4) :: il_jmax 2443 2444 INTEGER(i4), DIMENSION(2,2) :: il_xghost0 2445 INTEGER(i4), DIMENSION(2,2) :: il_yghost0 2446 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 2447 INTEGER(i4), DIMENSION(2,2) :: il_yghost1 2448 2449 TYPE(TVAR) :: tl_lon0 2450 TYPE(TVAR) :: tl_lat0 2451 TYPE(TVAR) :: tl_lon1 2452 TYPE(TVAR) :: tl_lat1 2453 2454 CHARACTER(LEN= 1) :: cl_point0 2455 CHARACTER(LEN= 1) :: cl_point1 2456 2540 INTEGER(i4), DIMENSION(2) :: il_ill 2541 INTEGER(i4), DIMENSION(2) :: il_ilr 2542 INTEGER(i4), DIMENSION(2) :: il_iul 2543 INTEGER(i4), DIMENSION(2) :: il_iur 2544 2545 INTEGER(i4), DIMENSION(2,2) :: il_xghost0 2546 INTEGER(i4), DIMENSION(2,2) :: il_yghost0 2547 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 2548 INTEGER(i4), DIMENSION(2,2) :: il_yghost1 2549 2550 TYPE(TVAR) :: tl_lon0 2551 TYPE(TVAR) :: tl_lat0 2552 TYPE(TVAR) :: tl_lon1 2553 TYPE(TVAR) :: tl_lat1 2554 2457 2555 ! loop indices 2458 INTEGER(i4) :: ji2459 INTEGER(i4) :: jj2460 2556 !---------------------------------------------------------------- 2461 2557 ! init … … 2465 2561 il_rho(:)=1 2466 2562 IF( PRESENT(id_rho) ) il_rho(:)=id_rho(:) 2563 2564 ll_even(:)=(/ (MOD(id_rho(jp_I),2)==0), (MOD(id_rho(jp_J),2)==0) /) 2467 2565 2468 2566 cl_point0='T' … … 2520 2618 CALL grid_del_ghost(tl_lon0, il_xghost0(:,:)) 2521 2619 CALL grid_del_ghost(tl_lat0, il_xghost0(:,:)) 2522 2620 2523 2621 ! "global" coarse grid indice 2524 2622 il_imin0=1 … … 2563 2661 ! get indices for each corner 2564 2662 !1- search lower left corner indices 2565 dl_lon1 _ll=tl_lon1%d_value( il_imin1, il_jmin1, 1, 1 )2566 dl_lat1 _ll=tl_lat1%d_value( il_imin1, il_jmin1, 1, 1 )2567 2568 IF( dl_lon1 _ll== tl_lon1%d_fill .OR. &2569 & dl_lat1 _ll== tl_lat1%d_fill )THEN2663 dl_lon1=tl_lon1%d_value( il_imin1, il_jmin1, 1, 1 ) 2664 dl_lat1=tl_lat1%d_value( il_imin1, il_jmin1, 1, 1 ) 2665 2666 IF( dl_lon1 == tl_lon1%d_fill .OR. & 2667 & dl_lat1 == tl_lat1%d_fill )THEN 2570 2668 CALL logger_error("GRID GET COARSE INDEX: lower left corner "//& 2571 2669 & "point is FillValue. remove ghost cell "//& 2572 2670 & "before running grid_get_coarse_index.") 2573 2671 ENDIF 2672 2673 !!!!! i-direction !!!!! 2674 IF( ll_even(jp_I) )THEN 2675 ! even 2676 SELECT CASE(TRIM(cl_point1)) 2677 CASE('F','U') 2678 dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmin1,1,1) - & 2679 & tl_lon1%d_value(il_imin1 ,il_jmin1,1,1) ) / & 2680 & 2. 2681 CASE DEFAULT 2682 dl_dlon=0 2683 END SELECT 2684 ELSE 2685 ! odd 2686 dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmin1,1,1) - & 2687 & tl_lon1%d_value(il_imin1 ,il_jmin1,1,1) ) / & 2688 & 2. 2689 ENDIF 2690 2691 !!!!! j-direction !!!!! 2692 IF( ll_even(jp_J) )THEN 2693 ! even 2694 SELECT CASE(TRIM(cl_point1)) 2695 CASE('F','V') 2696 dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmin1+1,1,1) - & 2697 & tl_lat1%d_value(il_imin1,il_jmin1 ,1,1) ) / & 2698 & 2. 2699 CASE DEFAULT 2700 dl_dlat=0 2701 END SELECT 2702 ELSE 2703 ! odd 2704 dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmin1+1,1,1) - & 2705 & tl_lat1%d_value(il_imin1,il_jmin1 ,1,1) ) / & 2706 & 2. 2707 ENDIF 2708 2709 dl_lon1 = dl_lon1 + dl_dlon 2710 dl_lat1 = dl_lat1 + dl_dlat 2711 2574 2712 ! look for closest point on coarse grid 2575 2713 il_ill(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & … … 2579 2717 & il_jmin0:il_jmax0, & 2580 2718 & 1,1), & 2581 & dl_lon1_ll, dl_lat1_ll ) 2582 2583 ! coarse grid point should be south west of fine grid domain 2584 ji = il_ill(1) 2585 jj = il_ill(2) 2586 2587 IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_ll) > dp_delta )THEN 2588 IF(tl_lon0%d_value(ji,jj,1,1) > dl_lon1_ll )THEN 2589 il_ill(1)=il_ill(1)-1 2590 IF( il_ill(1) <= 0 )THEN 2591 IF( tl_lon0%i_ew >= 0 )THEN 2592 il_ill(1)=tl_lon0%t_dim(jp_I)%i_len-tl_lon0%i_ew 2593 ELSE 2594 CALL logger_error("GRID GET COARSE INDEX: error "//& 2595 & "computing lower left corner "//& 2596 & "index for longitude") 2597 ENDIF 2598 ENDIF 2599 ENDIF 2600 ENDIF 2601 IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_ll) > dp_delta )THEN 2602 IF(tl_lat0%d_value(ji,jj,1,1) > dl_lat1_ll )THEN 2603 il_ill(2)=il_ill(2)-1 2604 IF( il_ill(2)-1 <= 0 )THEN 2605 CALL logger_error("GRID GET COARSE INDEX: error "//& 2606 & "computing lower left corner "//& 2607 & "index for latitude") 2608 ENDIF 2609 ENDIF 2610 ENDIF 2719 & dl_lon1, dl_lat1, 'll' ) 2720 2611 2721 2612 2722 !2- search upper left corner indices 2613 dl_lon1 _ul=tl_lon1%d_value( il_imin1, il_jmax1, 1, 1 )2614 dl_lat1 _ul=tl_lat1%d_value( il_imin1, il_jmax1, 1, 1 )2615 2616 IF( dl_lon1 _ul== tl_lon1%d_fill .OR. &2617 & dl_lat1 _ul== tl_lat1%d_fill )THEN2723 dl_lon1=tl_lon1%d_value( il_imin1, il_jmax1, 1, 1 ) 2724 dl_lat1=tl_lat1%d_value( il_imin1, il_jmax1, 1, 1 ) 2725 2726 IF( dl_lon1 == tl_lon1%d_fill .OR. & 2727 & dl_lat1 == tl_lat1%d_fill )THEN 2618 2728 CALL logger_error("GRID GET COARSE INDEX: upper left corner "//& 2619 2729 & "point is FillValue. remove ghost cell "//& 2620 2730 & "running grid_get_coarse_index.") 2621 2731 ENDIF 2732 2733 !!!!! i-direction !!!!! 2734 IF( ll_even(jp_I) )THEN 2735 ! even 2736 SELECT CASE(TRIM(cl_point1)) 2737 CASE('F','U') 2738 dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmax1,1,1) - & 2739 & tl_lon1%d_value(il_imin1 ,il_jmax1,1,1) ) / & 2740 & 2. 2741 CASE DEFAULT 2742 dl_dlon=0 2743 END SELECT 2744 ELSE 2745 ! odd 2746 dl_dlon= ( tl_lon1%d_value(il_imin1+1,il_jmax1,1,1) - & 2747 & tl_lon1%d_value(il_imin1 ,il_jmax1,1,1) ) / & 2748 & 2. 2749 ENDIF 2750 2751 !!!!! j-direction !!!!! 2752 IF( ll_even(jp_J) )THEN 2753 ! even 2754 SELECT CASE(TRIM(cl_point1)) 2755 CASE('F','V') 2756 dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmax1 ,1,1) - & 2757 & tl_lat1%d_value(il_imin1,il_jmax1-1,1,1) ) / & 2758 & 2. 2759 CASE DEFAULT 2760 dl_dlat=0 2761 END SELECT 2762 ELSE 2763 ! odd 2764 dl_dlat= ( tl_lat1%d_value(il_imin1,il_jmax1 ,1,1) - & 2765 & tl_lat1%d_value(il_imin1,il_jmax1-1,1,1) ) / & 2766 & 2. 2767 ENDIF 2768 2769 dl_lon1 = dl_lon1 + dl_dlon 2770 dl_lat1 = dl_lat1 - dl_dlat 2771 2622 2772 ! look for closest point on coarse grid 2623 2773 il_iul(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & … … 2627 2777 & il_jmin0:il_jmax0, & 2628 2778 & 1,1), & 2629 & dl_lon1_ul, dl_lat1_ul ) 2630 2631 ! coarse grid point should be north west of fine grid domain 2632 ji = il_iul(1) 2633 jj = il_iul(2) 2634 2635 IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_ul) > dp_delta )THEN 2636 IF(tl_lon0%d_value(ji,jj,1,1) > dl_lon1_ul )THEN 2637 il_iul(1)=il_iul(1)-1 2638 IF( il_iul(1) <= 0 )THEN 2639 IF( tl_lon0%i_ew >= 0 )THEN 2640 il_iul(1)=tl_lon0%t_dim(jp_I)%i_len-tl_lon0%i_ew 2641 ELSE 2642 CALL logger_error("GRID GET COARSE INDEX: error "//& 2643 & "computing upper left corner "//& 2644 & "index for longitude") 2645 ENDIF 2646 ENDIF 2647 ENDIF 2648 ENDIF 2649 2650 IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_ul) > dp_delta )THEN 2651 IF(tl_lat0%d_value(ji,jj,1,1) < dl_lat1_ul )THEN 2652 il_iul(2)=il_iul(2)+1 2653 IF( il_ill(2) > tl_lat0%t_dim(jp_J)%i_len )THEN 2654 CALL logger_error("GRID GET COARSE INDEX: error "//& 2655 & "computing upper left corner "//& 2656 & "index for latitude") 2657 ENDIF 2658 ENDIF 2659 ENDIF 2779 & dl_lon1, dl_lat1, 'ul' ) 2660 2780 2661 2781 !3- search lower right corner indices 2662 dl_lon1 _lr=tl_lon1%d_value( il_imax1, il_jmin1, 1, 1 )2663 dl_lat1 _lr=tl_lat1%d_value( il_imax1, il_jmin1, 1, 1 )2664 2665 IF( dl_lon1 _lr== tl_lon1%d_fill .OR. &2666 & dl_lat1 _lr== tl_lat1%d_fill )THEN2782 dl_lon1=tl_lon1%d_value( il_imax1, il_jmin1, 1, 1 ) 2783 dl_lat1=tl_lat1%d_value( il_imax1, il_jmin1, 1, 1 ) 2784 2785 IF( dl_lon1 == tl_lon1%d_fill .OR. & 2786 & dl_lat1 == tl_lat1%d_fill )THEN 2667 2787 CALL logger_error("GRID GET COARSE INDEX: lower right corner "//& 2668 2788 & "point is FillValue. remove ghost cell "//& 2669 2789 & "running grid_get_coarse_index.") 2670 2790 ENDIF 2791 2792 !!!!! i-direction !!!!! 2793 IF( ll_even(jp_I) )THEN 2794 ! even 2795 SELECT CASE(TRIM(cl_point1)) 2796 CASE('F','U') 2797 dl_dlon= ( tl_lon1%d_value(il_imax1 ,il_jmin1,1,1) - & 2798 & tl_lon1%d_value(il_imax1-1,il_jmin1,1,1) ) / & 2799 & 2. 2800 CASE DEFAULT 2801 dl_dlon=0 2802 END SELECT 2803 ELSE 2804 ! odd 2805 dl_dlon= ( tl_lon1%d_value(il_imax1 ,il_jmin1,1,1) - & 2806 & tl_lon1%d_value(il_imax1-1,il_jmin1,1,1) ) / & 2807 & 2. 2808 ENDIF 2809 2810 !!!!! j-direction !!!!! 2811 IF( ll_even(jp_J) )THEN 2812 ! even 2813 SELECT CASE(TRIM(cl_point1)) 2814 CASE('F','V') 2815 dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmin1+1,1,1) - & 2816 & tl_lat1%d_value(il_imax1,il_jmin1 ,1,1) ) / & 2817 & 2. 2818 CASE DEFAULT 2819 dl_dlat=0 2820 END SELECT 2821 ELSE 2822 ! odd 2823 dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmin1+1,1,1) - & 2824 & tl_lat1%d_value(il_imax1,il_jmin1 ,1,1) ) / & 2825 & 2. 2826 ENDIF 2827 2828 dl_lon1 = dl_lon1 - dl_dlon 2829 dl_lat1 = dl_lat1 + dl_dlat 2830 2671 2831 ! look for closest point on coarse grid 2672 2832 il_ilr(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & … … 2676 2836 & il_jmin0:il_jmax0, & 2677 2837 & 1,1), & 2678 & dl_lon1_lr, dl_lat1_lr ) 2679 2680 ! coarse grid point should be south east of fine grid domain 2681 ji = il_ilr(1) 2682 jj = il_ilr(2) 2683 IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_lr) > dp_delta )THEN 2684 IF( tl_lon0%d_value(ji,jj,1,1) < dl_lon1_lr )THEN 2685 il_ilr(1)=il_ilr(1)+1 2686 IF( il_ilr(1) > tl_lon0%t_dim(jp_I)%i_len )THEN 2687 IF( tl_lon0%i_ew >= 0 )THEN 2688 il_ilr(1)=tl_lon0%i_ew+1 2689 ELSE 2690 CALL logger_error("GRID GET COARSE INDEX: error "//& 2691 & "computing lower right corner "//& 2692 & "index for longitude") 2693 ENDIF 2694 ENDIF 2695 ENDIF 2696 ENDIF 2697 IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_lr) > dp_delta )THEN 2698 IF( tl_lat0%d_value(ji,jj,1,1) > dl_lat1_lr )THEN 2699 il_ilr(2)=il_ilr(2)-1 2700 IF( il_ilr(2) <= 0 )THEN 2701 CALL logger_error("GRID GET COARSE INDEX: error "//& 2702 & "computing lower right corner "//& 2703 & "index for latitude") 2704 ENDIF 2705 ENDIF 2706 ENDIF 2838 & dl_lon1, dl_lat1, 'lr' ) 2707 2839 2708 2840 !4- search upper right corner indices 2709 dl_lon1 _ur=tl_lon1%d_value( il_imax1, il_jmax1, 1, 1 )2710 dl_lat1 _ur=tl_lat1%d_value( il_imax1, il_jmax1, 1, 1 )2711 2712 IF( dl_lon1 _ur== tl_lon1%d_fill .OR. &2713 & dl_lat1 _ur== tl_lat1%d_fill )THEN2841 dl_lon1=tl_lon1%d_value( il_imax1, il_jmax1, 1, 1 ) 2842 dl_lat1=tl_lat1%d_value( il_imax1, il_jmax1, 1, 1 ) 2843 2844 IF( dl_lon1 == tl_lon1%d_fill .OR. & 2845 & dl_lat1 == tl_lat1%d_fill )THEN 2714 2846 CALL logger_error("GRID GET COARSE INDEX: upper right corner "//& 2715 2847 & "point is FillValue. remove ghost cell "//& 2716 & " running grid_get_coarse_index.")2848 & "before running grid_get_coarse_index.") 2717 2849 ENDIF 2850 2851 !!!!! i-direction !!!!! 2852 IF( ll_even(jp_I) )THEN 2853 ! even 2854 SELECT CASE(TRIM(cl_point1)) 2855 CASE('F','U') 2856 dl_dlon= ( tl_lon1%d_value(il_imax1 ,il_jmax1,1,1) - & 2857 & tl_lon1%d_value(il_imax1-1,il_jmax1,1,1) ) / & 2858 & 2. 2859 CASE DEFAULT 2860 dl_dlon=0 2861 END SELECT 2862 ELSE 2863 ! odd 2864 dl_dlon= ( tl_lon1%d_value(il_imax1 ,il_jmax1,1,1) - & 2865 & tl_lon1%d_value(il_imax1-1,il_jmax1,1,1) ) / & 2866 & 2. 2867 ENDIF 2868 2869 !!!!! j-direction !!!!! 2870 IF( ll_even(jp_J) )THEN 2871 ! even 2872 SELECT CASE(TRIM(cl_point1)) 2873 CASE('F','V') 2874 dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmax1 ,1,1) - & 2875 & tl_lat1%d_value(il_imax1,il_jmax1-1,1,1) ) / & 2876 & 2. 2877 CASE DEFAULT 2878 dl_dlat=0 2879 END SELECT 2880 ELSE 2881 ! odd 2882 dl_dlat= ( tl_lat1%d_value(il_imax1,il_jmax1 ,1,1) - & 2883 & tl_lat1%d_value(il_imax1,il_jmax1-1,1,1) ) / & 2884 & 2. 2885 ENDIF 2886 2887 dl_lon1 = dl_lon1 - dl_dlon 2888 dl_lat1 = dl_lat1 - dl_dlat 2889 2718 2890 ! look for closest point on coarse grid 2719 2891 il_iur(:)= grid_get_closest(tl_lon0%d_value(il_imin0:il_imax0, & … … 2723 2895 & il_jmin0:il_jmax0, & 2724 2896 & 1,1), & 2725 & dl_lon1_ur, dl_lat1_ur ) 2726 2727 ! coarse grid point should be north east fine grid domain 2728 ji = il_iur(1) 2729 jj = il_iur(2) 2730 IF( ABS(tl_lon0%d_value(ji,jj,1,1)-dl_lon1_ur) > dp_delta )THEN 2731 IF( tl_lon0%d_value(ji,jj,1,1) < dl_lon1_ur )THEN 2732 il_iur(1)=il_iur(1)+1 2733 IF( il_iur(1) > tl_lon0%t_dim(jp_I)%i_len )THEN 2734 IF( tl_lon0%i_ew >= 0 )THEN 2735 il_iur(1)=tl_lon0%i_ew+1 2736 ELSE 2737 CALL logger_error("GRID GET COARSE INDEX: error "//& 2738 & "computing upper right corner "//& 2739 & "index for longitude") 2740 ENDIF 2741 ENDIF 2742 ENDIF 2743 ENDIF 2744 IF( ABS(tl_lat0%d_value(ji,jj,1,1)-dl_lat1_ur) > dp_delta )THEN 2745 IF( tl_lat0%d_value(ji,jj,1,1) < dl_lat1_ur )THEN 2746 il_iur(2)=il_iur(2)+1 2747 IF( il_iur(2) > tl_lat0%t_dim(jp_J)%i_len )THEN 2748 CALL logger_error("GRID GET COARSE INDEX: error "//& 2749 & "computing upper right corner "//& 2750 & "index for latitude") 2751 ENDIF 2752 ENDIF 2753 ENDIF 2897 & dl_lon1, dl_lat1, 'ur' ) 2754 2898 2755 2899 ! coarse grid indices … … 2798 2942 ! 2799 2943 !> @author J.Paul 2800 !> - November, 2013- Initial Version2944 !> @date November, 2013 - Initial Version 2801 2945 ! 2802 2946 !> @param[in] td_lon longitude structure … … 2857 3001 END FUNCTION grid_is_global 2858 3002 !------------------------------------------------------------------- 2859 !> @brief This function return coarsegrid indices of the closest point2860 !> from fine gridpoint (lon1,lat1)3003 !> @brief This function return grid indices of the closest point 3004 !> from point (lon1,lat1) 2861 3005 !> 2862 3006 !> @details … … 2865 3009 !> of longitude and latitude, before running this function 2866 3010 !> 3011 !> if you add cd_pos argument, you could choice to return closest point at 3012 !> - lower left (ll) of the point 3013 !> - lower right (lr) of the point 3014 !> - upper left (ul) of the point 3015 !> - upper right (ur) of the point 3016 !> - lower (lo) of the point 3017 !> - upper (up) of the point 3018 !> - left (le) of the point 3019 !> - right (ri) of the point 3020 !> 2867 3021 !> @author J.Paul 2868 !> - November, 2013- Initial Version 3022 !> @date November, 2013 - Initial Version 3023 !> @date February, 2015 3024 !> - change dichotomy method to manage ORCA grid 3025 !> @date February, 2016 3026 !> - add optional use of relative position 2869 3027 ! 2870 3028 !> @param[in] dd_lon0 coarse grid array of longitude … … 2872 3030 !> @param[in] dd_lon1 fine grid longitude 2873 3031 !> @param[in] dd_lat1 fine grid latitude 3032 !> @param[in] cd_pos relative position of grid point from point 3033 !> @param[in] dd_fill fill value 2874 3034 !> @return coarse grid indices of closest point of fine grid point 2875 !> 2876 !------------------------------------------------------------------- 2877 FUNCTION grid_get_closest( dd_lon0, dd_lat0, dd_lon1, dd_lat1 ) 3035 !------------------------------------------------------------------- 3036 FUNCTION grid_get_closest( dd_lon0, dd_lat0, dd_lon1, dd_lat1, cd_pos, dd_fill ) 2878 3037 IMPLICIT NONE 2879 3038 ! Argument … … 2882 3041 REAL(dp), INTENT(IN) :: dd_lon1 2883 3042 REAL(dp), INTENT(IN) :: dd_lat1 3043 CHARACTER(LEN=*), INTENT(IN), OPTIONAL :: cd_pos 3044 REAL(dp), INTENT(IN), OPTIONAL :: dd_fill 2884 3045 2885 3046 ! function … … 2929 3090 2930 3091 ll_north=.FALSE. 2931 ll_continue=.TRUE. 2932 2933 ! look for meridian 0°/360° 2934 il_jmid = il_jinf + INT(il_shape(2)/2) 2935 il_ind(:) = MAXLOC( dl_lon0(:,il_jmid), dl_lon0(:,il_jmid) <= 360._dp ) 2936 2937 il_imid=il_ind(1) 2938 2939 IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. & 2940 & dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN 2941 2942 il_iinf = il_imid ; il_isup = il_imid 2943 il_jinf = il_jmid ; il_jsup = il_jmid 2944 2945 ll_continue=.FALSE. 2946 2947 ELSE 2948 IF( dl_lon1 < dl_lon0(il_isup,il_jmid) .AND. & 2949 & il_imid /= il_isup )THEN 2950 2951 ! point east 2952 il_iinf = il_imid 2953 2954 ELSE IF( dl_lon1 > dl_lon0(il_iinf,il_jmid) .AND. & 2955 & il_imid /= il_iinf )THEN 2956 2957 ! point west 2958 il_isup = il_imid 2959 2960 ENDIF 3092 ll_continue=.FALSE. 3093 3094 ! avoid to use fillvalue for reduce domain on first time 3095 IF( PRESENT(dd_fill) )THEN 3096 DO WHILE( ALL(dl_lon0(il_isup,:) == dd_fill) ) 3097 il_isup=il_isup-1 3098 ENDDO 3099 DO WHILE( ALL(dl_lon0(il_iinf,:) == dd_fill) ) 3100 il_iinf=il_iinf+1 3101 ENDDO 3102 DO WHILE( ALL(dd_lat0(:,il_jsup) == dd_fill) ) 3103 il_jsup=il_jsup-1 3104 ENDDO 3105 DO WHILE( ALL(dd_lat0(:,il_jinf) == dd_fill) ) 3106 il_jinf=il_jinf+1 3107 ENDDO 2961 3108 2962 3109 il_shape(1)= il_isup - il_iinf + 1 2963 3110 il_shape(2)= il_jsup - il_jinf + 1 2964 3111 2965 il_imid = il_iinf + INT(il_shape(1)/2) 3112 ENDIF 3113 3114 ! special case for north ORCA grid 3115 IF( dd_lat1 > 19. .AND. dl_lon1 < 74. )THEN 3116 ll_north=.TRUE. 3117 ENDIF 3118 3119 IF( .NOT. ll_north )THEN 3120 ! look for meridian 0°/360° 2966 3121 il_jmid = il_jinf + INT(il_shape(2)/2) 2967 2968 ! exit if too close from north fold (safer) 2969 IF( dd_lat0(il_imid,il_jmid) > 50.0 ) ll_north=.TRUE. 2970 2971 ! exit when close enough of point 2972 IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE. 3122 il_ind(:) = MAXLOC( dl_lon0(il_iinf:il_isup,il_jmid), & 3123 & dl_lon0(il_iinf:il_isup,il_jmid) <= 360._dp ) 3124 3125 il_imid=il_ind(1) 3126 3127 IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. & 3128 & dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN 3129 3130 il_iinf = il_imid ; il_isup = il_imid 3131 il_jinf = il_jmid ; il_jsup = il_jmid 3132 3133 ELSE 3134 IF( ALL(dl_lon0(il_isup,il_jinf:il_jsup) > dl_lon1 ) .AND. & 3135 & il_imid /= il_isup )THEN 3136 ! 0 < lon1 < lon0(isup) 3137 ! point east 3138 il_iinf = il_imid+1 3139 ll_continue=.TRUE. 3140 3141 ELSE IF( ALL(dl_lon0(il_iinf,il_jinf:il_jsup) < dl_lon1 ) .AND. & 3142 & il_imid /= il_iinf )THEN 3143 ! lon0(iinf) < lon1 < 360 3144 ! point west 3145 il_isup = il_imid 3146 ll_continue=.TRUE. 3147 3148 ENDIF 3149 3150 il_shape(1)= il_isup - il_iinf + 1 3151 il_shape(2)= il_jsup - il_jinf + 1 3152 3153 il_imid = il_iinf + INT(il_shape(1)/2) 3154 il_jmid = il_jinf + INT(il_shape(2)/2) 3155 3156 ! exit when close enough of point 3157 IF( ANY(il_shape(:) < 10 ) ) ll_continue=.FALSE. 3158 ENDIF 2973 3159 ENDIF 2974 3160 … … 2976 3162 DO WHILE( ll_continue .AND. .NOT. ll_north ) 2977 3163 3164 ll_continue=.FALSE. 2978 3165 IF( dl_lon1 == dl_lon0(il_imid,il_jmid) .AND. & 2979 3166 & dd_lat1 == dd_lat0(il_imid,il_jmid) )THEN … … 2982 3169 il_jinf = il_jmid ; il_jsup = il_jmid 2983 3170 2984 ll_continue=.FALSE.2985 2986 3171 ELSE 2987 IF( dl_lon1 > dl_lon0(il_imid,il_jmid) )THEN3172 IF( ALL(dl_lon0(il_imid,il_jinf:il_jsup) < dl_lon1) )THEN 2988 3173 2989 3174 ! point east 2990 3175 il_iinf = il_imid 3176 ll_continue=.TRUE. 2991 3177 2992 ELSE IF( dl_lon1 < dl_lon0(il_imid,il_jmid) )THEN3178 ELSE IF( ALL(dl_lon0(il_imid,il_jinf:il_jsup) > dl_lon1) )THEN 2993 3179 2994 3180 ! point west 2995 3181 il_isup = il_imid 3182 ll_continue=.TRUE. 2996 3183 2997 3184 ENDIF 2998 3185 2999 IF( dd_lat1 > dd_lat0(il_imid,il_jmid) )THEN3186 IF( ALL(dd_lat0(il_iinf:il_isup,il_jmid) < dd_lat1) )THEN 3000 3187 3001 3188 ! point north 3002 3189 il_jinf = il_jmid 3003 3004 ELSE IF(dd_lat1 < dd_lat0(il_imid,il_jmid) )THEN 3190 ll_continue=.TRUE. 3191 3192 ELSE IF( ALL(dd_lat0(il_iinf:il_isup,il_jmid) > dd_lat1) )THEN 3005 3193 3006 3194 ! point south 3007 3195 il_jsup = il_jmid 3196 ll_continue=.TRUE. 3008 3197 3009 3198 ENDIF … … 3014 3203 il_imid = il_iinf + INT(il_shape(1)/2) 3015 3204 il_jmid = il_jinf + INT(il_shape(2)/2) 3016 3017 ! exit if too close from north fold (safer)3018 IF( dd_lat0(il_imid,il_jmid) > 50.0 ) ll_north=.TRUE.3019 3205 3020 3206 ! exit when close enough of point … … 3034 3220 & dl_lon1, dd_lat1 ) 3035 3221 3222 IF( PRESENT(cd_pos) )THEN 3223 ! 3224 SELECT CASE(TRIM(cd_pos)) 3225 CASE('le') 3226 WHERE( dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 ) 3227 dl_dist(:,:)=NF90_FILL_DOUBLE 3228 END WHERE 3229 CASE('ri') 3230 WHERE( dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 ) 3231 dl_dist(:,:)=NF90_FILL_DOUBLE 3232 END WHERE 3233 CASE('up') 3234 WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 ) 3235 dl_dist(:,:)=NF90_FILL_DOUBLE 3236 END WHERE 3237 CASE('lo') 3238 WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 ) 3239 dl_dist(:,:)=NF90_FILL_DOUBLE 3240 END WHERE 3241 CASE('ll') 3242 WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 .OR. & 3243 & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 ) 3244 dl_dist(:,:)=NF90_FILL_DOUBLE 3245 END WHERE 3246 CASE('lr') 3247 WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 .OR. & 3248 & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) > dd_lat1 ) 3249 dl_dist(:,:)=NF90_FILL_DOUBLE 3250 END WHERE 3251 CASE('ul') 3252 WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) > dl_lon1 .OR. & 3253 & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 ) 3254 dl_dist(:,:)=NF90_FILL_DOUBLE 3255 END WHERE 3256 CASE('ur') 3257 WHERE( dl_lon0(il_iinf:il_isup,il_jinf:il_jsup) < dl_lon1 .OR. & 3258 & dd_lat0(il_iinf:il_isup,il_jinf:il_jsup) < dd_lat1 ) 3259 dl_dist(:,:)=NF90_FILL_DOUBLE 3260 END WHERE 3261 END SELECT 3262 ENDIF 3036 3263 grid_get_closest(:)=MINLOC(dl_dist(:,:),dl_dist(:,:)/=NF90_FILL_DOUBLE) 3037 3264 … … 3049 3276 ! 3050 3277 !> @author J.Paul 3051 !> - November, 2013- Initial Version3278 !> @date November, 2013 - Initial Version 3052 3279 ! 3053 3280 !> @param[in] dd_lon grid longitude array … … 3055 3282 !> @param[in] dd_lonA longitude of point A 3056 3283 !> @param[in] dd_latA latitude of point A 3284 !> @param[in] dd_fill 3057 3285 !> @return array of distance between point A and grid points. 3058 3286 !------------------------------------------------------------------- 3059 FUNCTION grid_distance(dd_lon, dd_lat, dd_lonA, dd_latA )3287 FUNCTION grid_distance(dd_lon, dd_lat, dd_lonA, dd_latA ) 3060 3288 IMPLICIT NONE 3061 3289 ! Argument … … 3110 3338 DO ji=1,il_shape(1) 3111 3339 IF( dl_lon(ji,jj) == dl_lonA .AND. & 3112 & dl_lat(ji,jj) == dl_la TA )THEN3340 & dl_lat(ji,jj) == dl_latA )THEN 3113 3341 grid_distance(ji,jj)=0.0 3114 3342 ELSE 3115 3343 dl_tmp= SIN(dl_latA)*SIN(dl_lat(ji,jj)) + & 3116 & COS(dl_latA)*COS(dl_lat(ji,jj))*COS(dl_lon(ji,jj)-dl_lonA) 3344 & COS(dl_latA)*COS(dl_lat(ji,jj)) * & 3345 & COS(dl_lon(ji,jj)-dl_lonA) 3117 3346 ! check to avoid mistake with ACOS 3118 3347 IF( dl_tmp < -1.0 ) dl_tmp = -1.0 … … 3136 3365 ! 3137 3366 !> @author J.Paul 3138 !> - September, 2014- Initial Version3367 !> @date September, 2014 - Initial Version 3139 3368 !> @date October, 2014 3140 3369 !> - work on mpp file structure instead of file structure … … 3170 3399 3171 3400 ! local variable 3172 INTEGER(i4) :: il_imin0 3173 INTEGER(i4) :: il_jmin0 3174 INTEGER(i4) :: il_imax0 3175 INTEGER(i4) :: il_jmax0 3401 INTEGER(i4) :: il_imin0 3402 INTEGER(i4) :: il_jmin0 3403 INTEGER(i4) :: il_imax0 3404 INTEGER(i4) :: il_jmax0 3405 INTEGER(i4) :: il_ind 3176 3406 3177 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho3407 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 3178 3408 3179 INTEGER(i4), DIMENSION(2,2) :: il_xghost03180 INTEGER(i4), DIMENSION(2,2) :: il_xghost13181 3182 CHARACTER(LEN= 1) :: cl_point3183 CHARACTER(LEN=lc) :: cl_name3409 INTEGER(i4), DIMENSION(2,2) :: il_xghost0 3410 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 3411 3412 CHARACTER(LEN= 1) :: cl_point 3413 CHARACTER(LEN=lc) :: cl_name 3184 3414 3185 3415 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon0 … … 3188 3418 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lat1 3189 3419 3190 TYPE(TVAR) :: tl_lon03191 TYPE(TVAR) :: tl_lat03192 TYPE(TVAR) :: tl_lon13193 TYPE(TVAR) :: tl_lat13194 3195 TYPE(TMPP) :: tl_coord03196 TYPE(TMPP) :: tl_coord13420 TYPE(TVAR) :: tl_lon0 3421 TYPE(TVAR) :: tl_lat0 3422 TYPE(TVAR) :: tl_lon1 3423 TYPE(TVAR) :: tl_lat1 3424 3425 TYPE(TMPP) :: tl_coord0 3426 TYPE(TMPP) :: tl_coord1 3197 3427 3198 3428 ! loop indices … … 3227 3457 ! read coarse longitue and latitude 3228 3458 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 3459 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 3460 IF( il_ind == 0 )THEN 3461 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3462 & TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". & 3463 & try to use longitude.") 3464 WRITE(cl_name,*) 'longitude' 3465 ENDIF 3229 3466 tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 3467 3230 3468 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 3469 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 3470 IF( il_ind == 0 )THEN 3471 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3472 & TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". & 3473 & try to use latitude.") 3474 WRITE(cl_name,*) 'latitude' 3475 ENDIF 3231 3476 tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 3232 3477 … … 3267 3512 ! read fine longitue and latitude 3268 3513 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 3514 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 3515 IF( il_ind == 0 )THEN 3516 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3517 & TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". & 3518 & try to use longitude.") 3519 WRITE(cl_name,*) 'longitude' 3520 ENDIF 3269 3521 tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 3522 3270 3523 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 3524 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 3525 IF( il_ind == 0 )THEN 3526 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3527 & TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". & 3528 & try to use latitude.") 3529 WRITE(cl_name,*) 'latitude' 3530 ENDIF 3271 3531 tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 3272 3532 … … 3297 3557 & il_imax0, il_jmax0, & 3298 3558 & dl_lon1(:,:), dl_lat1(:,:),& 3299 & id_rho(:) )3559 & id_rho(:), cl_point ) 3300 3560 3301 3561 DEALLOCATE(dl_lon0, dl_lat0) … … 3318 3578 ! 3319 3579 !> @author J.Paul 3320 !> - September, 2014- Initial Version3580 !> @date September, 2014 - Initial Version 3321 3581 !> @date October, 2014 3322 3582 !> - work on mpp file structure instead of file structure … … 3354 3614 3355 3615 ! local variable 3356 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 3616 INTEGER(i4) :: il_ind 3617 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 3618 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 3357 3619 3358 INTEGER(i4), DIMENSION(2,2) :: il_xghost1 3359 3360 CHARACTER(LEN= 1) :: cl_point 3361 CHARACTER(LEN=lc) :: cl_name 3620 CHARACTER(LEN= 1) :: cl_point 3621 CHARACTER(LEN=lc) :: cl_name 3362 3622 3363 3623 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon1 3364 3624 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lat1 3365 3625 3366 TYPE(TVAR) :: tl_lon13367 TYPE(TVAR) :: tl_lat13368 3369 TYPE(TMPP) :: tl_coord13626 TYPE(TVAR) :: tl_lon1 3627 TYPE(TVAR) :: tl_lat1 3628 3629 TYPE(TMPP) :: tl_coord1 3370 3630 ! loop indices 3371 3631 !---------------------------------------------------------------- … … 3397 3657 ! read fine longitue and latitude 3398 3658 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 3659 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 3660 IF( il_ind == 0 )THEN 3661 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3662 & TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". & 3663 & try to use longitude.") 3664 WRITE(cl_name,*) 'longitude' 3665 ENDIF 3399 3666 tl_lon1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 3667 3400 3668 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 3669 il_ind=var_get_id(tl_coord1%t_proc(1)%t_var(:), cl_name) 3670 IF( il_ind == 0 )THEN 3671 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3672 & TRIM(cl_name)//" in file "//TRIM(tl_coord1%c_name)//". & 3673 & try to use latitude.") 3674 WRITE(cl_name,*) 'latitude' 3675 ENDIF 3401 3676 tl_lat1=iom_mpp_read_var(tl_coord1, TRIM(cl_name)) 3402 3677 … … 3427 3702 & id_imax0, id_jmax0, & 3428 3703 & dl_lon1(:,:), dl_lat1(:,:),& 3429 & id_rho(:) )3704 & id_rho(:), cl_point ) 3430 3705 3431 3706 DEALLOCATE(dl_lon1, dl_lat1) … … 3446 3721 ! 3447 3722 !> @author J.Paul 3448 !> - September, 2014- Initial Version3723 !> @date September, 2014 - Initial Version 3449 3724 !> @date October, 2014 3450 3725 !> - work on mpp file structure instead of file structure … … 3483 3758 3484 3759 ! local variable 3485 INTEGER(i4) :: il_imin0 3486 INTEGER(i4) :: il_jmin0 3487 INTEGER(i4) :: il_imax0 3488 INTEGER(i4) :: il_jmax0 3760 INTEGER(i4) :: il_imin0 3761 INTEGER(i4) :: il_jmin0 3762 INTEGER(i4) :: il_imax0 3763 INTEGER(i4) :: il_jmax0 3764 INTEGER(i4) :: il_ind 3489 3765 3490 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho3766 INTEGER(i4), DIMENSION(:), ALLOCATABLE :: il_rho 3491 3767 3492 INTEGER(i4), DIMENSION(2,2) :: il_xghost03493 3494 CHARACTER(LEN= 1) :: cl_point3495 CHARACTER(LEN=lc) :: cl_name3768 INTEGER(i4), DIMENSION(2,2) :: il_xghost0 3769 3770 CHARACTER(LEN= 1) :: cl_point 3771 CHARACTER(LEN=lc) :: cl_name 3496 3772 3497 3773 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon0 3498 3774 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lat0 3499 3775 3500 TYPE(TVAR) :: tl_lon03501 TYPE(TVAR) :: tl_lat03502 3503 TYPE(TMPP) :: tl_coord03776 TYPE(TVAR) :: tl_lon0 3777 TYPE(TVAR) :: tl_lat0 3778 3779 TYPE(TMPP) :: tl_coord0 3504 3780 ! loop indices 3505 3781 !---------------------------------------------------------------- 3506 3782 ! init 3507 3783 grid__get_fine_offset_fc(:,:)=-1 3508 3509 3784 ALLOCATE(il_rho(ip_maxdim)) 3510 3785 il_rho(:)=1 … … 3528 3803 CALL iom_mpp_open(tl_coord0) 3529 3804 3530 ! read coarse longitu e and latitude3805 ! read coarse longitude and latitude 3531 3806 WRITE(cl_name,*) 'longitude_'//TRIM(cl_point) 3807 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 3808 IF( il_ind == 0 )THEN 3809 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3810 & TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". & 3811 & try to use longitude.") 3812 WRITE(cl_name,*) 'longitude' 3813 ENDIF 3532 3814 tl_lon0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 3815 3533 3816 WRITE(cl_name,*) 'latitude_'//TRIM(cl_point) 3817 il_ind=var_get_id(tl_coord0%t_proc(1)%t_var(:), cl_name) 3818 IF( il_ind == 0 )THEN 3819 CALL logger_warn("GRID GET FINE OFFSET: no variable "//& 3820 & TRIM(cl_name)//" in file "//TRIM(tl_coord0%c_name)//". & 3821 & try to use latitude.") 3822 WRITE(cl_name,*) 'latitude' 3823 ENDIF 3534 3824 tl_lat0=iom_mpp_read_var(tl_coord0, TRIM(cl_name)) 3535 3825 3536 3826 ! close mpp files 3537 3827 CALL iom_mpp_close(tl_coord0) … … 3539 3829 CALL grid_del_ghost(tl_lon0, il_xghost0(:,:)) 3540 3830 CALL grid_del_ghost(tl_lat0, il_xghost0(:,:)) 3831 3541 3832 3542 3833 ALLOCATE(dl_lon0(tl_lon0%t_dim(jp_I)%i_len, & … … 3561 3852 il_jmax0=id_jmax0-il_xghost0(jp_J,1) 3562 3853 3563 3564 3854 !3- compute 3565 3855 grid__get_fine_offset_fc(:,:)=grid_get_fine_offset(& … … 3568 3858 & il_imax0, il_jmax0, & 3569 3859 & dd_lon1(:,:), dd_lat1(:,:),& 3570 & id_rho(:) )3860 & id_rho(:), cl_point ) 3571 3861 3572 3862 DEALLOCATE(dl_lon0, dl_lat0) … … 3585 3875 ! 3586 3876 !> @author J.Paul 3587 !> - November, 2013 - Initial Version 3588 !> @date September, 2014 - rename from grid_get_fine_offset 3589 ! 3877 !> @date November, 2013 - Initial Version 3878 !> @date September, 2014 3879 !> - rename from grid_get_fine_offset 3880 !> @date May, 2015 3881 !> - improve way to find offset 3882 !> @date July, 2015 3883 !> - manage case close to greenwich meridian 3884 !> @date February, 2016 3885 !> - use grid_get_closest to assess offset 3886 !> - use delta (lon or lat) 3887 !> - manage cases for T,U,V or F point, with even or odd refinment 3888 !> - check lower left(upper right) fine grid point inside lower left(upper 3889 !> right) coarse grid cell. 3890 !> 3891 !> @todo check case close from North fold. 3892 !> 3590 3893 !> @param[in] dd_lon0 coarse grid longitude array 3591 3894 !> @param[in] dd_lat0 coarse grid latitude array … … 3597 3900 !> @param[in] dd_lat1 fine grid latitude array 3598 3901 !> @param[in] id_rho array of refinement factor 3902 !> @param[in] cd_point Arakawa grid point 3599 3903 !> @return offset array (/ (/i_offset_left,i_offset_right/),(/j_offset_lower,j_offset_upper/) /) 3600 3904 !------------------------------------------------------------------- 3601 3905 FUNCTION grid__get_fine_offset_cc( dd_lon0, dd_lat0, & 3602 3906 & id_imin0, id_jmin0, id_imax0, id_jmax0, & 3603 & dd_lon1, dd_lat1, id_rho )3907 & dd_lon1, dd_lat1, id_rho, cd_point ) 3604 3908 IMPLICIT NONE 3605 3909 ! Argument 3606 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lon0 3607 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lat0 3608 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lon1 3609 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lat1 3610 3611 INTEGER(i4), INTENT(IN) :: id_imin0 3612 INTEGER(i4), INTENT(IN) :: id_jmin0 3613 INTEGER(i4), INTENT(IN) :: id_imax0 3614 INTEGER(i4), INTENT(IN) :: id_jmax0 3615 3616 INTEGER(i4), DIMENSION(:) , INTENT(IN) :: id_rho 3910 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lon0 3911 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lat0 3912 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lon1 3913 REAL(dp) , DIMENSION(:,:), INTENT(IN) :: dd_lat1 3914 3915 INTEGER(i4) , INTENT(IN) :: id_imin0 3916 INTEGER(i4) , INTENT(IN) :: id_jmin0 3917 INTEGER(i4) , INTENT(IN) :: id_imax0 3918 INTEGER(i4) , INTENT(IN) :: id_jmax0 3919 3920 INTEGER(i4) , DIMENSION(:) , INTENT(IN) :: id_rho 3921 CHARACTER(LEN=*) , INTENT(IN), OPTIONAL :: cd_point 3617 3922 3618 3923 ! function … … 3620 3925 3621 3926 ! local variable 3622 INTEGER(i4), DIMENSION(2) :: il_shape0 3623 INTEGER(i4), DIMENSION(2) :: il_shape1 3927 CHARACTER(LEN= 1) :: cl_point 3928 3929 INTEGER(i4) :: i1 3930 INTEGER(i4) :: i2 3931 INTEGER(i4) :: j1 3932 INTEGER(i4) :: j2 3933 3934 INTEGER(i4), DIMENSION(2) :: il_shape0 3935 INTEGER(i4), DIMENSION(2) :: il_shape1 3936 3937 INTEGER(i4), DIMENSION(2) :: il_ind 3938 3624 3939 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon0 3625 3940 REAL(dp) , DIMENSION(:,:), ALLOCATABLE :: dl_lon1 3941 3942 REAL(dp) :: dl_lonmax0 3943 REAL(dp) :: dl_latmax0 3944 REAL(dp) :: dl_lonmin0 3945 REAL(dp) :: dl_latmin0 3946 3947 REAL(dp) :: dl_lon0F 3948 REAL(dp) :: dl_lat0F 3949 REAL(dp) :: dl_dlon 3950 REAL(dp) :: dl_dlat 3951 3952 LOGICAL , DIMENSION(2) :: ll_even 3953 LOGICAL :: ll_greenwich 3626 3954 3627 3955 ! loop indices 3628 INTEGER(i4) :: ji3629 INTEGER(i4) :: jj3630 3631 3956 INTEGER(i4) :: ii 3632 3957 INTEGER(i4) :: ij … … 3640 3965 CALL logger_fatal("GRID GET FINE OFFSET: dimension of fine "//& 3641 3966 & "longitude and latitude differ") 3642 ENDIF 3967 ENDIF 3968 3969 ll_even(:)=(/ (MOD(id_rho(jp_I),2)==0), (MOD(id_rho(jp_J),2)==0) /) 3970 3971 cl_point='T' 3972 IF( PRESENT(cd_point) ) cl_point=TRIM(fct_upper(cd_point)) 3643 3973 3644 3974 il_shape0(:)=SHAPE(dd_lon0(:,:)) 3645 3975 ALLOCATE( dl_lon0(il_shape0(1),il_shape0(2)) ) 3646 3976 3977 il_shape1(:)=SHAPE(dd_lon1(:,:)) 3978 ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) ) 3979 3647 3980 dl_lon0(:,:)=dd_lon0(:,:) 3648 3981 WHERE( dd_lon0(:,:) < 0 ) dl_lon0(:,:)=dd_lon0(:,:)+360. 3649 3982 3650 il_shape1(:)=SHAPE(dd_lon1(:,:))3651 ALLOCATE( dl_lon1(il_shape1(1),il_shape1(2)) )3652 3653 3983 dl_lon1(:,:)=dd_lon1(:,:) 3654 WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360. 3984 WHERE( dd_lon1(:,:) < 0 ) dl_lon1(:,:)=dd_lon1(:,:)+360. 3655 3985 3656 3986 ! init 3657 3987 grid__get_fine_offset_cc(:,:)=-1 3658 3659 IF( il_shape1(1) > 1 )THEN 3660 3661 ! look for i-direction left offset 3662 IF( dl_lon1(1,1) < dl_lon0(id_imin0+1,id_jmin0) )THEN 3663 DO ji=1,id_rho(jp_I)+2 3664 IF( dl_lon1(ji,1) > dl_lon0(id_imin0+1,id_jmin0) - dp_delta )THEN 3665 grid__get_fine_offset_cc(1,1)=(id_rho(jp_I)+1)-ji 3666 EXIT 3667 ENDIF 3668 ENDDO 3988 ll_greenwich=.FALSE. 3989 3990 IF( il_shape1(jp_J) == 1 )THEN 3991 3992 grid__get_fine_offset_cc(jp_J,:)=((id_rho(jp_J)-1)/2) 3993 3994 !!! work on i-direction 3995 !!! look for i-direction left offset 3996 i1=1 ; i2=MIN((id_rho(jp_I)+2),il_shape1(jp_I)) 3997 j1=1 ; j2=1 3998 3999 ! check if cross greenwich meridien 4000 IF( minval(dl_lon0(id_imin0:id_imin0+1,id_jmin0))<5. .OR. & 4001 & maxval(dl_lon0(id_imin0:id_imin0+1,id_jmin0))>355. )THEN 4002 ! close to greenwich meridien 4003 ll_greenwich=.TRUE. 4004 ! 0:360 => -180:180 4005 WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0) > 180. ) 4006 dl_lon0(id_imin0:id_imin0+1,id_jmin0) = & 4007 & dl_lon0(id_imin0:id_imin0+1,id_jmin0)-360. 4008 END WHERE 4009 4010 WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 4011 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 4012 END WHERE 4013 ENDIF 4014 4015 ! max lognitude of the left cell 4016 dl_lonmax0=dl_lon0(id_imin0+1,id_jmin0) 4017 IF( dl_lon1(1,1) < dl_lonmax0 )THEN 4018 4019 !!!!! i-direction !!!!! 4020 IF( ll_even(jp_I) )THEN 4021 ! even 4022 SELECT CASE(TRIM(cl_point)) 4023 CASE('F','U') 4024 dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0) - & 4025 & dl_lon0(id_imin0 ,id_jmin0) ) / & 4026 & ( 2.*id_rho(jp_I) ) 4027 CASE DEFAULT 4028 dl_dlon=0 4029 END SELECT 4030 ELSE 4031 ! odd 4032 dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0) - & 4033 & dl_lon0(id_imin0 ,id_jmin0) ) / & 4034 & ( 2.*id_rho(jp_I) ) 4035 ENDIF 4036 4037 dl_lon0F= dl_lon0(id_imin0+1,id_jmin0) + dl_dlon 4038 dl_lat0F= dd_lat0(id_imin0+1,id_jmin0) 4039 4040 il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 4041 & dl_lon0F, dl_lat0F, 'le' ) 4042 4043 ii=il_ind(1) 4044 4045 !!!!! i-direction !!!!! 4046 IF( ll_even(jp_I) )THEN 4047 ! even 4048 SELECT CASE(TRIM(cl_point)) 4049 CASE('T','V') 4050 grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii 4051 CASE DEFAULT !'F','U' 4052 grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 4053 END SELECT 4054 ELSE 4055 ! odd 4056 grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 4057 ENDIF 4058 3669 4059 ELSE 3670 4060 CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 3671 & " not match fine grid lower left corner.") 3672 ENDIF 3673 3674 ! look for i-direction right offset 3675 IF( dl_lon1(il_shape1(1),1) > dl_lon0(id_imax0-1,id_jmin0) )THEN 3676 DO ji=1,id_rho(jp_I)+2 3677 ii=il_shape1(1)-ji+1 3678 IF( dl_lon1(ii,1) < dl_lon0(id_imax0-1,id_jmin0) + dp_delta )THEN 3679 grid__get_fine_offset_cc(1,2)=(id_rho(jp_I)+1)-ji 3680 EXIT 3681 ENDIF 3682 ENDDO 4061 & " not match fine grid left corner.") 4062 ENDIF 4063 4064 IF( ll_greenwich )THEN 4065 ! close to greenwich meridien 4066 ll_greenwich=.FALSE. 4067 ! -180:180 => 0:360 4068 WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0) < 0. ) 4069 dl_lon0(id_imin0:id_imin0+1,id_jmin0) = & 4070 & dl_lon0(id_imin0:id_imin0+1,id_jmin0)+360. 4071 END WHERE 4072 4073 WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 4074 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 4075 END WHERE 4076 ENDIF 4077 4078 !!!!!! look for i-direction right offset !!!!!! 4079 i1=MAX(1,il_shape1(jp_I)-(id_rho(jp_I)+2)+1) ; i2=il_shape1(jp_I) 4080 j1=1 ; j2=1 4081 4082 ! check if cross greenwich meridien 4083 IF( minval(dl_lon0(id_imax0-1:id_imax0,id_jmin0))<5. .OR. & 4084 & maxval(dl_lon0(id_imax0-1:id_imax0,id_jmin0))>355. )THEN 4085 ! close to greenwich meridien 4086 ll_greenwich=.TRUE. 4087 ! 0:360 => -180:180 4088 WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmin0) > 180. ) 4089 dl_lon0(id_imax0-1:id_imax0,id_jmin0) = & 4090 & dl_lon0(id_imax0-1:id_imax0,id_jmin0)-360. 4091 END WHERE 4092 4093 WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 4094 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 4095 END WHERE 4096 ENDIF 4097 4098 ! min lognitude of the right cell 4099 dl_lonmin0=dl_lon0(id_imax0-1,id_jmin0) 4100 IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > dl_lonmin0 )THEN 4101 4102 !!!!! i-direction !!!!! 4103 IF( ll_even(jp_I) )THEN 4104 ! even 4105 SELECT CASE(TRIM(cl_point)) 4106 CASE('F','U') 4107 dl_dlon= ( dl_lon0(id_imax0 ,id_jmin0) - & 4108 & dl_lon0(id_imax0-1,id_jmin0) ) / & 4109 & ( 2.*id_rho(jp_I) ) 4110 CASE DEFAULT 4111 dl_dlon=0 4112 END SELECT 4113 ELSE 4114 ! odd 4115 dl_dlon= ( dl_lon0(id_imax0 ,id_jmin0) - & 4116 & dl_lon0(id_imax0-1,id_jmin0) ) / & 4117 & ( 2.*id_rho(jp_I) ) 4118 ENDIF 4119 4120 dl_lon0F= dl_lon0(id_imax0-1,id_jmin0) - dl_dlon 4121 dl_lat0F= dd_lat0(id_imax0-1,id_jmin0) 4122 4123 il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 4124 & dl_lon0F, dl_lat0F, 'ri' ) 4125 4126 ii=(MIN(il_shape1(jp_I),(id_rho(jp_I)+2))-il_ind(1)+1) 4127 4128 !!!!! i-direction !!!!! 4129 IF( ll_even(jp_I) )THEN 4130 ! even 4131 SELECT CASE(TRIM(cl_point)) 4132 CASE('T','V') 4133 grid__get_fine_offset_cc(jp_I,2)=id_rho(jp_I)-ii 4134 CASE DEFAULT !'F','U' 4135 grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 4136 END SELECT 4137 ELSE 4138 ! odd 4139 grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 4140 ENDIF 4141 3683 4142 ELSE 3684 4143 CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 3685 & " not match fine grid lower right corner.") 3686 ENDIF 3687 3688 ELSE 3689 grid__get_fine_offset_cc(1,:)=((id_rho(jp_I)-1)/2) 3690 ENDIF 3691 3692 IF( il_shape1(2) > 1 )THEN 3693 3694 ! look for j-direction lower offset 3695 IF( dd_lat1(1,1) < dd_lat0(id_imin0,id_jmin0+1) )THEN 3696 DO jj=1,id_rho(jp_J)+2 3697 IF( dd_lat1(1,jj) > dd_lat0(id_imin0,id_jmin0+1) - dp_delta )THEN 3698 grid__get_fine_offset_cc(2,1)=(id_rho(jp_J)+1)-jj 3699 EXIT 3700 ENDIF 3701 ENDDO 4144 & " not match fine grid right corner.") 4145 ENDIF 4146 4147 IF( ll_greenwich )THEN 4148 ! close to greenwich meridien 4149 ll_greenwich=.FALSE. 4150 ! -180:180 => 0:360 4151 WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmin0) < 0. ) 4152 dl_lon0(id_imax0-1:id_imax0,id_jmin0) = & 4153 & dl_lon0(id_imax0-1:id_imax0,id_jmin0)+360. 4154 END WHERE 4155 4156 WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 4157 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 4158 END WHERE 4159 ENDIF 4160 4161 ELSEIF( il_shape1(jp_I) == 1 )THEN 4162 4163 grid__get_fine_offset_cc(jp_I,:)=((id_rho(jp_I)-1)/2) 4164 4165 !!! work on j-direction 4166 !!! look for j-direction lower offset 4167 i1=1 ; i2=1 4168 j1=1 ; j2=MIN((id_rho(jp_J)+2),il_shape1(jp_J)) 4169 4170 4171 ! max latitude of the lower cell 4172 dl_latmax0=dd_lat0(id_imin0,id_jmin0+1) 4173 IF( dd_lat1(1,1) < dl_latmax0 )THEN 4174 4175 IF( ll_even(jp_J) )THEN 4176 ! even 4177 SELECT CASE(TRIM(cl_point)) 4178 CASE('F','V') 4179 dl_dlat= ( dd_lat0(id_imin0,id_jmin0+1) - & 4180 & dd_lat0(id_imin0,id_jmin0 ) ) / & 4181 & ( 2.*id_rho(jp_J) ) 4182 CASE DEFAULT 4183 dl_dlat=0 4184 END SELECT 4185 ELSE 4186 ! odd 4187 dl_dlat= ( dd_lat0(id_imin0,id_jmin0+1) - & 4188 & dd_lat0(id_imin0,id_jmin0 ) ) / & 4189 & ( 2.*id_rho(jp_J) ) 4190 ENDIF 4191 4192 dl_lon0F= dl_lon0(id_imin0,id_jmin0+1) 4193 dl_lat0F= dd_lat0(id_imin0,id_jmin0+1) + dl_dlat 4194 4195 il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 4196 & dl_lon0F, dl_lat0F, 'lo' ) 4197 4198 ij=il_ind(2) 4199 4200 !!!!! i-direction !!!!! 4201 IF( ll_even(jp_I) )THEN 4202 ! even 4203 SELECT CASE(TRIM(cl_point)) 4204 CASE('T','V') 4205 grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii 4206 CASE DEFAULT !'F','U' 4207 grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 4208 END SELECT 4209 ELSE 4210 ! odd 4211 grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 4212 ENDIF 4213 3702 4214 ELSE 3703 4215 CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 3704 & " not match fine grid upper left corner.") 3705 ENDIF 3706 3707 ! look for j-direction upper offset 3708 IF( dd_lat1(1,il_shape1(2)) > dd_lat0(id_imin0,id_jmax0-1) )THEN 3709 DO jj=1,id_rho(jp_J)+2 3710 ij=il_shape1(2)-jj+1 3711 IF( dd_lat1(1,ij) < dd_lat0(id_imin0,id_jmax0-1) + dp_delta )THEN 3712 grid__get_fine_offset_cc(2,2)=(id_rho(jp_J)+1)-jj 3713 EXIT 3714 ENDIF 3715 ENDDO 4216 & " not match fine grid lower corner.") 4217 ENDIF 4218 4219 !!! look for j-direction upper offset 4220 i1=1 ; i2=1 4221 j1=MAX(1,il_shape1(jp_J)-(id_rho(jp_J)+2)+1) ; j2=il_shape1(jp_J) 4222 4223 ! min latitude of the upper cell 4224 dl_latmin0=dd_lat0(id_imin0,id_jmax0-1) 4225 IF( dd_lat1(il_shape1(jp_I),il_shape1(jp_J)) > dl_latmin0 )THEN 4226 4227 IF( ll_even(jp_J) )THEN 4228 ! even 4229 SELECT CASE(TRIM(cl_point)) 4230 CASE('F','V') 4231 dl_dlat= ( dd_lat0(id_imin0,id_jmax0 ) - & 4232 & dd_lat0(id_imin0,id_jmax0-1) ) / & 4233 & ( 2.*id_rho(jp_J) ) 4234 CASE DEFAULT 4235 dl_dlat=0 4236 END SELECT 4237 ELSE 4238 ! odd 4239 dl_dlat= ( dd_lat0(id_imin0,id_jmax0 ) - & 4240 & dd_lat0(id_imin0,id_jmax0-1) ) / & 4241 & ( 2*id_rho(jp_J) ) 4242 ENDIF 4243 4244 dl_lon0F= dl_lon0(id_imin0,id_jmax0-1) 4245 dl_lat0F= dd_lat0(id_imin0,id_jmax0-1) - dl_dlat 4246 4247 il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 4248 & dl_lon0F, dl_lat0F, 'up' ) 4249 4250 ij=(MIN(il_shape1(jp_J),(id_rho(jp_J)+2))-il_ind(2)+1) 4251 4252 !!!!! j-direction !!!!! 4253 IF( ll_even(jp_J) )THEN 4254 ! even 4255 SELECT CASE(TRIM(cl_point)) 4256 CASE('T','U') 4257 grid__get_fine_offset_cc(jp_J,2)=id_rho(jp_J)-ij 4258 CASE DEFAULT !'F','V' 4259 grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 4260 END SELECT 4261 ELSE 4262 ! odd 4263 grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 4264 ENDIF 4265 3716 4266 ELSE 3717 4267 CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do "//& 4268 & " not match fine grid upper corner.") 4269 ENDIF 4270 4271 ELSE ! il_shape1(1) > 1 .AND. il_shape1(2) > 1 4272 4273 !!!!!! look for lower left offset !!!!!! 4274 i1=1 ; i2=MIN((id_rho(jp_I)+2),il_shape1(jp_I)) 4275 j1=1 ; j2=MIN((id_rho(jp_J)+2),il_shape1(jp_J)) 4276 4277 ! check if cross greenwich meridien 4278 IF( minval(dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1))<5. .OR. & 4279 & maxval(dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1))>355. )THEN 4280 ! close to greenwich meridien 4281 ll_greenwich=.TRUE. 4282 ! 0:360 => -180:180 4283 WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) > 180. ) 4284 dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) = & 4285 & dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1)-360. 4286 END WHERE 4287 4288 WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 4289 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 4290 END WHERE 4291 ENDIF 4292 4293 ! max longitude of the lower left cell 4294 dl_lonmax0=MAX(dl_lon0(id_imin0+1,id_jmin0),dl_lon0(id_imin0+1,id_jmin0+1)) 4295 ! max latitude of the lower left cell 4296 dl_latmax0=MAX(dd_lat0(id_imin0,id_jmin0+1),dd_lat0(id_imin0+1,id_jmin0+1)) 4297 IF( dl_lon1(1,1) < dl_lonmax0 .AND. & 4298 & dd_lat1(1,1) < dl_latmax0 )THEN 4299 4300 !!!!! i-direction !!!!! 4301 IF( ll_even(jp_I) )THEN 4302 ! even 4303 SELECT CASE(TRIM(cl_point)) 4304 CASE('F','U') 4305 dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0+1) - & 4306 & dl_lon0(id_imin0 ,id_jmin0+1) ) / & 4307 & ( 2.*id_rho(jp_I) ) 4308 CASE DEFAULT 4309 dl_dlon=0 4310 END SELECT 4311 ELSE 4312 ! odd 4313 dl_dlon= ( dl_lon0(id_imin0+1,id_jmin0+1) - & 4314 & dl_lon0(id_imin0 ,id_jmin0+1) ) / & 4315 & ( 2.*id_rho(jp_I) ) 4316 ENDIF 4317 4318 !!!!! j-direction !!!!! 4319 IF( ll_even(jp_J) )THEN 4320 ! even 4321 SELECT CASE(TRIM(cl_point)) 4322 CASE('F','V') 4323 dl_dlat= ( dd_lat0(id_imin0+1,id_jmin0+1) - & 4324 & dd_lat0(id_imin0+1,id_jmin0 ) ) / & 4325 & ( 2.*id_rho(jp_J) ) 4326 CASE DEFAULT 4327 dl_dlat=0 4328 END SELECT 4329 ELSE 4330 ! odd 4331 dl_dlat= ( dd_lat0(id_imin0+1,id_jmin0+1) - & 4332 & dd_lat0(id_imin0+1,id_jmin0 ) ) / & 4333 & ( 2.*id_rho(jp_J) ) 4334 ENDIF 4335 4336 dl_lon0F= dl_lon0(id_imin0+1,id_jmin0+1) + dl_dlon 4337 dl_lat0F= dd_lat0(id_imin0+1,id_jmin0+1) + dl_dlat 4338 4339 il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 4340 & dl_lon0F, dl_lat0F, 'll' ) 4341 4342 ii=il_ind(1) 4343 ij=il_ind(2) 4344 4345 !!!!! i-direction !!!!! 4346 IF( ll_even(jp_I) )THEN 4347 ! even 4348 SELECT CASE(TRIM(cl_point)) 4349 CASE('T','V') 4350 grid__get_fine_offset_cc(jp_I,1)=id_rho(jp_I)-ii 4351 CASE DEFAULT !'F','U' 4352 grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 4353 END SELECT 4354 ELSE 4355 ! odd 4356 grid__get_fine_offset_cc(jp_I,1)=(id_rho(jp_I)+1)-ii 4357 ENDIF 4358 4359 !!!!! j-direction !!!!! 4360 IF( ll_even(jp_J) )THEN 4361 ! even 4362 SELECT CASE(TRIM(cl_point)) 4363 CASE('T','U') 4364 grid__get_fine_offset_cc(jp_J,1)=id_rho(jp_J)-ij 4365 CASE DEFAULT !'F','V' 4366 grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij 4367 END SELECT 4368 ELSE 4369 ! odd 4370 grid__get_fine_offset_cc(jp_J,1)=(id_rho(jp_J)+1)-ij 4371 ENDIF 4372 4373 ELSE 4374 CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do"//& 4375 & " not match fine grid lower left corner.") 4376 ENDIF 4377 4378 IF( ll_greenwich )THEN 4379 ! close to greenwich meridien 4380 ll_greenwich=.FALSE. 4381 ! -180:180 => 0:360 4382 WHERE( dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) < 0. ) 4383 dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1) = & 4384 & dl_lon0(id_imin0:id_imin0+1,id_jmin0:id_jmin0+1)+360. 4385 END WHERE 4386 4387 WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 4388 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 4389 END WHERE 4390 ENDIF 4391 4392 !!!!!! look for upper right offset !!!!!! 4393 i1=MAX(1,il_shape1(jp_I)-(id_rho(jp_I)+2)+1) ; i2=il_shape1(jp_I) 4394 j1=MAX(1,il_shape1(jp_J)-(id_rho(jp_J)+2)+1) ; j2=il_shape1(jp_J) 4395 4396 ! check if cross greenwich meridien 4397 IF( minval(dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0))<5. .OR. & 4398 & maxval(dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0))>355. )THEN 4399 ! close to greenwich meridien 4400 ll_greenwich=.TRUE. 4401 ! 0:360 => -180:180 4402 WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) > 180. ) 4403 dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) = & 4404 & dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0)-360. 4405 END WHERE 4406 4407 WHERE( dl_lon1(i1:i2,j1:j2) > 180. ) 4408 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)-360. 4409 END WHERE 4410 ENDIF 4411 4412 ! min latitude of the upper right cell 4413 dl_lonmin0=MIN(dl_lon0(id_imax0-1,id_jmax0-1),dl_lon0(id_imax0-1,id_jmax0)) 4414 ! min latitude of the upper right cell 4415 dl_latmin0=MIN(dd_lat0(id_imax0-1,id_jmax0-1),dd_lat0(id_imax0,id_jmax0-1)) 4416 IF( dl_lon1(il_shape1(jp_I),il_shape1(jp_J)) > dl_lonmin0 .AND. & 4417 & dd_lat1(il_shape1(jp_I),il_shape1(jp_J)) > dl_latmin0 )THEN 4418 4419 !!!!! i-direction !!!!! 4420 IF( ll_even(jp_I) )THEN 4421 ! even 4422 SELECT CASE(TRIM(cl_point)) 4423 CASE('F','U') 4424 dl_dlon= ( dl_lon0(id_imax0 ,id_jmax0-1) - & 4425 & dl_lon0(id_imax0-1,id_jmax0-1) ) / & 4426 & ( 2.*id_rho(jp_I) ) 4427 CASE DEFAULT 4428 dl_dlon=0 4429 END SELECT 4430 ELSE 4431 ! odd 4432 dl_dlon= ( dl_lon0(id_imax0 ,id_jmax0-1) - & 4433 & dl_lon0(id_imax0-1,id_jmax0-1) ) / & 4434 & ( 2*id_rho(jp_I) ) 4435 ENDIF 4436 4437 !!!!! j-direction !!!!! 4438 IF( ll_even(jp_J) )THEN 4439 ! even 4440 SELECT CASE(TRIM(cl_point)) 4441 CASE('F','V') 4442 dl_dlat= ( dd_lat0(id_imax0-1,id_jmax0 ) - & 4443 & dd_lat0(id_imax0-1,id_jmax0-1) ) / & 4444 & ( 2.*id_rho(jp_J) ) 4445 CASE DEFAULT 4446 dl_dlat=0 4447 END SELECT 4448 ELSE 4449 ! odd 4450 dl_dlat= ( dd_lat0(id_imax0-1,id_jmax0 ) - & 4451 & dd_lat0(id_imax0-1,id_jmax0-1) ) / & 4452 & ( 2*id_rho(jp_J) ) 4453 ENDIF 4454 4455 dl_lon0F= dl_lon0(id_imax0-1,id_jmax0-1) - dl_dlon 4456 dl_lat0F= dd_lat0(id_imax0-1,id_jmax0-1) - dl_dlat 4457 4458 il_ind(:)=grid_get_closest( dl_lon1(i1:i2,j1:j2), dd_lat1(i1:i2,j1:j2), & 4459 & dl_lon0F, dl_lat0F, 'ur' ) 4460 4461 ii=(MIN(il_shape1(jp_I),(id_rho(jp_I)+2))-il_ind(1)+1) 4462 ij=(MIN(il_shape1(jp_J),(id_rho(jp_J)+2))-il_ind(2)+1) 4463 4464 !!!!! i-direction !!!!! 4465 IF( ll_even(jp_I) )THEN 4466 ! even 4467 SELECT CASE(TRIM(cl_point)) 4468 CASE('T','V') 4469 grid__get_fine_offset_cc(jp_I,2)=id_rho(jp_I)-ii 4470 CASE DEFAULT !'F','U' 4471 grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 4472 END SELECT 4473 ELSE 4474 ! odd 4475 grid__get_fine_offset_cc(jp_I,2)=(id_rho(jp_I)+1)-ii 4476 ENDIF 4477 4478 !!!!! j-direction !!!!! 4479 IF( ll_even(jp_J) )THEN 4480 ! even 4481 SELECT CASE(TRIM(cl_point)) 4482 CASE('T','U') 4483 grid__get_fine_offset_cc(jp_J,2)=id_rho(jp_J)-ij 4484 CASE DEFAULT !'F','V' 4485 grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 4486 END SELECT 4487 ELSE 4488 ! odd 4489 grid__get_fine_offset_cc(jp_J,2)=(id_rho(jp_J)+1)-ij 4490 ENDIF 4491 4492 ELSE 4493 CALL logger_error("GRID GET FINE OFFSET: coarse grid indices do"//& 3718 4494 & " not match fine grid upper right corner.") 3719 4495 ENDIF 3720 ELSE 3721 grid__get_fine_offset_cc(2,:)=((id_rho(jp_J)-1)/2) 4496 4497 IF( ll_greenwich )THEN 4498 ! close to greenwich meridien 4499 ll_greenwich=.FALSE. 4500 ! -180:180 => 0:360 4501 WHERE( dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) < 0. ) 4502 dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0) = & 4503 & dl_lon0(id_imax0-1:id_imax0,id_jmax0-1:id_jmax0)+360. 4504 END WHERE 4505 4506 WHERE( dl_lon1(i1:i2,j1:j2) < 0. ) 4507 dl_lon1(i1:i2,j1:j2)=dl_lon1(i1:i2,j1:j2)+360. 4508 END WHERE 4509 ENDIF 4510 3722 4511 ENDIF 3723 4512 3724 4513 DEALLOCATE( dl_lon0 ) 3725 4514 DEALLOCATE( dl_lon1 ) 4515 4516 IF( ANY(grid__get_fine_offset_cc(:,:)==-1) )THEN 4517 CALL logger_fatal("GRID GET FINE OFFSET: can not found "//& 4518 & " offset between coarse and fine grid.") 4519 ENDIF 3726 4520 3727 4521 END FUNCTION grid__get_fine_offset_cc … … 3732 4526 ! 3733 4527 !> @author J.Paul 3734 !> -November, 2013- Initial Version4528 !> @date November, 2013- Initial Version 3735 4529 !> @date October, 2014 3736 4530 !> - work on mpp file structure instead of file structure 3737 ! 4531 !> @date February, 2016 4532 !> - use F-point to check coincidence for even refinment 4533 !> - use F-point estimation, if can not read it. 4534 !> 3738 4535 !> @param[in] td_coord0 coarse grid coordinate file structure 3739 4536 !> @param[in] td_coord1 fine grid coordinate file structure … … 3742 4539 !> @param[in] id_jmin0 coarse grid lower left corner j-indice of fine grid domain 3743 4540 !> @param[in] id_jmax0 coarse grid upper right corner j-indice of fine grid domain 3744 !> @param[in] id_rho array of refinement factor (default 1)4541 !> @param[in] id_rho array of refinement factor 3745 4542 !------------------------------------------------------------------- 3746 4543 SUBROUTINE grid_check_coincidence( td_coord0, td_coord1, & … … 3760 4557 3761 4558 ! local variable 3762 INTEGER(i4) :: il_imid13763 INTEGER(i4) :: il_jmid14559 INTEGER(i4) :: il_imid1 4560 INTEGER(i4) :: il_jmid1 3764 4561 3765 INTEGER(i4) :: il_ew0 3766 INTEGER(i4) :: il_ew1 3767 3768 INTEGER(i4) :: il_imin1 3769 INTEGER(i4) :: il_imax1 3770 INTEGER(i4) :: il_jmin1 3771 INTEGER(i4) :: il_jmax1 3772 3773 INTEGER(i4), DIMENSION(2) :: il_indC 3774 INTEGER(i4), DIMENSION(2) :: il_indF 3775 INTEGER(i4), DIMENSION(2) :: il_iind 3776 INTEGER(i4), DIMENSION(2) :: il_jind 3777 3778 REAL(dp) :: dl_lon0 3779 REAL(dp) :: dl_lat0 3780 REAL(dp) :: dl_lon1 3781 REAL(dp) :: dl_lat1 3782 3783 REAL(dp) :: dl_lon1p 3784 REAL(dp) :: dl_lat1p 3785 3786 LOGICAL :: ll_coincidence 3787 3788 TYPE(TVAR) :: tl_lon0 3789 TYPE(TVAR) :: tl_lat0 3790 TYPE(TVAR) :: tl_lon1 3791 TYPE(TVAR) :: tl_lat1 3792 3793 TYPE(TMPP) :: tl_coord0 3794 TYPE(TMPP) :: tl_coord1 3795 3796 TYPE(TDOM) :: tl_dom0 4562 INTEGER(i4) :: il_ew0 4563 INTEGER(i4) :: il_ew1 4564 4565 INTEGER(i4) :: il_ind 4566 4567 INTEGER(i4) :: il_imin1 4568 INTEGER(i4) :: il_imax1 4569 INTEGER(i4) :: il_jmin1 4570 INTEGER(i4) :: il_jmax1 4571 4572 INTEGER(i4), DIMENSION(2) :: il_ind0 4573 INTEGER(i4), DIMENSION(2) :: il_ind1 4574 4575 INTEGER(i4), DIMENSION(2) :: il_ill1 4576 INTEGER(i4), DIMENSION(2) :: il_ilr1 4577 INTEGER(i4), DIMENSION(2) :: il_iul1 4578 INTEGER(i4), DIMENSION(2) :: il_iur1 4579 4580 REAL(dp) :: dl_lon0F 4581 REAL(dp) :: dl_lat0F 4582 REAL(dp) :: dl_lon0 4583 REAL(dp) :: dl_lat0 4584 REAL(dp) :: dl_lon1F 4585 REAL(dp) :: dl_lat1F 4586 REAL(dp) :: dl_lon1 4587 REAL(dp) :: dl_lat1 4588 4589 REAL(dp) :: dl_delta 4590 4591 LOGICAL :: ll_coincidence 4592 LOGICAL :: ll_even 4593 LOGICAL :: ll_grid0F 4594 LOGICAL :: ll_grid1F 4595 4596 TYPE(TVAR) :: tl_lon0 4597 TYPE(TVAR) :: tl_lat0 4598 TYPE(TVAR) :: tl_lon0F 4599 TYPE(TVAR) :: tl_lat0F 4600 TYPE(TVAR) :: tl_lon1 4601 TYPE(TVAR) :: tl_lat1 4602 TYPE(TVAR) :: tl_lon1F 4603 TYPE(TVAR) :: tl_lat1F 4604 4605 TYPE(TMPP) :: tl_coord0 4606 TYPE(TMPP) :: tl_coord1 4607 4608 TYPE(TDOM) :: tl_dom0 3797 4609 3798 4610 ! loop indices … … 3803 4615 ll_coincidence=.TRUE. 3804 4616 4617 ll_even=.FALSE. 4618 IF( MOD(id_rho(jp_I)*id_rho(jp_J),2) == 0 )THEN 4619 ll_even=.TRUE. 4620 ENDIF 4621 3805 4622 ! copy structure 3806 4623 tl_coord0=mpp_copy(td_coord0) … … 3815 4632 3816 4633 ! read variable value on domain 3817 tl_lon0=iom_dom_read_var(tl_coord0,'longitude',tl_dom0) 3818 tl_lat0=iom_dom_read_var(tl_coord0,'latitude' ,tl_dom0) 4634 il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'longitude_T') 4635 IF( il_ind /= 0 )THEN 4636 tl_lon0=iom_dom_read_var(tl_coord0,'longitude_T',tl_dom0) 4637 ELSE 4638 tl_lon0=iom_dom_read_var(tl_coord0,'longitude',tl_dom0) 4639 ENDIF 4640 4641 il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'latitude_T') 4642 IF( il_ind /= 0 )THEN 4643 tl_lat0=iom_dom_read_var(tl_coord0,'latitude_T' ,tl_dom0) 4644 ELSE 4645 tl_lat0=iom_dom_read_var(tl_coord0,'latitude' ,tl_dom0) 4646 ENDIF 4647 4648 IF( ll_even )THEN 4649 ! look for variable value on domain for F point 4650 il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'longitude_F') 4651 IF( il_ind /= 0 )THEN 4652 tl_lon0F=iom_dom_read_var(tl_coord0,'longitude_F',tl_dom0) 4653 ENDIF 4654 4655 il_ind=var_get_index(tl_coord0%t_proc(1)%t_var(:), 'latitude_F') 4656 IF( il_ind /= 0 )THEN 4657 tl_lat0F=iom_dom_read_var(tl_coord0,'latitude_F' ,tl_dom0) 4658 ENDIF 4659 4660 ll_grid0F=.FALSE. 4661 IF( ASSOCIATED(tl_lon0F%d_value) .AND. & 4662 & ASSOCIATED(tl_lat0F%d_value) )THEN 4663 ll_grid0F=.TRUE. 4664 ENDIF 4665 4666 ENDIF 3819 4667 3820 4668 ! close mpp files … … 3832 4680 3833 4681 ! read fine longitue and latitude 3834 tl_lon1=iom_mpp_read_var(tl_coord1,'longitude') 3835 tl_lat1=iom_mpp_read_var(tl_coord1,'latitude') 4682 il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), TRIM(tl_lon0%c_longname)) 4683 IF( il_ind /= 0 )THEN 4684 tl_lon1=iom_mpp_read_var(tl_coord1,TRIM(tl_lon0%c_longname)) 4685 ELSE 4686 tl_lon1=iom_mpp_read_var(tl_coord1,'longitude') 4687 ENDIF 4688 il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), TRIM(tl_lat0%c_longname)) 4689 IF( il_ind /= 0 )THEN 4690 tl_lat1=iom_mpp_read_var(tl_coord1,TRIM(tl_lat0%c_longname)) 4691 ELSE 4692 tl_lat1=iom_mpp_read_var(tl_coord1,'latitude') 4693 ENDIF 3836 4694 4695 IF( ll_even )THEN 4696 4697 ! look for variable value on domain for F point 4698 il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), 'longitude_F') 4699 IF( il_ind /= 0 )THEN 4700 tl_lon1F=iom_mpp_read_var(tl_coord1,'longitude_F') 4701 ENDIF 4702 4703 il_ind=var_get_index(tl_coord1%t_proc(1)%t_var(:), 'latitude_F') 4704 IF( il_ind /= 0 )THEN 4705 tl_lat1F=iom_mpp_read_var(tl_coord1,'latitude_F') 4706 ENDIF 4707 4708 ll_grid1F=.FALSE. 4709 IF( ASSOCIATED(tl_lon1F%d_value) .AND. & 4710 & ASSOCIATED(tl_lat1F%d_value) )THEN 4711 ll_grid1F=.TRUE. 4712 ENDIF 4713 4714 ENDIF 4715 3837 4716 ! close mpp files 3838 CALL iom_ dom_close(tl_coord1)4717 CALL iom_mpp_close(tl_coord1) 3839 4718 ! clean structure 3840 4719 CALL mpp_clean(tl_coord1) … … 3898 4777 IF( .NOT. ll_coincidence )THEN 3899 4778 CALL logger_fatal("GRID CHECK COINCIDENCE: no coincidence "//& 3900 & "between fine grid and coarse grid . invalid domain" )4779 & "between fine grid and coarse grid: invalid domain." ) 3901 4780 ENDIF 3902 4781 … … 3912 4791 3913 4792 ! select closest point on coarse grid 3914 il_ind C(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),&4793 il_ind0(:)=grid_get_closest(tl_lon0%d_value(:,:,1,1),& 3915 4794 & tl_lat0%d_value(:,:,1,1),& 3916 4795 & dl_lon1, dl_lat1 ) 3917 4796 3918 IF( ANY(il_ind C(:)==0) )THEN4797 IF( ANY(il_ind0(:)==0) )THEN 3919 4798 CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//& 3920 & "coarse grid indices. invalid domain" ) 3921 ENDIF 3922 3923 dl_lon0=tl_lon0%d_value(il_indC(1),il_indC(2),1,1) 3924 dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2),1,1) 3925 3926 ! look for closest fine grid point from selected coarse grid point 3927 il_iind(:)=MAXLOC( tl_lon1%d_value(:,:,1,1), & 3928 & tl_lon1%d_value(:,:,1,1) <= dl_lon0 ) 3929 3930 il_jind(:)=MAXLOC( tl_lat1%d_value(:,:,1,1), & 3931 & tl_lat1%d_value(:,:,1,1) <= dl_lat0 ) 3932 3933 il_indF(1)=il_iind(1) 3934 il_indF(2)=il_jind(2) 3935 3936 IF( ANY(il_indF(:)==0) )THEN 3937 CALL logger_fatal("GRID CHECK COINCIDENCE: can not find valid "//& 3938 & "fine grid indices. invalid domain" ) 3939 ENDIF 3940 3941 dl_lon1=tl_lon1%d_value(il_indF(1),il_indF(2),1,1) 3942 dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2),1,1) 3943 3944 ! check i-direction refinement factor 3945 DO ji=1,MIN(3,il_imid1) 3946 3947 IF( il_indF(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN 3948 CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//& 3949 & " to check i-direction refinement factor ") 3950 EXIT 3951 ELSE 3952 dl_lon1=tl_lon1%d_value(il_indF(1)+ji*id_rho(jp_I),il_indF(2),1,1) 3953 dl_lon0=tl_lon0%d_value(il_indC(1)+ji,il_indC(2),1,1) 3954 3955 dl_lon1p=tl_lon1%d_value(il_indF(1)+ji*id_rho(jp_I)+1,il_indF(2),1,1) 3956 3957 SELECT CASE(MOD(id_rho(jp_I),2)) 3958 3959 CASE(0) 3960 3961 IF( dl_lon1 >= dl_lon0 .OR. dl_lon0 >= dl_lon1p )THEN 3962 ll_coincidence=.FALSE. 3963 CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//& 3964 & "i-direction refinement factor ("//& 3965 & TRIM(fct_str(id_rho(jp_I)))//& 3966 & ") between fine grid and coarse grid ") 3967 ENDIF 3968 3969 CASE DEFAULT 3970 4799 & "coarse grid indices: invalid domain." ) 4800 ENDIF 4801 4802 IF( .NOT. ll_even )THEN 4803 ! case odd refinment in both direction 4804 ! work on T-point 4805 4806 dl_lon0=tl_lon0%d_value(il_ind0(1),il_ind0(2),1,1) 4807 dl_lat0=tl_lat0%d_value(il_ind0(1),il_ind0(2),1,1) 4808 4809 il_ind1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),& 4810 & tl_lat1%d_value(:,:,1,1),& 4811 & dl_lon0, dl_lat0 ) 4812 4813 ! check i-direction refinement factor 4814 DO ji=0,MIN(3,il_imid1) 4815 4816 IF( il_ind1(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN 4817 CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//& 4818 & " to check i-direction refinement factor ") 4819 EXIT 4820 ELSE 4821 dl_lon0=tl_lon0%d_value(il_ind0(1)+ji ,il_ind0(2),1,1) 4822 dl_lon1=tl_lon1%d_value(il_ind1(1)+ji*id_rho(jp_I),il_ind1(2),1,1) 4823 4824 ! assume there could be little difference due to interpolation 3971 4825 IF( ABS(dl_lon1 - dl_lon0) > dp_delta )THEN 3972 4826 ll_coincidence=.FALSE. … … 3976 4830 & ") between fine grid and coarse grid ") 3977 4831 ENDIF 3978 3979 END SELECT 3980 ENDIF 3981 3982 ENDDO 3983 3984 ! check j-direction refinement factor 3985 DO jj=1,MIN(3,il_jmid1) 3986 3987 IF( il_indF(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN 3988 CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//& 3989 & " to check j-direction refinement factor ") 3990 EXIT 3991 ELSE 3992 dl_lat1=tl_lat1%d_value(il_indF(1),il_indF(2)+jj*id_rho(jp_J),1,1) 3993 dl_lat0=tl_lat0%d_value(il_indC(1),il_indC(2)+jj,1,1) 3994 3995 dl_lat1p=tl_lat1%d_value(il_indF(1),il_indF(2)+jj*id_rho(jp_J)+1,1,1) 3996 3997 SELECT CASE(MOD(id_rho(jp_J),2)) 3998 3999 CASE(0) 4000 4001 IF( dl_lat1 >= dl_lat0 .OR. dl_lat0 >= dl_lat1p )THEN 4002 ll_coincidence=.FALSE. 4003 CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//& 4004 & "j-direction refinement factor ("//& 4005 & TRIM(fct_str(id_rho(jp_J)))//& 4006 & ") between fine grid and coarse grid ") 4007 ENDIF 4008 4009 CASE DEFAULT 4010 4832 ENDIF 4833 4834 ENDDO 4835 4836 ! check j-direction refinement factor 4837 DO jj=0,MIN(3,il_jmid1) 4838 4839 IF( il_ind1(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN 4840 CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//& 4841 & " to check j-direction refinement factor ") 4842 EXIT 4843 ELSE 4844 dl_lat0=tl_lat0%d_value(il_ind0(1),il_ind0(2)+jj ,1,1) 4845 dl_lat1=tl_lat1%d_value(il_ind1(1),il_ind1(2)+jj*id_rho(jp_J),1,1) 4846 4847 ! assume there could be little difference due to interpolation 4011 4848 IF( ABS(dl_lat1-dl_lat0) > dp_delta )THEN 4012 4849 ll_coincidence=.FALSE. … … 4016 4853 & ") between fine grid and coarse grid ") 4017 4854 ENDIF 4018 4019 END SELECT 4020 ENDIF 4021 4022 ENDDO 4855 ENDIF 4856 4857 ENDDO 4858 4859 ELSE 4860 ! case even refinment at least in one direction 4861 ! work on F-point 4862 4863 dl_delta=dp_delta 4864 ! look for lower left fine point in coarse cell. 4865 IF( ll_grid0F )THEN 4866 4867 ! lower left corner of coarse cell 4868 dl_lon0F=tl_lon0F%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1) 4869 dl_lat0F=tl_lat0F%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1) 4870 4871 ELSE 4872 4873 ! approximate lower left corner of coarse cell (with T point) 4874 dl_lon0F=( tl_lon0%d_value(il_ind0(1) ,il_ind0(2) ,1,1) + & 4875 & tl_lon0%d_value(il_ind0(1) ,il_ind0(2)-1,1,1) + & 4876 & tl_lon0%d_value(il_ind0(1)-1,il_ind0(2) ,1,1) + & 4877 & tl_lon0%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1) ) * 0.25 4878 4879 dl_lat0F=( tl_lat0%d_value(il_ind0(1) ,il_ind0(2) ,1,1) + & 4880 & tl_lat0%d_value(il_ind0(1) ,il_ind0(2)-1,1,1) + & 4881 & tl_lat0%d_value(il_ind0(1)-1,il_ind0(2) ,1,1) + & 4882 & tl_lat0%d_value(il_ind0(1)-1,il_ind0(2)-1,1,1) ) * 0.25 4883 4884 ! as we use approximation of F-point we relax condition 4885 dl_delta=100*dp_delta 4886 4887 ENDIF 4888 4889 IF( ll_grid1F )THEN 4890 4891 il_ind1(:)=grid_get_closest(tl_lon1F%d_value(:,:,1,1),& 4892 & tl_lat1F%d_value(:,:,1,1),& 4893 & dl_lon0F, dl_lat0F ) 4894 4895 ELSE 4896 4897 il_ill1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),& 4898 & tl_lat1%d_value(:,:,1,1),& 4899 & dl_lon0F, dl_lat0F, 'll' ) 4900 4901 il_ilr1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),& 4902 & tl_lat1%d_value(:,:,1,1),& 4903 & dl_lon0F, dl_lat0F, 'lr' ) 4904 4905 il_iul1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),& 4906 & tl_lat1%d_value(:,:,1,1),& 4907 & dl_lon0F, dl_lat0F, 'ul' ) 4908 4909 il_iur1(:)=grid_get_closest(tl_lon1%d_value(:,:,1,1),& 4910 & tl_lat1%d_value(:,:,1,1),& 4911 & dl_lon0F, dl_lat0F, 'ur' ) 4912 4913 ! as we use approximation of F-point we relax condition 4914 dl_delta=100*dp_delta 4915 4916 ENDIF 4917 4918 ! check i-direction refinement factor 4919 DO ji=0,MIN(3,il_imid1) 4920 4921 IF( il_ind1(1)+ji*id_rho(jp_I)+1 > tl_lon1%t_dim(1)%i_len )THEN 4922 CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//& 4923 & " to check i-direction refinement factor ") 4924 EXIT 4925 ELSE 4926 IF( ll_grid0F )THEN 4927 dl_lon0F=tl_lon0F%d_value(il_ind0(1)+ji-1, il_ind0(2)-1,1,1) 4928 ELSE 4929 dl_lon0F= 0.25 * & 4930 & ( tl_lon0%d_value(il_ind0(1)+ji , il_ind0(2) ,1,1) + & 4931 & tl_lon0%d_value(il_ind0(1)+ji-1, il_ind0(2) ,1,1) + & 4932 & tl_lon0%d_value(il_ind0(1)+ji , il_ind0(2)-1,1,1) + & 4933 & tl_lon0%d_value(il_ind0(1)+ji-1, il_ind0(2)-1,1,1) ) 4934 ENDIF 4935 4936 IF( ll_grid1F )THEN 4937 dl_lon1F= tl_lon1F%d_value( il_ind1(1)+ji*id_rho(jp_I), & 4938 & il_ind1(2),1,1) 4939 ELSE 4940 dl_lon1F= 0.25 * & 4941 & ( tl_lon1%d_value( il_ill1(1)+ji*id_rho(jp_I), & 4942 & il_ill1(2),1,1) + & 4943 & tl_lon1%d_value( il_ilr1(1)+ji*id_rho(jp_I), & 4944 & il_ilr1(2),1,1) + & 4945 & tl_lon1%d_value( il_iul1(1)+ji*id_rho(jp_I), & 4946 & il_iul1(2),1,1) + & 4947 & tl_lon1%d_value( il_iur1(1)+ji*id_rho(jp_I), & 4948 & il_iur1(2),1,1) ) 4949 4950 ENDIF 4951 4952 ! assume there could be little difference due to interpolation 4953 IF( ABS(dl_lon1F - dl_lon0F) > dl_delta )THEN 4954 ll_coincidence=.FALSE. 4955 CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//& 4956 & "i-direction refinement factor ("//& 4957 & TRIM(fct_str(id_rho(jp_I)))//& 4958 & ") between fine grid and coarse grid ") 4959 ENDIF 4960 ENDIF 4961 4962 ENDDO 4963 4964 ! check j-direction refinement factor 4965 DO jj=0,MIN(3,il_jmid1) 4966 4967 IF( il_ind1(2)+jj*id_rho(jp_J)+1 > tl_lat1%t_dim(2)%i_len )THEN 4968 CALL logger_warn("GRID CHECK COINCIDENCE: domain to small "//& 4969 & " to check j-direction refinement factor ") 4970 EXIT 4971 ELSE 4972 IF( ll_grid0F )THEN 4973 dl_lat0F=tl_lat0F%d_value(il_ind0(1)-1, il_ind0(2)+jj-1,1,1) 4974 ELSE 4975 dl_lat0F= 0.25 * & 4976 & ( tl_lat0%d_value(il_ind0(1) , il_ind0(2)+jj ,1,1) + & 4977 & tl_lat0%d_value(il_ind0(1)-1, il_ind0(2)+jj ,1,1) + & 4978 & tl_lat0%d_value(il_ind0(1) , il_ind0(2)+jj-1,1,1) + & 4979 & tl_lat0%d_value(il_ind0(1)-1, il_ind0(2)+jj-1,1,1) ) 4980 ENDIF 4981 4982 IF( ll_grid1F )THEN 4983 dl_lat1F= tl_lat1F%d_value( il_ind1(1), & 4984 & il_ind1(2)+jj*id_rho(jp_J),1,1) 4985 ELSE 4986 dl_lat1F= 0.25 * & 4987 & ( tl_lat1%d_value( il_ill1(1), & 4988 & il_ill1(2)+jj*id_rho(jp_J),1,1) + & 4989 & tl_lat1%d_value( il_ilr1(1), & 4990 & il_ilr1(2)+jj*id_rho(jp_J),1,1) + & 4991 & tl_lat1%d_value( il_iul1(1), & 4992 & il_iul1(2)+jj*id_rho(jp_J),1,1) + & 4993 & tl_lat1%d_value( il_iur1(1), & 4994 & il_iur1(2)+jj*id_rho(jp_J),1,1) ) 4995 4996 ENDIF 4997 4998 ! assume there could be little difference due to interpolation 4999 IF( ABS(dl_lat1F - dl_lat0F) > dl_delta )THEN 5000 ll_coincidence=.FALSE. 5001 CALL logger_debug("GRID CHECK COINCIDENCE: invalid "//& 5002 & "i-direction refinement factor ("//& 5003 & TRIM(fct_str(id_rho(jp_I)))//& 5004 & ") between fine grid and coarse grid ") 5005 ENDIF 5006 ENDIF 5007 5008 ENDDO 5009 ENDIF 4023 5010 4024 5011 ! clean … … 4042 5029 !> 4043 5030 !> @author J.Paul 4044 !> - November, 2013- Initial Version5031 !> @date November, 2013 - Initial Version 4045 5032 ! 4046 5033 !> @param[in] dd_lon0 array of coarse grid longitude … … 4103 5090 dl_lon1 = dd_lon1(il_imin1, il_jmin1) 4104 5091 dl_lat1 = dd_lat1(il_imin1, il_jmin1) 4105 4106 5092 4107 5093 IF( (ABS(dl_lon1-dl_lon0)>dp_delta) .AND. (dl_lon1 < dl_lon0 ) .OR. & … … 4202 5188 ! 4203 5189 !> @author J.Paul 4204 !> - November, 2013- Initial Version5190 !> @date November, 2013 - Initial Version 4205 5191 ! 4206 5192 !> @param[in] dd_lat0 array of coarse grid latitude … … 4272 5258 !> 4273 5259 !> @author J.Paul 4274 !> - November, 2013-Initial version5260 !> @date November, 2013 - Initial version 4275 5261 ! 4276 5262 !> @param[inout] td_var array of variable structure … … 4348 5334 !> 4349 5335 !> @author J.Paul 4350 !> - November, 2013-Initial version5336 !> @date November, 2013 - Initial version 4351 5337 ! 4352 5338 !> @param[inout] td_var array of variable structure … … 4374 5360 IF( ALL(td_var%t_dim(1:2)%l_use) )THEN 4375 5361 4376 CALL logger_warn( "GRID DEL GHOST: dimension change in variable "//& 4377 & TRIM(td_var%c_name) ) 5362 IF( ANY(id_ghost(:,:)/=0) )THEN 5363 CALL logger_warn( "GRID DEL GHOST: dimension change in variable "//& 5364 & TRIM(td_var%c_name) ) 5365 ENDIF 4378 5366 4379 5367 ! copy variable … … 4425 5413 !> 4426 5414 !> @author J.Paul 4427 !> - September, 2014- Initial Version5415 !> @date September, 2014 - Initial Version 4428 5416 ! 4429 5417 !> @param[in] td_var variable sturcture … … 4555 5543 !> 4556 5544 !> @author J.Paul 4557 !> -September, 2014 - Initial Version5545 !> @date September, 2014 - Initial Version 4558 5546 !> @date October, 2014 4559 5547 !> - work on mpp file structure instead of file structure … … 4590 5578 4591 5579 ! copy structure 4592 tl_mpp=mpp_copy(td_mpp) 4593 4594 IF( tl_mpp%i_perio < 0 )THEN 4595 ! compute NEMO periodicity index 4596 CALL grid_get_info(tl_mpp) 4597 ENDIF 5580 tl_mpp=mpp_copy(td_mpp) 5581 5582 CALL logger_info("GRID GET FINE GHOST perio"//TRIM(fct_str(tl_mpp%i_perio))) 5583 IF( tl_mpp%i_perio < 0 )THEN 5584 ! compute NEMO periodicity index 5585 CALL grid_get_info(tl_mpp) 5586 ENDIF 4598 5587 4599 5588 SELECT CASE(tl_mpp%i_perio) … … 4627 5616 !> 4628 5617 !> @author J.Paul 4629 !> - November, 2013- Initial Version5618 !> @date November, 2013 - Initial Version 4630 5619 ! 4631 5620 !> @param[in] td_var variable strucutre … … 4694 5683 il_tmp(jim:jip,jjm:jjp)=1 4695 5684 END WHERE 5685 4696 5686 ENDIF 4697 5687 ENDDO … … 4720 5710 !> 4721 5711 !> @details 4722 !> the minimum size (n bumber of point) of closed sea to be kept could be5712 !> the minimum size (number of point) of closed sea to be kept could be 4723 5713 !> sepcify with id_minsize. 4724 5714 !> By default only the biggest sea is preserve. 4725 5715 !> 4726 5716 !> @author J.Paul 4727 !> - November, 2013- Initial Version5717 !> @date November, 2013 - Initial Version 4728 5718 !> 4729 5719 !> @param[inout] td_var variable structure … … 4782 5772 4783 5773 END SUBROUTINE grid_fill_small_dom 5774 !------------------------------------------------------------------- 5775 !> @brief This subroutine fill small domain inside bigger one. 5776 !> 5777 !> @details 5778 !> the minimum size (number of point) of domain sea to be kept could be 5779 !> is sepcified with id_minsize. 5780 !> smaller domain are included in the one they are embedded. 5781 !> 5782 !> @author J.Paul 5783 !> @date Ferbruay, 2015 - Initial Version 5784 !> 5785 !> @param[inout] id_mask domain mask (from grid_split_domain) 5786 !> @param[in] id_minsize minimum size of sea to be kept 5787 !------------------------------------------------------------------- 5788 SUBROUTINE grid_fill_small_msk(id_mask, id_minsize) 5789 IMPLICIT NONE 5790 ! Argument 5791 INTEGER(i4), DIMENSION(:,:), INTENT(INOUT) :: id_mask 5792 INTEGER(i4), INTENT(IN ) :: id_minsize 5793 5794 ! local variable 5795 INTEGER(i4) :: il_ndom 5796 INTEGER(i4) :: il_minsize 5797 INTEGER(i4) :: il_msk 5798 5799 INTEGER(i4) :: jim 5800 INTEGER(i4) :: jjm 5801 INTEGER(i4) :: jip 5802 INTEGER(i4) :: jjp 5803 5804 INTEGER(i4), DIMENSION(2) :: il_shape 5805 INTEGER(i4), DIMENSION(:,:), ALLOCATABLE :: il_tmp 5806 5807 ! loop indices 5808 INTEGER(i4) :: ii 5809 INTEGER(i4) :: ij 5810 5811 INTEGER(i4) :: ji 5812 INTEGER(i4) :: jj 5813 !---------------------------------------------------------------- 5814 5815 il_shape(:)=SHAPE(id_mask(:,:)) 5816 il_ndom=MINVAL(id_mask(:,:)) 5817 5818 ALLOCATE( il_tmp(il_shape(1),il_shape(2)) ) 5819 il_tmp(:,:)=0 5820 DO ji=-1,il_ndom,-1 5821 WHERE( id_mask(:,:)==ji ) 5822 il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji 5823 END WHERE 5824 ENDDO 5825 5826 DO WHILE( id_minsize > MINVAL(il_tmp(:,:)) ) 5827 5828 DO jj=1,il_shape(2) 5829 DO ji=1,il_shape(1) 5830 5831 IF( il_tmp(ji,jj) < il_minsize )THEN 5832 jim=MAX(1,ji-1) ; jip=MIN(il_shape(1),ji+1) 5833 jjm=MAX(1,jj-1) ; jjp=MIN(il_shape(2),jj+1) 5834 5835 il_msk=0 5836 DO ij=jjm,jjp 5837 DO ii=jim,jip 5838 IF( id_mask(ii,ij) /= id_mask(ji,jj) )THEN 5839 IF( il_msk == 0 )THEN 5840 il_msk=id_mask(ii,ij) 5841 ELSEIF( il_msk /= id_mask(ii,ij) )THEN 5842 CALL logger_error("GRID FILL SMALL MSK: "//& 5843 & "small domain not embedded in bigger one"//& 5844 & ". point should be between two different"//& 5845 & " domain.") 5846 ENDIF 5847 ENDIF 5848 ENDDO 5849 ENDDO 5850 IF( il_msk /= 0 ) id_mask(ji,jj)=il_msk 5851 5852 ENDIF 5853 5854 ENDDO 5855 ENDDO 5856 5857 5858 il_tmp(:,:)=0 5859 DO ji=-1,il_ndom,-1 5860 WHERE( id_mask(:,:)==ji ) 5861 il_tmp(:,:)=SUM(id_mask(:,:),id_mask(:,:)==ji)/ji 5862 END WHERE 5863 ENDDO 5864 5865 ENDDO 5866 5867 DEALLOCATE( il_tmp ) 5868 5869 5870 END SUBROUTINE grid_fill_small_msk 4784 5871 END MODULE grid 4785 5872
Note: See TracChangeset
for help on using the changeset viewer.
