- Timestamp:
- 2024-01-22T17:05:12+01:00 (5 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing_native_flow.f90
r8365 r8394 18 18 PUBLIC :: routing_flow_xios_initialize, routing_flow_set, routing_flow_get, routing_flow_main 19 19 PUBLIC :: routing_flow_initialize, routing_flow_finalize, routing_flow_clear, routing_flow_make_mean 20 PUBLIC :: routing_flow_diags 20 21 PUBLIC :: compute_coastline 21 22 … … 45 46 46 47 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_reservoir_r(:) !! Water amount in the stream reservoir (kg) - (local routing grid) 47 !$OMP THREADPRIVATE(stream_reservoir_r) 48 !$OMP THREADPRIVATE(stream_reservoir_r) 49 50 51 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: fast_diag(:) !! Diag on water amount in the fast reservoir (kg/m^2) - (orchidee grid) 52 !$OMP THREADPRIVATE(fast_diag) 53 54 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET:: slow_diag(:) !! Diag on water amount in the slow reservoir (kg/m^2) - (orchidee grid) 55 !$OMP THREADPRIVATE(slow_diag) 56 57 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: stream_diag(:) !! Diag on water amount in the stream reservoir (kg/m^2) - (orchidee grid) 58 !$OMP THREADPRIVATE(stream_diag) 59 60 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED,TARGET :: hydrographs_diag(:) !! Diag on water amount in the stream reservoir (kg/m^2) - (orchidee grid) 61 !$OMP THREADPRIVATE(hydrographs_diag) 62 48 63 49 64 REAL(r_std),SAVE,ALLOCATABLE,PUBLIC,PROTECTED :: riverflow_mean(:) !! Water amount in the stream reservoir (kg) - (local routing grid) … … 252 267 253 268 SUBROUTINE routing_flow_initialize(kjit, rest_id, nbpt_, dt_routing, contfrac, nbpt_r_, riverflow, coastalflow) 269 USE grid, ONLY : area 270 USE xios 254 271 IMPLICIT NONE 255 272 INTEGER ,INTENT(IN) :: kjit … … 273 290 ! 274 291 ALLOCATE (lakeinflow_mean(nbpt), stat=ier) 275 IF (ier /= 0) CALL ipslerr_p(3,'routing_ simple_init_1','Pb in allocate for lakeinflow_mean','','')292 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize','Pb in allocate for lakeinflow_mean','','') 276 293 var_name = 'lakeinflow' 277 294 CALL ioconf_setatt_p('UNITS', 'Kg/dt') … … 281 298 282 299 ALLOCATE (riverflow_mean(nbpt), stat=ier) 283 IF (ier /= 0) CALL ipslerr_p(3,'routing_ simple_init_1','Pb in allocate for riverflow_mean','','')300 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize','Pb in allocate for riverflow_mean','','') 284 301 var_name = 'riverflow' 285 302 CALL ioconf_setatt_p('UNITS', 'Kg/dt') … … 289 306 290 307 ALLOCATE (coastalflow_mean(nbpt), stat=ier) 291 IF (ier /= 0) CALL ipslerr_p(3,'routing_ simple_init_1','Pb in allocate for coastalflow_mean','','')308 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize','Pb in allocate for coastalflow_mean','','') 292 309 var_name = 'coastalflow' 293 310 CALL ioconf_setatt_p('UNITS', 'Kg/dt') … … 299 316 coastalflow(:) = coastalflow_mean(:) 300 317 301 END SUBROUTINE routing_flow_initialize 318 CALL routing_flow_initialize_diag(kjit, rest_id, contfrac) 319 320 END SUBROUTINE routing_flow_initialize 321 322 SUBROUTINE routing_flow_initialize_diag(kjit, rest_id, contfrac) 323 USE xios 324 USE grid, ONLY : area 325 IMPLICIT NONE 326 INTEGER ,INTENT(IN) :: kjit 327 INTEGER ,INTENT(IN) :: rest_id 328 REAL(r_std), INTENT(IN) :: contfrac(nbpt) !! fraction of land 329 INTEGER :: ier 330 CHARACTER(LEN=80) :: var_name !! To store variables names for I/O (unitless) 331 332 333 REAL(r_std) :: contfrac_mpi(nbp_mpi) 334 REAL(r_std) :: area_mpi(nbp_mpi) 335 REAL(r_std) :: fast_diag_mpi(nbp_mpi) 336 REAL(r_std) :: slow_diag_mpi(nbp_mpi) 337 REAL(r_std) :: stream_diag_mpi(nbp_mpi) 338 REAL(r_std) :: hydrographs_diag_mpi(nbp_mpi) 339 340 CALL gather_omp(contfrac, contfrac_mpi) 341 CALL gather_omp(area, area_mpi) 342 343 IF (is_omp_root) THEN 344 ALLOCATE(fast_diag(nbpt), stat=ier) 345 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for fast_diag','','') 346 ALLOCATE(slow_diag(nbpt), stat=ier) 347 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for slow_diag','','') 348 ALLOCATE(stream_diag(nbpt), stat=ier) 349 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for stream_diag','','') 350 351 352 CALL xios_send_field("routing_fast_diag0_r", fast_reservoir_r) 353 CALL xios_recv_field("routing_fast_diag0", fast_diag_mpi) 354 CALL xios_send_field("routing_slow_diag0_r", slow_reservoir_r) 355 CALL xios_recv_field("routing_slow_diag0", slow_diag_mpi) 356 CALL xios_send_field("routing_stream_diag0_r", stream_reservoir_r) 357 CALL xios_recv_field("routing_stream_diag0", stream_diag_mpi) 358 359 fast_diag_mpi=fast_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2 360 slow_diag_mpi=slow_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2 361 stream_diag_mpi=stream_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2 362 ENDIF 363 364 CALL scatter_omp(fast_diag_mpi,fast_diag) 365 CALL scatter_omp(slow_diag_mpi,slow_diag) 366 CALL scatter_omp(stream_diag_mpi,stream_diag) 367 368 ALLOCATE(hydrographs_diag(nbpt), stat=ier) 369 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow_initialize_diag','Pb in allocate for hydrographs_diag','','') 370 371 var_name = 'hydrographs' 372 CALL ioconf_setatt_p('UNITS', 'kg/dt_sechiba') 373 CALL ioconf_setatt_p('LONG_NAME','Hydrograph at outlow of grid') 374 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs_diag, "gather", nbp_glo, index_g) 375 CALL setvar_p (hydrographs_diag, val_exp, 'NO_KEYWORD', zero) 376 377 END SUBROUTINE routing_flow_initialize_diag 302 378 303 379 … … 1636 1712 1637 1713 REAL,ALLOCATABLE :: station_lonlat(:,:) 1638 INTEGER,ALLOCATABLE :: station_prec(:) 1639 1714 REAL,ALLOCATABLE :: station_prec(:) 1715 REAL,ALLOCATABLE :: station_area(:) 1716 INTEGER :: station_index_prec, station_index_prec_i, station_index_prec_j 1717 INTEGER :: station_index_area, station_index_area_i, station_index_area_j 1718 LOGICAL :: has_lonlat 1719 LOGICAL :: has_area 1640 1720 INTEGER ::i,j,r,k 1641 1721 … … 1651 1731 1652 1732 CHARACTER(LEN=256) :: str_station_ind 1653 REAL :: lon_a,lat_a,lon_b,lat_b,dist,max_area, max_max_area, min_dist 1654 INTEGER :: min_dist_index 1733 REAL :: lon_a,lat_a,lon_b,lat_b,dist,max_area, max_max_area, min_dist, min_area, min_min_area 1734 INTEGER :: min_dist_index, min_dist_index_i, min_dist_index_j 1655 1735 INTEGER :: ierr 1656 1736 … … 1666 1746 ALLOCATE(station_index(nb_station)) 1667 1747 ALLOCATE(station_prec(nb_station)) 1668 1748 ALLOCATE(station_area(nb_station)) 1749 1669 1750 DO k=1,nb_station 1670 1751 WRITE(str_station_ind,*) k 1671 1752 str_station_ind=ADJUSTL(str_station_ind) 1672 1753 CALL getin("station"//TRIM(str_station_ind)//"_id",station(k)) 1754 station_lonlat(k,:) = 1.1e10 1755 has_lonlat=.TRUE. 1673 1756 CALL getin("station"//TRIM(str_station_ind)//"_coor",station_lonlat(k,:)) 1757 IF (station_lonlat(k,1)>1e10 .OR. station_lonlat(k,2)>1e10) has_lonlat=.FALSE. 1674 1758 station_prec(k)=50000 1675 1759 CALL getin("station"//TRIM(str_station_ind)//"_prec",station_prec(k)) 1676 1760 station_area(k) = -1 1761 has_area=.TRUE. 1762 CALL getin("station"//TRIM(str_station_ind)//"_area",station_area(k)) 1763 IF (station_area(k)==-1) has_area=.FALSE. 1764 1677 1765 lon_a = station_lonlat(k,1)*Pi/180. 1678 1766 lat_a = station_lonlat(k,2)*Pi/180. … … 1680 1768 max_area=0 1681 1769 min_dist=HUGE(min_dist) 1770 min_area=HUGE(min_area) 1771 1682 1772 DO j=1,nj 1683 1773 DO i=1,ni … … 1691 1781 min_dist=dist 1692 1782 min_dist_index = r 1783 min_dist_index_i = i 1784 min_dist_index_j = j 1693 1785 ENDIF 1694 1786 1695 1787 IF (basins_area_r(r) > max_area) THEN 1696 1788 max_area = basins_area_r(r) 1697 station_index(k) = r 1789 station_index_prec = r 1790 station_index_prec_i = i 1791 station_index_prec_j = j 1698 1792 ENDIF 1793 1794 IF (has_area) THEN 1795 IF (ABS(basins_area_r(r)-station_area(k)*1e6) < min_area) THEN 1796 min_area = ABS(basins_area_r(r)-station_area(k)*1e6) 1797 station_index_area = r 1798 station_index_area_i = i 1799 station_index_area_j = j 1800 ENDIF 1801 ENDIF 1802 1699 1803 ENDIF 1700 1804 ENDIF … … 1703 1807 1704 1808 CALL MPI_ALLREDUCE(max_area, max_max_area, 1, MPI_REAL_ORCH, MPI_MAX, MPI_COMM_ORCH, ierr) 1705 IF (max_area /= max_max_area) station_index(k)=-1 1706 IF (max_area == 0) station_index(k)=-1 1707 1708 ! PRINT*,"Station ",TRIM(station(k))," old coor :",lon(min_dist_index),lat(min_dist_index)," basin_area (km2)", basins_area_r(min_dist_index)/1000/1000 1709 ! PRINT*,"Station ",TRIM(station(k)),"new coor ",lon(station_index(k)),lat(station_index(k))," basin_area(km2) ", basins_area_r(station_index(k))/1000/1000 1710 ENDDO 1809 IF (max_area /= max_max_area) station_index_prec=-1 1810 IF (max_area == 0) station_index_prec=-1 1811 1812 IF (has_area) THEN 1813 CALL MPI_ALLREDUCE(min_area, min_min_area, 1, MPI_REAL_ORCH, MPI_MIN, MPI_COMM_ORCH, ierr) 1814 IF (min_area /= min_min_area) station_index_area=-1 1815 IF (min_area == HUGE(min_area)) station_index_area=-1 1816 ELSE 1817 station_index_area=-1 1818 ENDIF 1819 1820 IF (station_index_prec/=-1) THEN 1821 PRINT*,"Found station ",TRIM(station(k))," based on coordinates lon-lat",lon(min_dist_index_i),lat(min_dist_index_j), & 1822 " with basin_area (km2) ", basins_area_r(min_dist_index)/1000/1000 1823 PRINT*,"=> at coordinate ",lon(station_index_prec_i),lat(station_index_prec_j)," basin_area(km2) ", & 1824 basins_area_r(station_index_prec)/1000/1000 1825 ENDIF 1826 1827 IF (station_index_area/=-1) THEN 1828 PRINT*,"Found station ",TRIM(station(k))," based on basin area (km2)", station_area(k) 1829 PRINT*,"=> at coordinate ",lon(station_index_area_i),lat(station_index_area_j)," basin_area(km2) ", & 1830 basins_area_r(station_index_area)/1000/1000 1831 ENDIF 1832 1833 IF (station_index_area/=-1 ) THEN 1834 station_index(k) = station_index_area 1835 ELSE 1836 station_index(k) = station_index_prec 1837 ENDIF 1838 1839 ENDDO 1711 1840 1712 1841 CALL xios_get_calendar_type(calendar_type) … … 1749 1878 END SUBROUTINE initialize_stations 1750 1879 1751 SUBROUTINE routing_flow_main(dt_routing ) !1880 SUBROUTINE routing_flow_main(dt_routing, contfrac) ! 1752 1881 1753 1882 USE xios … … 1760 1889 !! 0.1 Input variables 1761 1890 REAL(r_std), INTENT (in) :: dt_routing !! Routing time step (s) 1891 REAL(r_std), INTENT (in) :: contfrac(nbpt) !! Routing time step (s) 1762 1892 1763 1893 !! 0.4 Local variables … … 1769 1899 REAL(r_std) :: coastalflow(nbp_mpi) 1770 1900 REAL(r_std) :: lakeinflow(nbp_mpi) 1771 REAL(r_std) :: fast_diag_mpi(nbp_mpi)1772 REAL(r_std) :: slow_diag_mpi(nbp_mpi)1773 REAL(r_std) :: stream_diag_mpi(nbp_mpi)1774 1901 REAL(r_std) :: area_mpi(nbp_mpi) ! cell area 1902 REAL(r_std) :: contfrac_mpi(nbp_mpi) ! cell area 1775 1903 REAL(r_std) :: flow_coast(nbp_mpi) 1776 1904 REAL(r_std) :: flow_lake(nbp_mpi) 1777 1905 1906 ! diags 1907 REAL(r_std) :: fast_diag_mpi(nbp_mpi) 1908 REAL(r_std) :: slow_diag_mpi(nbp_mpi) 1909 REAL(r_std) :: stream_diag_mpi(nbp_mpi) 1910 REAL(r_std) :: hydrographs_diag_mpi(nbp_mpi) 1911 1778 1912 ! from input model -> routing_grid 1779 1913 REAL(r_std) :: runoff_r(nbpt_r) !! Grid-point runoff (kg/m^2/dt) … … 1823 1957 CALL gather_omp(drainage_mean, drainage) 1824 1958 CALL gather_omp(area, area_mpi) 1825 1959 CALL gather_omp(contfrac, contfrac_mpi) 1826 1960 1827 1961 IF (is_omp_root) THEN … … 2092 2226 CALL xios_send_field(TRIM(station(k)),value) 2093 2227 ENDDO 2094 CALL xios_orchidee_change_context("orchidee") 2228 2229 CALL xios_orchidee_change_context("orchidee") !! return on orchidee context 2230 2231 2232 !!! for orchidee history diag 2233 2234 CALL xios_send_field("routing_fast_diag_r", fast_reservoir_r) 2235 CALL xios_recv_field("routing_fast_diag", fast_diag_mpi) 2236 CALL xios_send_field("routing_slow_diag_r", slow_reservoir_r) 2237 CALL xios_recv_field("routing_slow_diag", slow_diag_mpi) 2238 CALL xios_send_field("routing_stream_diag_r", stream_reservoir_r) 2239 CALL xios_recv_field("routing_stream_diag", stream_diag_mpi) 2240 CALL xios_send_field("routing_hydrographs_diag_r", hydrographs_r+lakeinflow_r+coastalflow_r+riverflow_r) 2241 CALL xios_recv_field("routing_hydrographs_diag", hydrographs_diag_mpi) 2242 2243 fast_diag_mpi=fast_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2 2244 slow_diag_mpi=slow_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2 2245 stream_diag_mpi=stream_diag_mpi/(area_mpi*contfrac_mpi) !! kg => kg/m^2 2246 hydrographs_diag_mpi = hydrographs_diag_mpi/1000. 2247 2095 2248 ENDIF ! is_omp_root 2096 2249 2097 CALL scatter_omp(riverflow,riverflow_mean) 2098 CALL scatter_omp(coastalflow,coastalflow_mean) 2099 CALL scatter_omp(lakeinflow,lakeinflow_mean) 2100 2250 CALL scatter_omp(riverflow, riverflow_mean) 2251 CALL scatter_omp(coastalflow, coastalflow_mean) 2252 CALL scatter_omp(lakeinflow, lakeinflow_mean) 2253 CALL scatter_omp(fast_diag_mpi, fast_diag) 2254 CALL scatter_omp(slow_diag_mpi, slow_diag) 2255 CALL scatter_omp(stream_diag_mpi, stream_diag) 2256 CALL scatter_omp(hydrographs_diag_mpi, hydrographs_diag) 2257 2101 2258 CALL routing_flow_reset_mean 2102 2259 2103 2260 END SUBROUTINE routing_flow_main 2261 2262 SUBROUTINE routing_flow_diags(dt_sechiba) 2263 USE xios 2264 IMPLICIT NONE 2265 REAL(r_std) :: dt_sechiba 2266 2267 CALL xios_orchidee_send_field("fastr",fast_diag) 2268 CALL xios_orchidee_send_field("slowr",slow_diag) 2269 CALL xios_orchidee_send_field("streamr",stream_diag) 2270 CALL xios_orchidee_send_field("hydrographs",hydrographs_diag*dt_sechiba) 2271 2272 END SUBROUTINE routing_flow_diags 2104 2273 2105 2274 … … 2132 2301 CALL xios_send_field("stream_reservoir_restart",stream_reservoir_r) 2133 2302 ENDIF 2134 2303 2135 2304 CALL restput_p (rest_id, 'riverflow', nbp_glo, 1, 1, kjit, riverflow_mean, 'scatter', nbp_glo, index_g) 2136 2305 CALL restput_p (rest_id, 'coastalflow', nbp_glo, 1, 1, kjit, coastalflow_mean, 'scatter', nbp_glo, index_g) 2137 2306 CALL restput_p (rest_id, 'lakeinflow', nbp_glo, 1, 1, kjit, lakeinflow_mean, 'scatter', nbp_glo, index_g) 2138 2307 CALL routing_flow_finalize_mean(kjit, rest_id) 2308 2309 CALL routing_flow_finalize_diag(kjit, rest_id) 2139 2310 2140 2311 CALL xios_orchidee_change_context("orchidee_routing_out") … … 2144 2315 END SUBROUTINE routing_flow_finalize 2145 2316 2317 SUBROUTINE routing_flow_finalize_diag(kjit, rest_id) 2318 USE ioipsl 2319 IMPLICIT NONE 2320 INTEGER, INTENT(IN) :: kjit 2321 INTEGER, INTENT(IN) :: rest_id 2322 2323 CALL restput_p (rest_id, 'hydrographs', nbp_glo, 1, 1, kjit, hydrographs_diag, 'scatter', nbp_glo, index_g) 2324 2325 END SUBROUTINE routing_flow_finalize_diag 2146 2326 2147 2327 !! ================================================================================================================================
Note: See TracChangeset
for help on using the changeset viewer.
