Changeset 7709 for branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90
- Timestamp:
- 2022-07-20T11:30:43+02:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/ORCHIDEE_2_2/ORCHIDEE/src_sechiba/routing.f90
r7576 r7709 12 12 !!\n DESCRIPTION: None 13 13 !! 14 !! RECENT CHANGE(S): None15 !! 14 !! RECENT CHANGE(S): July 2022: New irrigation scheme. Here new irrigation offer with more information from maps, and 15 !! calculation of water withdrawal. 16 16 !! REFERENCE(S) : 17 17 !! … … 226 226 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: reinfiltration_mean !! Mean water flow which returns to the grid box (kg/m^2/dt) 227 227 !$OMP THREADPRIVATE(reinfiltration_mean) 228 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigdeficit_mean !! Mean irrigation deficit. 229 !! This is between irrigation and water requirement (kg/m^2/dt) 230 !$OMP THREADPRIVATE(irrigdeficit_mean) 228 231 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigation_mean !! Mean irrigation flux. 229 232 !! This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt) 230 233 !$OMP THREADPRIVATE(irrigation_mean) 234 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_gw_source_mean !! Mean groundwater irrigation flux. 235 !! This is the water taken from the GW reservoir only (kg/m^2/dt) 236 !$OMP THREADPRIVATE(irrig_gw_source_mean) 237 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_fast_source_mean !! Mean irrigation flux from fast. 238 !! This is the water taken from the fast reservoir only (kg/m^2/dt) 239 !$OMP THREADPRIVATE(irrig_fast_source_mean) 240 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrig_str_source_mean !! Mean streamflow irrigation flux. 241 !! This is the water taken from the streamflow reservoir only (kg/m^2/dt) 242 !$OMP THREADPRIVATE(irrig_str_source_mean) 243 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: irrigadduct_mean !! Irrigation that comes from adduction. It includes water from basins inside the grid cell 244 ! and water from nearby grid cells 245 !$OMP THREADPRIVATE(irrigadduct_mean) 231 246 REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:) :: riverflow_mean !! Mean Outflow of the major rivers. 232 247 !! The flux will be located on the continental grid but this should be a coastal point (kg/dt) … … 291 306 neighbours, resolution, contfrac, stempdiag, & 292 307 returnflow, reinfiltration, irrigation, riverflow, & 293 coastalflow, flood_frac, flood_res ) 294 308 coastalflow, flood_frac, flood_res, soiltile, irrig_frac, & 309 veget_max, irrigated_next) ! 310 295 311 IMPLICIT NONE 296 312 … … 309 325 REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) 310 326 REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 327 REAL(r_std), INTENT(in) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless) 328 REAL(r_std), INTENT(in) :: veget_max(nbpt,nvm) !! Maximal fraction of vegetation (unitless;0-1) ! 329 REAL(r_std), INTENT(in) :: irrigated_next (nbpt) !! Dynamic irrig. area, calculated in slowproc and passed to routing! 330 REAL(r_std), INTENT(in) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 311 331 312 332 !! 0.2 Output variables … … 323 343 !! 0.3 Local variables 324 344 REAL(r_std), DIMENSION(nbp_glo):: mask_coast_glo !! Mask with coastal gridcells on global grid (1/0) 325 LOGICAL :: init_irrig!! Logical to initialize the irrigation (true/false)345 !LOGICAL :: init_irrig !! Logical to initialize the irrigation (true/false) 326 346 LOGICAL :: init_flood !! Logical to initialize the floodplains (true/false) 327 347 LOGICAL :: init_swamp !! Logical to initialize the swamps (true/false) … … 385 405 !! Initialisation of flags for irrigated land, flood plains and swamps 386 406 ! 387 init_irrig = .FALSE.388 407 IF ( do_irrigation ) THEN 389 IF (COUNT(irrigated .GE. undef_sechiba-1) > 0) init_irrig = .TRUE.408 irrigated(:) = irrigated_next(:) 390 409 END IF 391 410 … … 399 418 IF (COUNT(swamp .GE. undef_sechiba-1) > 0 ) init_swamp = .TRUE. 400 419 END IF 401 420 402 421 !! If we have irrigated land, flood plains or swamps then we need to interpolate the 0.5 degree 403 422 !! base data set to the resolution of the model. 404 405 IF ( init_ irrig .OR. init_flood .OR. init_swamp ) THEN406 CALL routing_ irrigmap(nbpt, index, lalo, neighbours, resolution, &407 contfrac, init_ irrig, irrigated, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)423 !IF ( init_irrig .OR. init_flood .OR. init_swamp ) THEN 424 IF ( init_flood .OR. init_swamp ) THEN 425 CALL routing_floodmap(nbpt, index, lalo, neighbours, resolution, & 426 contfrac, init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) 408 427 ENDIF 409 410 IF ( do_irrigation ) THEN 411 CALL xios_orchidee_send_field("irrigmap",irrigated) 412 413 IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) 414 IF ( .NOT. almaoutput ) THEN 415 CALL histwrite_p(hist_id, 'irrigmap', 1, irrigated, nbpt, index) 416 ELSE 417 CALL histwrite_p(hist_id, 'IrrigationMap', 1, irrigated, nbpt, index) 418 ENDIF 419 IF ( hist2_id > 0 ) THEN 420 IF ( .NOT. almaoutput ) THEN 421 CALL histwrite_p(hist2_id, 'irrigmap', 1, irrigated, nbpt, index) 422 ELSE 423 CALL histwrite_p(hist2_id, 'IrrigationMap', 1, irrigated, nbpt, index) 424 ENDIF 425 ENDIF 428 429 IF ( do_irrigation ) THEN 430 431 IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrigated : ', MINVAL(irrigated), MAXVAL(irrigated) 432 IF (printlev >= 3) WRITE(numout,*) 'Verification : range of irrig_frac : ', MINVAL(irrig_frac), MAXVAL(irrig_frac) 426 433 ENDIF 427 428 IF ( do_floodplains ) THEN 429 CALL xios_orchidee_send_field("floodmap",floodplains) 430 431 IF (printlev>=3) WRITE(numout,*) 'Verification : range of floodplains : ', MINVAL(floodplains), MAXVAL(floodplains) 432 IF ( .NOT. almaoutput ) THEN 433 CALL histwrite_p(hist_id, 'floodmap', 1, floodplains, nbpt, index) 434 ELSE 435 CALL histwrite_p(hist_id, 'FloodplainsMap', 1, floodplains, nbpt, index) 436 ENDIF 437 IF ( hist2_id > 0 ) THEN 438 IF ( .NOT. almaoutput ) THEN 439 CALL histwrite_p(hist2_id, 'floodmap', 1, floodplains, nbpt, index) 440 ELSE 441 CALL histwrite_p(hist2_id, 'FloodplainsMap', 1, floodplains, nbpt, index) 442 ENDIF 443 ENDIF 444 ENDIF 445 434 446 435 IF ( doswamps ) THEN 447 436 CALL xios_orchidee_send_field("swampmap",swamp) … … 558 547 & lalo, neighbours, resolution, contfrac, totfrac_nobio, veget_max, floodout, runoff, & 559 548 & drainage, transpot, precip_rain, humrel, k_litt, flood_frac, flood_res, & 560 & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id) 549 & stempdiag, reinf_slope, returnflow, reinfiltration, irrigation, riverflow, coastalflow, rest_id, hist_id, hist2_id,& 550 soiltile, root_deficit, irrigated_next, irrig_frac, fraction_aeirrig_sw) ! 561 551 562 552 IMPLICIT NONE … … 584 574 REAL(r_std), INTENT(in) :: stempdiag(nbpt,nslm) !! Diagnostic soil temperature profile 585 575 REAL(r_std), INTENT(in) :: reinf_slope(nbpt) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) 576 REAL(r_std), INTENT(in) :: root_deficit(nbpt) !! soil water deficit 577 REAL(r_std), INTENT(in) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless) 578 REAL(r_std), INTENT(in) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 579 REAL(r_std), INTENT(in) :: irrigated_next (nbpt)!! Dynamic irrig. area, calculated in slowproc and passed to routing 580 REAL(r_std), INTENT(in) :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac 586 581 587 582 !! 0.2 Output variables … … 670 665 ! the growing seasons we will give more weight to these areas. 671 666 ! 672 DO jv=2,nvm 673 DO ig=1,nbpt 674 humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing 675 vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing 676 ENDDO 677 ENDDO 667 ! New irrigation scheme uses mean of vegtot with jv 1 to nvm 668 ! Old scheme keeps jv 2 to nvm, even if possibly an error 669 IF ( .NOT. old_irrig_scheme ) THEN 670 DO jv=1,nvm 671 DO ig=1,nbpt 672 humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing 673 vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing 674 ENDDO 675 ENDDO 676 ELSE 677 DO jv=2,nvm 678 DO ig=1,nbpt 679 humrel_mean(ig) = humrel_mean(ig) + humrel(ig,jv)*veget_max(ig,jv)*dt_sechiba/dt_routing 680 vegtot_mean(ig) = vegtot_mean(ig) + veget_max(ig,jv)*dt_sechiba/dt_routing 681 ENDDO 682 ENDDO 683 ENDIF 684 ! Here updates irrigmap to irrigated_next from slowproc, every timestep 685 !irrigated_next is updated in slowproc when time comes 686 !irrig_frac was also updated in slowproc, here used as input variable 687 688 IF ( do_irrigation .AND. irrig_map_dynamic_flag ) THEN 689 irrigated(:) = irrigated_next(:) 690 ENDIF 678 691 ! 679 692 time_counter = time_counter + dt_sechiba … … 685 698 !! Computes the transport of water in the various reservoirs 686 699 ! 687 CALL routing_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, &700 CALL routing_flow(nbpt, dt_routing, lalo, floodout_mean, runoff_mean, drainage_mean, root_deficit, soiltile, & 688 701 & vegtot_mean, totnobio_mean, transpot_mean, precip_mean, humrel_mean, k_litt_mean, floodtemp, reinf_slope, & 689 & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrigation_mean, riverflow_mean, & 702 & lakeinflow_mean, returnflow_mean, reinfiltration_mean, irrig_frac, irrigation_mean, irrigdeficit_mean, & 703 & irrigadduct_mean, irrig_gw_source_mean, & ! 704 & irrig_fast_source_mean, irrig_str_source_mean, riverflow_mean, & ! 690 705 & coastalflow_mean, hydrographs, slowflow_diag, flood_frac, flood_res, & 691 & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag )706 & netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, fraction_aeirrig_sw) 692 707 ! 693 708 !! Responsible for storing the water in lakes … … 718 733 reinfiltration_mean(:) = reinfiltration_mean(:)/dt_routing*dt_sechiba 719 734 irrigation_mean(:) = irrigation_mean(:)/dt_routing*dt_sechiba 735 irrigdeficit_mean(:) = irrigdeficit_mean/dt_routing*dt_sechiba ! 736 irrigadduct_mean(:) = irrigadduct_mean(:)/dt_routing*dt_sechiba! 737 irrig_gw_source_mean(:) = irrig_gw_source_mean(:)/dt_routing*dt_sechiba ! 738 irrig_fast_source_mean(:) = irrig_fast_source_mean(:)/dt_routing*dt_sechiba ! 739 irrig_str_source_mean(:) = irrig_str_source_mean(:)/dt_routing*dt_sechiba ! 720 740 irrig_netereq(:) = irrig_netereq(:)/dt_routing*dt_sechiba 721 741 … … 766 786 ! Water fluxes converted from kg/m^2/dt_sechiba into kg/m^2/s 767 787 CALL xios_orchidee_send_field("irrigation",irrigation/dt_sechiba) 788 CALL xios_orchidee_send_field("irrig_deficit",irrigdeficit_mean/dt_sechiba)! 789 CALL xios_orchidee_send_field("irrig_adduct",irrigadduct_mean/dt_sechiba)! 790 CALL xios_orchidee_send_field("irrig_gw_source",irrig_gw_source_mean/dt_sechiba) ! 791 CALL xios_orchidee_send_field("irrig_fast_source",irrig_fast_source_mean/dt_sechiba) ! 792 CALL xios_orchidee_send_field("irrig_str_source",irrig_str_source_mean/dt_sechiba) ! 768 793 CALL xios_orchidee_send_field("netirrig",irrig_netereq/dt_sechiba) 769 794 CALL xios_orchidee_send_field("riversret",returnflow/dt_sechiba) … … 879 904 880 905 !! 0.2 Local variables 881 REAL(r_std), DIMENSION(1) :: tmp_day 906 REAL(r_std), DIMENSION(1) :: tmp_day 882 907 883 908 !_ ================================================================================================================================ … … 931 956 932 957 IF ( do_irrigation ) THEN 933 CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter', nbp_glo, index_g)958 !CALL restput_p (rest_id, 'irrigated', nbp_glo, 1, 1, kjit, irrigated, 'scatter', nbp_glo, index_g) 934 959 CALL restput_p (rest_id, 'irrigation', nbp_glo, 1, 1, kjit, irrigation_mean, 'scatter', nbp_glo, index_g) 960 CALL restput_p (rest_id, 'irrigdeficit', nbp_glo, 1, 1, kjit, irrigdeficit_mean, 'scatter', nbp_glo, index_g) 961 CALL restput_p (rest_id, 'irrigadduct', nbp_glo, 1, 1, kjit, irrigadduct_mean, 'scatter', nbp_glo, index_g) 962 CALL restput_p (rest_id, 'irrig_gw_source', nbp_glo, 1, 1, kjit, irrig_gw_source_mean, 'scatter', nbp_glo, index_g) ! 963 CALL restput_p (rest_id, 'irrig_fast_source', nbp_glo, 1, 1, kjit, irrig_fast_source_mean, 'scatter', nbp_glo, index_g) ! 964 CALL restput_p (rest_id, 'irrig_str_source', nbp_glo, 1, 1, kjit, irrig_str_source_mean, 'scatter', nbp_glo, index_g) ! 935 965 ENDIF 936 966 … … 1022 1052 ! 1023 1053 !Config Key = DO_FLOODINFILT 1024 !Config Desc = Should floodplains reinfiltrate into the soil 1054 !Config Desc = Should floodplains reinfiltrate into the soil 1025 1055 !Config If = RIVER_ROUTING 1026 1056 !Config Def = n 1027 1057 !Config Help = This parameters allows the user to ask the model 1028 !Config to take into account the flood plains reinfiltration 1029 !Config into the soil moisture. It then can go 1058 !Config to take into account the flood plains reinfiltration 1059 !Config into the soil moisture. It then can go 1030 1060 !Config back to the slow and fast reservoirs 1031 1061 !Config Units = [FLAG] … … 1039 1069 !Config Def = n 1040 1070 !Config Help = This parameters allows the user to ask the model 1041 !Config to take into account the swamps and return 1042 !Config the water into the bottom of the soil. It then can go 1043 !Config back to the atmopshere. This tried to simulate 1071 !Config to take into account the swamps and return 1072 !Config the water into the bottom of the soil. It then can go 1073 !Config back to the atmopshere. This tried to simulate 1044 1074 !Config internal deltas of rivers. 1045 1075 !Config Units = [FLAG] … … 1049 1079 ! 1050 1080 !Config Key = DO_PONDS 1051 !Config Desc = Should we include ponds 1081 !Config Desc = Should we include ponds 1052 1082 !Config If = RIVER_ROUTING 1053 1083 !Config Def = n 1054 1084 !Config Help = This parameters allows the user to ask the model 1055 !Config to take into account the ponds and return 1056 !Config the water into the soil moisture. It then can go 1057 !Config back to the atmopshere. This tried to simulate 1085 !Config to take into account the ponds and return 1086 !Config the water into the soil moisture. It then can go 1087 !Config back to the atmopshere. This tried to simulate 1058 1088 !Config little ponds especially in West Africa. 1059 1089 !Config Units = [FLAG] … … 1064 1094 1065 1095 !Config Key = SLOW_TCST 1066 !Config Desc = Time constant for the slow reservoir 1067 !Config If = RIVER_ROUTING 1068 !Config Def = 25.0 1069 !Config Help = This parameters allows the user to fix the 1096 !Config Desc = Time constant for the slow reservoir 1097 !Config If = RIVER_ROUTING 1098 !Config Def = 25.0 1099 !Config Help = This parameters allows the user to fix the 1070 1100 !Config time constant (in days) of the slow reservoir 1071 !Config in order to get better river flows for 1101 !Config in order to get better river flows for 1072 1102 !Config particular regions. 1073 1103 !Config Units = [days] … … 1078 1108 !> during the 1 degree NCEP Corrected by Cru (NCC) resolution simulations (Ngo-Duc et al., 2005, Ngo-Duc et al., 2006) and 1079 1109 !> generalized for all the basins of the world. The "slow reservoir" and the "fast reservoir" 1080 !> have the highest value in order to simulate the groundwater. 1110 !> have the highest value in order to simulate the groundwater. 1081 1111 !> The "stream reservoir", which represents all the water of the stream, has the lowest value. 1082 1112 !> Those figures are the same for all the basins of the world. … … 1087 1117 ! 1088 1118 !Config Key = FAST_TCST 1089 !Config Desc = Time constant for the fast reservoir 1090 !Config If = RIVER_ROUTING 1091 !Config Def = 3.0 1092 !Config Help = This parameters allows the user to fix the 1119 !Config Desc = Time constant for the fast reservoir 1120 !Config If = RIVER_ROUTING 1121 !Config Def = 3.0 1122 !Config Help = This parameters allows the user to fix the 1093 1123 !Config time constant (in days) of the fast reservoir 1094 !Config in order to get better river flows for 1124 !Config in order to get better river flows for 1095 1125 !Config particular regions. 1096 1126 !Config Units = [days] 1097 1127 CALL getin_p('FAST_TCST', fast_tcst) 1098 1128 1099 1129 !Config Key = STREAM_TCST 1100 !Config Desc = Time constant for the stream reservoir 1130 !Config Desc = Time constant for the stream reservoir 1101 1131 !Config If = RIVER_ROUTING 1102 1132 !Config Def = 0.24 1103 !Config Help = This parameters allows the user to fix the 1133 !Config Help = This parameters allows the user to fix the 1104 1134 !Config time constant (in days) of the stream reservoir 1105 !Config in order to get better river flows for 1135 !Config in order to get better river flows for 1106 1136 !Config particular regions. 1107 1137 !Config Units = [days] 1108 1138 CALL getin_p('STREAM_TCST', stream_tcst) 1109 1139 1110 1140 !Config Key = FLOOD_TCST 1111 !Config Desc = Time constant for the flood reservoir 1141 !Config Desc = Time constant for the flood reservoir 1112 1142 !Config If = RIVER_ROUTING 1113 1143 !Config Def = 4.0 1114 !Config Help = This parameters allows the user to fix the 1144 !Config Help = This parameters allows the user to fix the 1115 1145 !Config time constant (in days) of the flood reservoir 1116 !Config in order to get better river flows for 1146 !Config in order to get better river flows for 1117 1147 !Config particular regions. 1118 1148 !Config Units = [days] 1119 1149 CALL getin_p('FLOOD_TCST', flood_tcst) 1120 1150 1121 1151 !Config Key = SWAMP_CST 1122 !Config Desc = Fraction of the river that flows back to swamps 1152 !Config Desc = Fraction of the river that flows back to swamps 1123 1153 !Config If = RIVER_ROUTING 1124 1154 !Config Def = 0.2 1125 !Config Help = This parameters allows the user to fix the 1155 !Config Help = This parameters allows the user to fix the 1126 1156 !Config fraction of the river transport 1127 1157 !Config that flows to swamps 1128 1158 !Config Units = [-] 1129 1159 CALL getin_p('SWAMP_CST', swamp_cst) 1130 1160 1131 1161 !Config Key = FLOOD_BETA 1132 !Config Desc = Parameter to fix the shape of the floodplain 1162 !Config Desc = Parameter to fix the shape of the floodplain 1133 1163 !Config If = RIVER_ROUTING 1134 1164 !Config Def = 2.0 1135 1165 !Config Help = Parameter to fix the shape of the floodplain 1136 1166 !Config (>1 for convex edges, <1 for concave edges) 1137 !Config Units = [-] 1167 !Config Units = [-] 1138 1168 CALL getin_p("FLOOD_BETA", beta) 1139 1169 ! … … 1142 1172 !Config If = RIVER_ROUTING 1143 1173 !Config Def = 0.5 1144 !Config Help = 1145 !Config Units = [-] 1146 CALL getin_p("POND_BETAP", betap) 1174 !Config Help = 1175 !Config Units = [-] 1176 CALL getin_p("POND_BETAP", betap) 1147 1177 ! 1148 1178 !Config Key = FLOOD_CRI … … 1150 1180 !Config If = DO_FLOODPLAINS or DO_PONDS 1151 1181 !Config Def = 2000. 1152 !Config Help = 1153 !Config Units = [mm] 1182 !Config Help = 1183 !Config Units = [mm] 1154 1184 CALL getin_p("FLOOD_CRI", floodcri) 1155 1185 ! … … 1158 1188 !Config If = DO_FLOODPLAINS or DO_PONDS 1159 1189 !Config Def = 2000. 1160 !Config Help = 1161 !Config Units = [mm] 1190 !Config Help = 1191 !Config Units = [mm] 1162 1192 CALL getin_p("POND_CRI", pondcri) 1163 1193 … … 1166 1196 !Config If = RIVER_ROUTING 1167 1197 !Config Def = 7000 1168 !Config Help = 1169 !Config Units = [kg/m2(routing area)] 1198 !Config Help = 1199 !Config Units = [kg/m2(routing area)] 1170 1200 max_lake_reservoir = 7000 1171 1201 CALL getin_p("MAX_LAKE_RESERVOIR", max_lake_reservoir) … … 1205 1235 ELSE 1206 1236 ! Take the value from restart file 1207 time_counter = tmp_day(1) 1237 time_counter = tmp_day(1) 1208 1238 ENDIF 1209 1239 ENDIF 1210 1240 CALL bcast(time_counter) 1211 1241 1212 1242 1213 1243 ALLOCATE (routing_area_loc(nbpt,nbasmax), stat=ier) 1214 1244 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for routing_area_loc','','') … … 1365 1395 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_height, "gather", nbp_glo, index_g) 1366 1396 CALL setvar_p (flood_height, val_exp, 'NO_KEYWORD', zero) 1367 1397 1368 1398 ALLOCATE (pond_frac(nbpt), stat=ier) 1369 1399 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for pond_frac','','') … … 1373 1403 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., pond_frac, "gather", nbp_glo, index_g) 1374 1404 CALL setvar_p (pond_frac, val_exp, 'NO_KEYWORD', zero) 1375 1405 1376 1406 var_name = 'flood_frac' 1377 1407 CALL ioconf_setatt_p('UNITS', '-') … … 1379 1409 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., flood_frac, "gather", nbp_glo, index_g) 1380 1410 CALL setvar_p (flood_frac, val_exp, 'NO_KEYWORD', zero) 1381 1411 1382 1412 var_name = 'flood_res' 1383 1413 CALL ioconf_setatt_p('UNITS','mm') … … 1386 1416 CALL setvar_p (flood_res, val_exp, 'NO_KEYWORD', zero) 1387 1417 ! flood_res = zero 1388 1418 1389 1419 ALLOCATE (lake_reservoir(nbpt), stat=ier) 1390 1420 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for lake_reservoir','','') … … 1394 1424 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lake_reservoir, "gather", nbp_glo, index_g) 1395 1425 CALL setvar_p (lake_reservoir, val_exp, 'NO_KEYWORD', zero) 1396 1426 1397 1427 ALLOCATE (pond_reservoir(nbpt), stat=ier) 1398 1428 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for pond_reservoir','','') … … 1404 1434 ! 1405 1435 ! Map of irrigated areas 1406 ! 1407 IF ( do_irrigation ) THEN 1408 ALLOCATE (irrigated(nbpt), stat=ier) 1409 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigated','','') 1410 var_name = 'irrigated' 1411 CALL ioconf_setatt_p('UNITS', 'm^2') 1412 CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area') 1413 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigated, "gather", nbp_glo, index_g) 1414 CALL setvar_p (irrigated, val_exp, 'NO_KEYWORD', undef_sechiba) 1415 ENDIF 1416 1436 ! irrigated equal to irrigated_next from slowproc. Here, we initialize to zero. 1437 ! Values from slowproc given in routing_main 1438 ALLOCATE (irrigated(nbpt), stat=ier) 1439 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigated','','') 1440 var_name = 'irrigated' 1441 CALL ioconf_setatt_p('UNITS', 'm^2') 1442 CALL ioconf_setatt_p('LONG_NAME','Surface of irrigated area') 1443 irrigated(:) = zero 1444 1417 1445 IF ( do_floodplains ) THEN 1418 1446 ALLOCATE (floodplains(nbpt), stat=ier) … … 1443 1471 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., lakeinflow_mean, "gather", nbp_glo, index_g) 1444 1472 CALL setvar_p (lakeinflow_mean, val_exp, 'NO_KEYWORD', zero) 1445 1473 1446 1474 ALLOCATE (returnflow_mean(nbpt), stat=ier) 1447 1475 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for returnflow_mean','','') … … 1452 1480 CALL setvar_p (returnflow_mean, val_exp, 'NO_KEYWORD', zero) 1453 1481 returnflow(:) = returnflow_mean(:) 1454 1482 1455 1483 ALLOCATE (reinfiltration_mean(nbpt), stat=ier) 1456 1484 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for reinfiltration_mean','','') … … 1461 1489 CALL setvar_p (reinfiltration_mean, val_exp, 'NO_KEYWORD', zero) 1462 1490 reinfiltration(:) = reinfiltration_mean(:) 1463 1491 1464 1492 ALLOCATE (irrigation_mean(nbpt), stat=ier) 1465 1493 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigation_mean','','') 1494 ALLOCATE (irrigdeficit_mean(nbpt), stat=ier) 1495 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigdeficit_mean','','') 1496 ALLOCATE (irrigadduct_mean(nbpt), stat=ier) 1497 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrigadduct_mean','','') 1466 1498 ALLOCATE (irrig_netereq(nbpt), stat=ier) 1467 1499 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_netereq','','') 1500 ALLOCATE (irrig_gw_source_mean(nbpt), stat=ier) ! 1501 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_gw_source_mean','','')! 1502 ALLOCATE (irrig_fast_source_mean(nbpt), stat=ier) ! 1503 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_fast_source_mean','','')! 1504 ALLOCATE (irrig_str_source_mean(nbpt), stat=ier) ! 1505 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for irrig_str_source_mean','','')! 1468 1506 irrig_netereq(:) = zero 1469 1507 1470 1508 IF ( do_irrigation ) THEN 1471 1509 var_name = 'irrigation' … … 1474 1512 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigation_mean, "gather", nbp_glo, index_g) 1475 1513 CALL setvar_p (irrigation_mean, val_exp, 'NO_KEYWORD', zero) 1514 var_name = 'irrigdeficit' 1515 CALL ioconf_setatt_p('UNITS', 'Kg/dt') 1516 CALL ioconf_setatt_p('LONG_NAME','Irrigation deficit') 1517 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigdeficit_mean, "gather", nbp_glo, index_g) 1518 CALL setvar_p (irrigdeficit_mean, val_exp, 'NO_KEYWORD', zero) 1519 var_name = 'irrigadduct' 1520 CALL ioconf_setatt_p('UNITS', 'Kg/dt') 1521 CALL ioconf_setatt_p('LONG_NAME','Artificial irrigation flux from adduction') 1522 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrigadduct_mean, "gather", nbp_glo, index_g) 1523 CALL setvar_p (irrigadduct_mean, val_exp, 'NO_KEYWORD', zero) 1524 var_name = 'irrig_gw_source'! 1525 CALL ioconf_setatt_p('UNITS', 'Kg/dt') 1526 CALL ioconf_setatt_p('LONG_NAME','Irrigation from GW reservoir') 1527 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_gw_source_mean, "gather", nbp_glo, index_g) 1528 CALL setvar_p (irrig_gw_source_mean, val_exp, 'NO_KEYWORD', zero) 1529 var_name = 'irrig_fast_source'! 1530 CALL ioconf_setatt_p('UNITS', 'Kg/dt') 1531 CALL ioconf_setatt_p('LONG_NAME','Irrigation from fast reservoir') 1532 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_fast_source_mean, "gather", nbp_glo, index_g) 1533 CALL setvar_p (irrig_fast_source_mean, val_exp, 'NO_KEYWORD', zero) 1534 var_name = 'irrig_str_source'! 1535 CALL ioconf_setatt_p('UNITS', 'Kg/dt') 1536 CALL ioconf_setatt_p('LONG_NAME','Irrigation from stream reservoir') 1537 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., irrig_str_source_mean, "gather", nbp_glo, index_g) 1538 CALL setvar_p (irrig_str_source_mean, val_exp, 'NO_KEYWORD', zero) 1476 1539 ELSE 1477 1540 irrigation_mean(:) = zero 1541 irrig_gw_source_mean(:) = zero 1542 irrig_fast_source_mean(:) = zero 1543 irrig_str_source_mean(:) = zero 1544 irrigdeficit_mean(:) = zero 1545 irrigadduct_mean(:) = zero 1478 1546 ENDIF 1479 irrigation(:) = irrigation_mean(:) 1480 1547 irrigation(:) = irrigation_mean(:) 1548 1481 1549 ALLOCATE (riverflow_mean(nbpt), stat=ier) 1482 1550 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for riverflow_mean','','') … … 1487 1555 CALL setvar_p (riverflow_mean, val_exp, 'NO_KEYWORD', zero) 1488 1556 riverflow(:) = riverflow_mean(:) 1489 1557 1490 1558 ALLOCATE (coastalflow_mean(nbpt), stat=ier) 1491 1559 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for coastalflow_mean','','') … … 1496 1564 CALL setvar_p (coastalflow_mean, val_exp, 'NO_KEYWORD', zero) 1497 1565 coastalflow(:) = coastalflow_mean(:) 1498 1566 1499 1567 ! Locate it at the 2m level 1500 1568 ipn = MINLOC(ABS(diaglev-2)) … … 1503 1571 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for floodtemp','','') 1504 1572 floodtemp(:) = stempdiag(:,floodtemp_lev) 1505 1573 1506 1574 ALLOCATE(hydrographs(nbpt), stat=ier) 1507 1575 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for hydrographs','','') … … 1511 1579 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., hydrographs, "gather", nbp_glo, index_g) 1512 1580 CALL setvar_p (hydrographs, val_exp, 'NO_KEYWORD', zero) 1513 1581 1514 1582 ALLOCATE(slowflow_diag(nbpt), stat=ier) 1515 1583 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for slowflow_diag','','') … … 1526 1594 & pond_diag(nbpt), lake_diag(nbpt), stat=ier) 1527 1595 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for fast_diag,..','','') 1528 1596 1529 1597 fast_diag(:) = zero 1530 1598 slow_diag(:) = zero … … 1533 1601 pond_diag(:) = zero 1534 1602 lake_diag(:) = zero 1535 1603 1536 1604 DO ig=1,nbpt 1537 1605 totarea = zero … … 1566 1634 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., floodout_mean, "gather", nbp_glo, index_g) 1567 1635 CALL setvar_p (floodout_mean, val_exp, 'NO_KEYWORD', zero) 1568 1636 1569 1637 ALLOCATE (runoff_mean(nbpt), stat=ier) 1570 1638 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for runoff_mean','','') … … 1574 1642 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., runoff_mean, "gather", nbp_glo, index_g) 1575 1643 CALL setvar_p (runoff_mean, val_exp, 'NO_KEYWORD', zero) 1576 1644 1577 1645 ALLOCATE(drainage_mean(nbpt), stat=ier) 1578 1646 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for drainage_mean','','') … … 1582 1650 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., drainage_mean, "gather", nbp_glo, index_g) 1583 1651 CALL setvar_p (drainage_mean, val_exp, 'NO_KEYWORD', zero) 1584 1652 1585 1653 ALLOCATE(transpot_mean(nbpt), stat=ier) 1586 1654 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for transpot_mean','','') … … 1598 1666 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., precip_mean, "gather", nbp_glo, index_g) 1599 1667 CALL setvar_p (precip_mean, val_exp, 'NO_KEYWORD', zero) 1600 1668 1601 1669 ALLOCATE(humrel_mean(nbpt), stat=ier) 1602 1670 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for humrel_mean','','') … … 1606 1674 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., humrel_mean, "gather", nbp_glo, index_g) 1607 1675 CALL setvar_p (humrel_mean, val_exp, 'NO_KEYWORD', un) 1608 1676 1609 1677 ALLOCATE(k_litt_mean(nbpt), stat=ier) 1610 1678 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for k_litt_mean','','') … … 1614 1682 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., k_litt_mean, "gather", nbp_glo, index_g) 1615 1683 CALL setvar_p (k_litt_mean, val_exp, 'NO_KEYWORD', zero) 1616 1684 1617 1685 ALLOCATE(totnobio_mean(nbpt), stat=ier) 1618 1686 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for totnobio_mean','','') … … 1622 1690 CALL restget_p (rest_id, var_name, nbp_glo, 1, 1, kjit, .TRUE., totnobio_mean, "gather", nbp_glo, index_g) 1623 1691 CALL setvar_p (totnobio_mean, val_exp, 'NO_KEYWORD', zero) 1624 1692 1625 1693 ALLOCATE(vegtot_mean(nbpt), stat=ier) 1626 1694 IF (ier /= 0) CALL ipslerr_p(3,'routing_init','Pb in allocate for vegtot_mean','','') … … 1693 1761 IF (ALLOCATED(hydrodiag_loc)) DEALLOCATE(hydrodiag_loc) 1694 1762 IF (ALLOCATED(hydrodiag_glo)) DEALLOCATE(hydrodiag_glo) 1695 IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc) 1763 IF (ALLOCATED(hydroupbasin_loc)) DEALLOCATE(hydroupbasin_loc) 1696 1764 IF (ALLOCATED(hydroupbasin_glo)) DEALLOCATE(hydroupbasin_glo) 1697 1765 IF (ALLOCATED(hydrographs)) DEALLOCATE(hydrographs) 1698 1766 IF (ALLOCATED(slowflow_diag)) DEALLOCATE(slowflow_diag) 1699 1767 IF (ALLOCATED(irrigation_mean)) DEALLOCATE(irrigation_mean) 1768 IF (ALLOCATED(irrigdeficit_mean)) DEALLOCATE(irrigdeficit_mean)! 1769 IF (ALLOCATED(irrigadduct_mean)) DEALLOCATE(irrigadduct_mean)! 1770 IF (ALLOCATED(irrig_gw_source_mean)) DEALLOCATE(irrig_gw_source_mean) ! 1771 IF (ALLOCATED(irrig_fast_source_mean)) DEALLOCATE(irrig_fast_source_mean) ! 1772 IF (ALLOCATED(irrig_str_source_mean)) DEALLOCATE(irrig_str_source_mean) ! 1700 1773 IF (ALLOCATED(irrigated)) DEALLOCATE(irrigated) 1701 1774 IF (ALLOCATED(floodplains)) DEALLOCATE(floodplains) … … 1718 1791 !! 1719 1792 !! DESCRIPTION (definitions, functional, design, flags) : 1720 !! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an 1721 !! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes 1722 !! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of 1793 !! This will first compute the amount of water which flows out of each of the 3 reservoirs using the assumption of an 1794 !! exponential decrease of water in the reservoir (see Hagemann S and Dumenil L. (1998)). Then we compute the fluxes 1795 !! for floodplains and ponds. All this will then be used in order to update each of the basins : taking water out of 1723 1796 !! the up-stream basin and adding it to the down-stream one. 1724 !! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once 1797 !! As this step happens globaly we have to stop the parallel processing in order to exchange the information. Once 1725 1798 !! all reservoirs are updated we deal with irrigation. The final step is to compute diagnostic fluxes. Among them 1726 1799 !! the hydrographs of the largest rivers we have chosen to monitor. … … 1778 1851 !_ ================================================================================================================================ 1779 1852 1780 SUBROUTINE routing_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, &1853 SUBROUTINE routing_flow(nbpt, dt_routing, lalo, floodout, runoff, drainage, root_deficit, soiltile, & 1781 1854 & vegtot, totnobio, transpot_mean, precip, humrel, k_litt, floodtemp, reinf_slope, & 1782 & lakeinflow, returnflow, reinfiltration, irrigation, riverflow, & 1855 & lakeinflow, returnflow, reinfiltration, irrig_frac, irrigation, irrigdeficit, & 1856 & irrigadduct, irrig_gw_source, & 1857 & irrig_fast_source, irrig_str_source, riverflow, & ! 1783 1858 & coastalflow, hydrographs, slowflow_diag, flood_frac, flood_res, & 1784 netflow_stream_diag, netflow_fast_diag, netflow_slow_diag )1859 netflow_stream_diag, netflow_fast_diag, netflow_slow_diag, fraction_aeirrig_sw) ! 1785 1860 ! 1786 1861 IMPLICIT NONE … … 1793 1868 REAL(r_std), INTENT(in) :: floodout(nbpt) !! Grid-point flow out of floodplains (kg/m^2/dt) 1794 1869 REAL(r_std), INTENT(in) :: drainage(nbpt) !! Grid-point drainage (kg/m^2/dt) 1870 REAL(r_std), INTENT(in) :: root_deficit(nbpt) !! soil water deficit 1795 1871 REAL(r_std), INTENT(in) :: vegtot(nbpt) !! Potentially vegetated fraction (unitless;0-1) 1796 1872 REAL(r_std), INTENT(in) :: totnobio(nbpt) !! Other areas which can not have vegetation … … 1802 1878 REAL(r_std), INTENT(in) :: reinf_slope(nbpt) !! Coefficient which determines the reinfiltration ratio in the grid box due to flat areas (unitless;0-1) 1803 1879 REAL(r_std), INTENT(out) :: lakeinflow(nbpt) !! Water inflow to the lakes (kg/dt) 1880 REAL(r_std), INTENT(in) :: soiltile(nbpt,nstm) !! Fraction of each soil tile within vegtot (0-1, unitless) 1881 REAL(r_std), INTENT(in) :: irrig_frac(nbpt) !! Irrig. fraction interpolated in routing, and saved to pass to slowproc if irrigated_soiltile = .TRUE. 1882 REAL(r_std), INTENT(in) :: fraction_aeirrig_sw(nbpt) !! Fraction of area equipped for irrigation from surface water, of irrig_frac 1883 1804 1884 ! 1805 1885 !! OUTPUT VARIABLES … … 1808 1888 REAL(r_std), INTENT(out) :: reinfiltration(nbpt) !! Water flow from ponds and floodplains which returns to the grid box (kg/m^2/dt) 1809 1889 REAL(r_std), INTENT(out) :: irrigation(nbpt) !! Irrigation flux. This is the water taken from the reservoirs and beeing put into the upper layers of the soil (kg/m^2/dt_routing) 1890 REAL(r_std), INTENT(out) :: irrigdeficit(nbpt) !! Irrigation deficit. Difference btw irrig. demand and irrigation 1891 REAL(r_std), INTENT(out) :: irrigadduct(nbpt) !! Irrigation from adducted water. Included in irrigation 1892 1810 1893 REAL(r_std), INTENT(out) :: riverflow(nbpt) !! Outflow of the major rivers. The flux will be located on the continental grid but this should be a coastal point (kg/dt_routing) 1811 1894 REAL(r_std), INTENT(out) :: coastalflow(nbpt) !! Outflow on coastal points by small basins. This is the water which flows in a disperse way into the ocean (kg/dt_routing) … … 1818 1901 REAL(r_std), INTENT(out) :: netflow_fast_diag(nbpt) !! Input - Output flow to fast reservoir 1819 1902 REAL(r_std), INTENT(out) :: netflow_slow_diag(nbpt) !! Input - Output flow to slow reservoir 1903 1904 REAL(r_std), INTENT(out) :: irrig_gw_source(nbpt) !! Irrigation from groundwater per cell 1905 REAL(r_std), INTENT(out) :: irrig_fast_source(nbpt) !! Irrigation from fast reservoir per cell 1906 REAL(r_std), INTENT(out) :: irrig_str_source(nbpt) !! Irrigation from stramflow per cell 1820 1907 ! 1821 1908 !! LOCAL VARIABLES … … 1843 1930 REAL(r_std), DIMENSION(nbpt) :: totarea !! Total area of basin (m^2) 1844 1931 REAL(r_std), DIMENSION(nbpt) :: totflood !! Total amount of water in the floodplains reservoir (kg) 1845 REAL(r_std), DIMENSION(nbasmax) :: pond_excessflow !! 1932 REAL(r_std), DIMENSION(nbasmax) :: pond_excessflow !! 1846 1933 REAL(r_std) :: flow !! Outflow computation for the reservoirs (kg/dt) 1847 1934 REAL(r_std) :: floodindex !! Fraction of grid box area inundated (unitless;0-1) 1848 REAL(r_std) :: pondex !! 1935 REAL(r_std) :: pondex !! 1849 1936 REAL(r_std) :: flood_frac_pot !! Total fraction of the grid box which is flooded at optimum repartition (unitless;0-1) 1850 1937 REAL(r_std) :: stream_tot !! Total water amount in the stream reservoirs (kg) … … 1855 1942 REAL(r_std) :: total_lake_overflow !! Sum of lake_overflow over full grid (kg) 1856 1943 REAL(r_std), DIMENSION(8,nbasmax) :: streams_around !! Stream reservoirs of the neighboring grid boxes (kg) 1857 INTEGER(i_std), DIMENSION(8) :: igrd !! 1858 INTEGER(i_std), DIMENSION(2) :: ff !! 1859 INTEGER(i_std), DIMENSION(1) :: fi !! 1944 INTEGER(i_std), DIMENSION(8) :: igrd !! 1945 INTEGER(i_std), DIMENSION(2) :: ff !! 1946 INTEGER(i_std), DIMENSION(1) :: fi !! 1860 1947 INTEGER(i_std) :: ig, ib, ib2, ig2 !! Indices (unitless) 1861 1948 INTEGER(i_std) :: rtg, rtb, in !! Indices (unitless) 1862 INTEGER(i_std) :: ier !! Error handling 1949 INTEGER(i_std) :: ier !! Error handling 1950 REAL(r_std), DIMENSION(nbpt) :: Count_failure_slow !! Counter times slow reserv. does not fit irrigation demand 1951 REAL(r_std), DIMENSION(nbpt) :: Count_failure_fast !! Counter times fast reserv. does not fit irrigation demand 1952 REAL(r_std), DIMENSION(nbpt) :: Count_failure_stre !! Counter times stream reserv. does not fit irrigation demand 1953 LOGICAL :: IsFail_slow !! Logical to ask if slow reserv. does not fit irrigation demand 1954 LOGICAL :: IsFail_fast !! Logical to ask if fast reserv. does not fit irrigation demand 1955 LOGICAL :: IsFail_stre !! Logical to ask if stream reserv. does not fit irrigation demand 1863 1956 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: fast_flow_g !! Outflow from the fast reservoir (kg/dt) 1864 1957 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: slow_flow_g !! Outflow from the slow reservoir (kg/dt) … … 1872 1965 REAL(r_std), DIMENSION(nbpt, nbasmax) :: netflow_slow !! Input - Output flow to slow reservoir 1873 1966 1874 1967 REAL(r_std) :: pcent_vol_irrig !! Percentage of surface water for irrigation from total available water 1968 REAL(r_std) :: slow_wdr_dummy,fast_wdr_dummy,stre_wdr_dummy !! Dummy variables, real abstraction in each reservoir for irrigation 1969 REAL(r_std) :: pot_slow_wdr_dummy,pot_fast_wdr_dummy,pot_stre_wdr_dummy !! Dummy variables, potential abstraction in each reservoir for irrigatio 1875 1970 !! PARAMETERS 1876 1971 LOGICAL, PARAMETER :: check_reservoir = .FALSE. !! Logical to choose if we write informations when a negative amount of water is occurring in a reservoir (true/false) … … 1886 1981 totarea(:) = zero 1887 1982 totflood(:) = zero 1983 irrig_gw_source(:) = zero ! 1984 irrig_fast_source(:) = zero ! 1985 irrig_str_source(:) = zero ! 1986 Count_failure_slow(:) = zero ! 1987 Count_failure_fast(:) = zero ! 1988 Count_failure_stre(:) = zero ! 1888 1989 ! 1889 1990 ! Compute all the fluxes … … 1897 1998 ENDDO 1898 1999 ! 1899 !> The outflow fluxes from the three reservoirs are computed. 2000 !> The outflow fluxes from the three reservoirs are computed. 1900 2001 !> The outflow of volume of water Vi into the reservoir i is assumed to be linearly related to its volume. 1901 2002 !> The water travel simulated by the routing scheme is dependent on the water retention index topo_resid … … 1912 2013 IF ( route_tobasin(ig,ib) .GT. 0 ) THEN 1913 2014 ! 1914 ! Each of the fluxes is limited by the water in the reservoir and a small margin 2015 ! Each of the fluxes is limited by the water in the reservoir and a small margin 1915 2016 ! (min_reservoir) to avoid rounding errors. 1916 2017 ! … … 1923 2024 slow_flow(ig,ib) = MAX(flow, zero) 1924 2025 1925 flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst* & 2026 flow = MIN(stream_reservoir(ig,ib)/((topo_resid(ig,ib)/1000.)*stream_tcst* & 1926 2027 & MAX(un-SQRT(flood_frac_bas(ig,ib)),min_sechiba)*one_day/dt_routing),& 1927 2028 & stream_reservoir(ig,ib)-min_sechiba) … … 1944 2045 flow = MIN(floodout(ig)*totarea(ig)*pond_frac(ig)/flood_frac(ig), pond_reservoir(ig)+totflood(ig)) 1945 2046 pondex = MAX(flow - pond_reservoir(ig), zero) 1946 pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex) 2047 pond_reservoir(ig) = pond_reservoir(ig) - (flow - pondex) 1947 2048 ! 1948 2049 ! If demand was over reservoir size, we will take it out from floodplains … … 1981 2082 !- 1982 2083 !- Computing the drainage and outflow from floodplains 1983 !> Drainage from floodplains is depending on a averaged conductivity (k_litt) 2084 !> Drainage from floodplains is depending on a averaged conductivity (k_litt) 1984 2085 !> for saturated infiltration in the 'litter' layer. Flood_drainage will be 1985 2086 !> a component of the total reinfiltration that leaves the routing scheme. … … 1997 2098 DO ib=1,nbasmax 1998 2099 DO ig=1,nbpt 1999 flood_drainage(ig,ib) = zero 2100 flood_drainage(ig,ib) = zero 2000 2101 ENDDO 2001 2102 ENDDO … … 2045 2146 pond_drainage(ig,ib) = MIN(pond_reservoir(ig)*routing_area(ig,ib)/totarea(ig), & 2046 2147 & pond_frac(ig)*routing_area(ig,ib)*k_litt(ig)*dt_routing/one_day) 2047 fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib) 2148 fast_flow(ig,ib) = fast_flow(ig,ib) - pond_inflow(ig,ib) 2048 2149 ENDDO 2049 2150 ENDDO … … 2073 2174 ENDIF 2074 2175 IF (ier /= 0) CALL ipslerr_p(3,'routing_flow','Pb in allocate for fast_flow_g','','') 2075 2176 2076 2177 CALL gather(fast_flow,fast_flow_g) 2077 2178 CALL gather(slow_flow,slow_flow_g) … … 2092 2193 2093 2194 DEALLOCATE( fast_flow_g, slow_flow_g, stream_flow_g ) 2094 2195 2095 2196 CALL scatter(transport_glo,transport) 2096 2197 … … 2112 2213 DO ib=1,nbasmax 2113 2214 DO ig=1,nbpt 2114 potflood(ig,ib) = transport(ig,ib) 2215 potflood(ig,ib) = transport(ig,ib) 2115 2216 ! 2116 2217 IF ( tobeflooded(ig) > 0. .AND. potflood(ig,ib) > 0. .AND. floodtemp(ig) > tp_00 ) THEN … … 2123 2224 return_swamp(ig,ib) = swamp_cst * potflood(ig,ib) * floodindex 2124 2225 ! 2125 tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib) 2226 tobeflooded(ig) = tobeflooded(ig) - routing_area(ig,ib) 2126 2227 ! 2127 2228 ENDIF … … 2136 2237 IF (floodplains(ig) .GT. min_sechiba .AND. floodtemp(ig) .GT. tp_00) THEN 2137 2238 DO ib=1,nbasmax 2138 floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib) 2239 floods(ig,ib) = transport(ig,ib) - return_swamp(ig,ib) 2139 2240 ENDDO 2140 2241 ENDIF … … 2162 2263 ! 2163 2264 flood_reservoir(ig,ib) = flood_reservoir(ig,ib) + floods(ig,ib) - & 2164 & flood_flow(ig,ib) 2265 & flood_flow(ig,ib) 2165 2266 ! 2166 2267 pond_reservoir(ig) = pond_reservoir(ig) + pond_inflow(ig,ib) - pond_drainage(ig,ib) … … 2170 2271 WRITE(numout,*) "WARNING : negative flood reservoir at :", ig, ib, ". Problem is being corrected." 2171 2272 WRITE(numout,*) "flood_reservoir, floods, flood_flow : ", flood_reservoir(ig,ib), floods(ig,ib), & 2172 & flood_flow(ig,ib) 2273 & flood_flow(ig,ib) 2173 2274 ENDIF 2174 2275 stream_reservoir(ig,ib) = stream_reservoir(ig,ib) + flood_reservoir(ig,ib) … … 2229 2330 flood_frac_pot = (totflood(ig) / (totarea(ig)*floodcri/(beta+un)))**(beta/(beta+un)) 2230 2331 flood_frac(ig) = MIN(floodplains(ig) / totarea(ig), flood_frac_pot) 2231 ! Then we diagnose the fraction for each basin with the size of its flood_reservoir 2332 ! Then we diagnose the fraction for each basin with the size of its flood_reservoir 2232 2333 ! (flood_frac_bas may be > 1) 2233 2334 DO ib=1,nbasmax … … 2238 2339 ENDDO 2239 2340 ! We diagnose the maximum height of floodplain 2240 flood_height(ig) = (beta/(beta+1))*floodcri*(flood_frac(ig))**(un/beta) + totflood(ig)/(totarea(ig)*flood_frac(ig)) 2341 flood_height(ig) = (beta/(beta+1))*floodcri*(flood_frac(ig))**(un/beta) + totflood(ig)/(totarea(ig)*flood_frac(ig)) 2241 2342 ! And finally add the pond surface 2242 pond_frac(ig) = MIN(un-flood_frac(ig), ((betap+1)*pond_reservoir(ig) / (pondcri*totarea(ig)))**(betap/(betap+1)) ) 2343 pond_frac(ig) = MIN(un-flood_frac(ig), ((betap+1)*pond_reservoir(ig) / (pondcri*totarea(ig)))**(betap/(betap+1)) ) 2243 2344 flood_frac(ig) = flood_frac(ig) + pond_frac(ig) 2244 2345 ! … … 2265 2366 DO ig=1,nbpt 2266 2367 returnflow(ig) = returnflow(ig) + return_swamp(ig,ib) 2267 reinfiltration(ig) = reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib) 2368 reinfiltration(ig) = reinfiltration(ig) + pond_drainage(ig,ib) + flood_drainage(ig,ib) 2268 2369 ENDDO 2269 2370 ENDDO … … 2304 2405 ! 2305 2406 IF ( do_irrigation ) THEN 2306 DO ig=1,nbpt 2307 ! 2308 IF ((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. & 2309 & (runoff(ig) .LT. min_sechiba) ) THEN 2310 2311 irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - & 2407 DO ig=1,nbpt 2408 !It enters to the new irrigation module only if there is an irrigated fraction, if not irrig_netereq = zero for that cell 2409 IF ((irrig_frac(ig) .GT. min_sechiba) .AND. .NOT. old_irrig_scheme ) THEN 2410 2411 irrig_netereq(ig) = irrig_netereq(ig) + MIN( irrig_dosmax/3600*dt_routing, & 2412 root_deficit(ig) ) * soiltile(ig, irrig_st) * vegtot(ig) 2413 ! By definition, irrig_dosmax is in kg/m2 of soil tile/hour,dividing by 3600(seconds/hour) * DT_ROUTING ! 2414 ! = kg/m2 of soil tile/(routing timestep) 2415 ! irrig_netereq(kg/m2 of grid cell / routing timstep ) is equal to 2416 ! root_deficit (kg/m2 of soil tile) * soiltile*vegtot (fraction of soil tile at cell level) = kg/m2 of grid cell 2417 2418 IF (.NOT. irrigated_soiltile .AND. ( soiltile(ig,irrig_st) .GT. min_sechiba ) .AND. (vegtot(ig) .GT. min_sechiba) ) THEN 2419 ! Irrigated_soiltile asks if there is an independent soil tile for irrigated crops. If not, 2420 ! actual volume calculated for irrig_netereq assumed that the whole SOILTILE was irrigated, but in this case 2421 ! just a fraction of the irrig_st (irrigated soil tile, by default = 3) is actually irrigated, 2422 ! and irrig_netereq needs to be reduced by irrig_frac/( soiltile * vegtot ) (note that it is max = 1 thanks to irrig_frac calculation in l. 424) 2423 ! Demand(ST3)*irrig_frac/(soiltile(3)*vegtot) = irrig_netereq_In_ST3, then 2424 !irrig_netereq_In_ST3 * (soiltile(3)*vegtot) = irrig_netereq at grid scale = Demand(ST3)*irrig_frac. 2425 irrig_netereq(ig) = irrig_netereq(ig) * irrig_frac(ig) / ( soiltile(ig,irrig_st) * vegtot(ig) ) 2426 !irrig_netereq = kg/m2 of grid cell 2427 2428 ENDIF 2429 !Old irrigation scheme as in tag 2.0 2430 ELSE IF((vegtot(ig) .GT. min_sechiba) .AND. (humrel(ig) .LT. un-min_sechiba) .AND. & 2431 & (runoff(ig) .LT. min_sechiba) .AND. old_irrig_scheme) THEN 2432 2433 irrig_netereq(ig) = (irrigated(ig) / totarea(ig) ) * MAX(zero, transpot_mean(ig) - & 2312 2434 & (precip(ig)+reinfiltration(ig)) ) 2313 2314 ENDIF 2315 ! 2316 DO ib=1,nbasmax 2317 IF ( routing_area(ig,ib) .GT. 0 ) THEN 2318 2319 irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 2320 2321 irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 2322 & stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) ) 2323 2324 slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + & 2325 & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) 2326 2327 fast_reservoir(ig,ib) = MAX( zero, & 2328 & fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))) 2329 2330 stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) 2331 2332 irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 2333 2334 ENDIF 2335 ENDDO 2435 2436 ENDIF 2437 2438 DO ib=1,nbasmax 2439 IF ( routing_area(ig,ib) .GT. 0 ) THEN 2440 2441 IF (.NOT. old_irrig_scheme .AND. select_source_irrig) THEN 2442 2443 ! For irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir 2444 ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir 2445 2446 ! The new priorization scheme takes into account irrig. infrastructur according to GMIA map 2447 ! It also withdraw water according to availability, it means that it wont seek for all the water in then 2448 ! stream reservoir, even if this one could respond to the demand by itself 2449 2450 pot_slow_wdr_dummy = ( 1 - fraction_aeirrig_sw(ig)) * avail_reserve(3)*slow_reservoir(ig,ib) 2451 pot_fast_wdr_dummy = fraction_aeirrig_sw(ig) * avail_reserve(2)*fast_reservoir(ig,ib) 2452 pot_stre_wdr_dummy = fraction_aeirrig_sw(ig) * avail_reserve(1) * stream_reservoir(ig,ib) 2453 pcent_vol_irrig = zero 2454 IsFail_slow = .FALSE. ! 2455 IsFail_fast = .FALSE. ! 2456 IsFail_stre = .FALSE. ! 2457 irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 2458 2459 irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 2460 pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy) 2461 2462 !! additional IF to calculate pcent_vol_irrig, in the case the total avail. 2463 !! water is zero, I.E. when there is no water in surface and fraction_ae = 0, 2464 !! so GW is not taken into account 2465 !! Note on pcent_vol_irrig: It correspond to the fraction of available water in surface, 2466 !! considering environmental needs and irrigation equipement by source from map. It controls 2467 !! how the source of water withdrawl, especially when requirements < available water 2468 2469 IF ( (pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy) .GT. min_sechiba ) THEN 2470 2471 pcent_vol_irrig = ( pot_stre_wdr_dummy + pot_fast_wdr_dummy ) / & 2472 ( pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy ) 2473 2474 !Irrig_actual set to zero, because there is no available water. 2475 !Put to avoid negative values due to problems in the Min function 2476 irrig_actual(ig,ib) = MAX(irrig_actual(ig,ib), zero) 2477 !Already zero because pcent_vol_irrig initialized to zero 2478 !Put here to readability but not necessary 2479 !ELSE 2480 ! pcent_vol_irrig = zero 2481 2482 ENDIF 2483 2484 !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the 2485 !Volume used for irrigation that comes from GW 2486 ! Idem for irrig_fast_source and irrig_str_source 2487 2488 slow_wdr_dummy = slow_reservoir(ig,ib) 2489 slow_reservoir(ig,ib) = MAX( (un - ( un - fraction_aeirrig_sw(ig) ) * avail_reserve(3) ) * & 2490 slow_reservoir(ig,ib), slow_reservoir(ig,ib) + & 2491 MIN( - irrig_actual(ig,ib) * (un - pcent_vol_irrig ), & 2492 avail_reserve(2) * fraction_aeirrig_sw(ig) * fast_reservoir(ig,ib) + & 2493 MIN(zero, avail_reserve(1) * fraction_aeirrig_sw(ig) * stream_reservoir(ig,ib) - & 2494 pcent_vol_irrig * irrig_actual(ig,ib) ) ) ) 2495 2496 slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib) 2497 irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy 2498 2499 fast_wdr_dummy = fast_reservoir(ig,ib) 2500 fast_reservoir(ig,ib) = MAX( (un - avail_reserve(2) * fraction_aeirrig_sw(ig) ) * fast_reservoir(ig,ib) , & 2501 fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1) * fraction_aeirrig_sw(ig) * stream_reservoir(ig,ib) - & 2502 pcent_vol_irrig * irrig_actual(ig,ib) ) ) 2503 fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib) 2504 irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy 2505 2506 stre_wdr_dummy = stream_reservoir(ig,ib) 2507 stream_reservoir(ig,ib) = MAX((un - avail_reserve(1)* fraction_aeirrig_sw(ig) )*stream_reservoir(ig,ib), & 2508 stream_reservoir(ig,ib) - & 2509 pcent_vol_irrig * irrig_actual(ig,ib) ) 2510 stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib) 2511 irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy 2512 2513 irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 2514 2515 !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal 2516 !We assume that the pot. requirement = Needs * fraction of area equipped for SW/GW 2517 !In the case of surface. we also sustract the withdrawal from Fast/Stream, because both are 2518 ! considered as surface water 2519 IsFail_slow = ( ( irrig_needs(ig,ib)*(un - fraction_aeirrig_sw(ig)) ) > pot_slow_wdr_dummy ) 2520 IsFail_fast = ( ( irrig_needs(ig,ib)*fraction_aeirrig_sw(ig) - stre_wdr_dummy ) > pot_fast_wdr_dummy ) 2521 IsFail_stre = ( ( irrig_needs(ig,ib)*fraction_aeirrig_sw(ig) - fast_wdr_dummy ) > pot_stre_wdr_dummy ) 2522 2523 IF( IsFail_stre ) THEN 2524 Count_failure_stre(ig) = un 2525 ENDIF 2526 IF( IsFail_fast ) THEN 2527 Count_failure_fast(ig) = un 2528 ENDIF 2529 IF( IsFail_slow ) THEN 2530 Count_failure_slow(ig) = un 2531 ENDIF 2532 2533 ELSE IF (.NOT. old_irrig_scheme .AND. .NOT. select_source_irrig) THEN 2534 ! For irrig. scheme, available_reserve gives the amount of water available for irrigation in every reservoir 2535 ! --> avail_reserve is a vector of dimension=3, BY DEFINITION i=1 for streamflow, i=2 fast, and i=3 slow reservoir 2536 irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 2537 2538 pot_slow_wdr_dummy = avail_reserve(3)*slow_reservoir(ig,ib) 2539 pot_fast_wdr_dummy = avail_reserve(2)*fast_reservoir(ig,ib) 2540 pot_stre_wdr_dummy = avail_reserve(1)*stream_reservoir(ig,ib) 2541 IsFail_slow = .FALSE. ! 2542 IsFail_fast = .FALSE. ! 2543 IsFail_stre = .FALSE. ! 2544 2545 irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 2546 & pot_stre_wdr_dummy + pot_fast_wdr_dummy + pot_slow_wdr_dummy ) 2547 2548 !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the 2549 !Volume used for irrigation that comes from GW 2550 ! Idem for irrig_fast_source and irrig_str_source 2551 slow_wdr_dummy = slow_reservoir(ig,ib) 2552 slow_reservoir(ig,ib) = MAX( (1-avail_reserve(3) )*slow_reservoir(ig,ib), slow_reservoir(ig,ib) + & 2553 & MIN(zero, avail_reserve(2)*fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1)*stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) 2554 slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib) 2555 irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy 2556 2557 fast_wdr_dummy = fast_reservoir(ig,ib) 2558 fast_reservoir(ig,ib) = MAX( (1-avail_reserve(2) )*fast_reservoir(ig,ib) , & 2559 fast_reservoir(ig,ib) + MIN(zero, avail_reserve(1)*stream_reservoir(ig,ib)-irrig_actual(ig,ib))) 2560 fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib) 2561 irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy 2562 2563 stre_wdr_dummy = stream_reservoir(ig,ib) 2564 stream_reservoir(ig,ib) = MAX( (1-avail_reserve(1) )*stream_reservoir(ig,ib), stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) 2565 stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib) 2566 irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy 2567 2568 irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 2569 !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal 2570 ! Because it follows the old scheme, we do not separate between surface/gw, but consider that 2571 ! priority is given in this order: River, Fast and Slow reservoir. 2572 IsFail_slow = ( ( irrig_needs(ig,ib) - stre_wdr_dummy - fast_wdr_dummy ) > pot_slow_wdr_dummy ) 2573 IsFail_fast = ( ( irrig_needs(ig,ib) - stre_wdr_dummy ) > pot_fast_wdr_dummy ) 2574 IsFail_stre = ( irrig_needs(ig,ib) > pot_stre_wdr_dummy ) 2575 2576 IF( IsFail_stre ) THEN 2577 Count_failure_stre(ig) = un 2578 ENDIF 2579 IF( IsFail_fast ) THEN 2580 Count_failure_fast(ig) = un 2581 ENDIF 2582 IF( IsFail_slow ) THEN 2583 Count_failure_slow(ig) = un 2584 ENDIF 2585 2586 ELSE !Old irrigation scheme as in tag 2.0 2587 !Note for irrig_gw_source(ig): first we add the slow_reservoir volume. Then we substract the updated slow_reservoir. It should be the 2588 !Volume used for irrigation that comes from GW 2589 ! Idem for irrig_fast_source and irrig_str_source 2590 irrig_needs(ig,ib) = irrig_netereq(ig) * routing_area(ig,ib) 2591 2592 pot_slow_wdr_dummy = slow_reservoir(ig,ib) 2593 pot_fast_wdr_dummy = fast_reservoir(ig,ib) 2594 pot_stre_wdr_dummy = stream_reservoir(ig,ib) 2595 IsFail_slow = .FALSE. ! 2596 IsFail_fast = .FALSE. ! 2597 IsFail_stre = .FALSE. ! 2598 irrig_actual(ig,ib) = MIN(irrig_needs(ig,ib),& 2599 & stream_reservoir(ig,ib) + fast_reservoir(ig,ib) + slow_reservoir(ig,ib) ) 2600 2601 slow_wdr_dummy = slow_reservoir(ig,ib) 2602 slow_reservoir(ig,ib) = MAX(zero, slow_reservoir(ig,ib) + & 2603 & MIN(zero, fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib)))) 2604 slow_wdr_dummy = slow_wdr_dummy - slow_reservoir(ig,ib) 2605 irrig_gw_source(ig) = irrig_gw_source(ig) + slow_wdr_dummy 2606 2607 fast_wdr_dummy = fast_reservoir(ig,ib) 2608 fast_reservoir(ig,ib) = MAX( zero, & 2609 & fast_reservoir(ig,ib) + MIN(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib))) 2610 fast_wdr_dummy = fast_wdr_dummy - fast_reservoir(ig,ib) 2611 irrig_fast_source(ig) = irrig_fast_source(ig) + fast_wdr_dummy 2612 2613 stre_wdr_dummy = stream_reservoir(ig,ib) 2614 stream_reservoir(ig,ib) = MAX(zero, stream_reservoir(ig,ib)-irrig_actual(ig,ib) ) 2615 stre_wdr_dummy = stre_wdr_dummy - stream_reservoir(ig,ib) 2616 irrig_str_source(ig) = irrig_str_source(ig) + stre_wdr_dummy 2617 2618 irrig_deficit(ig,ib) = irrig_needs(ig,ib)-irrig_actual(ig,ib) 2619 !A reservoir is failing to give water for infiltration if pot. req > pot. withdrawal 2620 ! Because it follows the old scheme, we do not separate between surface/gw, but consider that 2621 ! priority is given in this order: River, Fast and Slow reservoir. 2622 IsFail_slow = ( ( irrig_needs(ig,ib) - stre_wdr_dummy - fast_wdr_dummy ) > pot_slow_wdr_dummy ) 2623 IsFail_fast = ( ( irrig_needs(ig,ib) - stre_wdr_dummy ) > pot_fast_wdr_dummy ) 2624 IsFail_stre = ( irrig_needs(ig,ib) > pot_stre_wdr_dummy ) 2625 2626 IF( IsFail_stre ) THEN 2627 Count_failure_stre(ig) = un 2628 ENDIF 2629 IF( IsFail_fast ) THEN 2630 Count_failure_fast(ig) = un 2631 ENDIF 2632 IF( IsFail_slow ) THEN 2633 Count_failure_slow(ig) = un 2634 ENDIF 2635 2636 ENDIF 2637 2638 ENDIF 2639 ENDDO 2336 2640 ! 2337 2641 ! Check if we cannot find the missing water in another basin of the same grid (stream reservoir only). … … 2342 2646 !> in another basin of the same grid. If there is water in the stream reservoir of this subbasin, we create some adduction 2343 2647 !> from that subbasin to the one where we need it for irrigation. 2344 !> 2345 DO ib=1,nbasmax2346 2347 stream_tot =SUM(stream_reservoir(ig,:))2348 2349 DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba)2350 2351 fi = MAXLOC(stream_reservoir(ig,:))2352 ib2 = fi(1)2353 2354 irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib),stream_reservoir(ig,ib2))2355 stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib)2356 irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib)2357 2358 stream_tot =SUM(stream_reservoir(ig,:))2359 2360 ENDDO2361 2362 ENDDO2363 ! 2364 2365 2366 2367 2368 !2369 !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes2370 !> to the one where we need it for irrigation.2371 2648 !> 2649 DO ib=1,nbasmax 2650 2651 stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:)) 2652 2653 DO WHILE ( irrig_deficit(ig,ib) > min_sechiba .AND. stream_tot > min_sechiba) 2654 2655 fi = MAXLOC(stream_reservoir(ig,:)) 2656 ib2 = fi(1) 2657 2658 irrig_adduct(ig,ib) = MIN(irrig_deficit(ig,ib), a_stream_adduction * stream_reservoir(ig,ib2)) 2659 stream_reservoir(ig,ib2) = stream_reservoir(ig,ib2)-irrig_adduct(ig,ib) 2660 irrig_deficit(ig,ib) = irrig_deficit(ig,ib)-irrig_adduct(ig,ib) 2661 2662 stream_tot = a_stream_adduction * SUM(stream_reservoir(ig,:)) 2663 2664 ENDDO 2665 2666 ENDDO 2667 ! 2668 ENDDO 2669 ! 2670 ! If we are at higher resolution we might need to look at neighboring grid boxes to find the streams 2671 ! which can feed irrigation 2672 ! 2673 !> At higher resolution (grid box smaller than 100x100km), we can import water from neighboring grid boxes 2674 !> to the one where we need it for irrigation. 2675 ! 2372 2676 IF (is_root_prc) THEN 2373 2677 ALLOCATE(irrig_deficit_glo(nbp_glo, nbasmax), stream_reservoir_glo(nbp_glo, nbasmax), & … … 2398 2702 ig2 = neighbours_g(ig,in) 2399 2703 IF (ig2 .GT. 0 ) THEN 2400 streams_around(in,:) = stream_reservoir_glo(ig2,:)2704 streams_around(in,:) = a_stream_adduction * stream_reservoir_glo(ig2,:) 2401 2705 igrd(in) = ig2 2402 2706 ENDIF … … 2409 2713 ib2=ff(2) 2410 2714 ! 2411 IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. stream_reservoir_glo(ig2,ib2) > zero ) THEN2412 adduction = MIN(irrig_deficit_glo(ig,ib), stream_reservoir_glo(ig2,ib2))2715 IF ( routing_area_glo(ig2,ib2) .GT. 0 .AND. a_stream_adduction * stream_reservoir_glo(ig2,ib2) > zero ) THEN 2716 adduction = MIN(irrig_deficit_glo(ig,ib), a_stream_adduction * stream_reservoir_glo(ig2,ib2)) 2413 2717 stream_reservoir_glo(ig2,ib2) = stream_reservoir_glo(ig2,ib2) - adduction 2414 2718 irrig_deficit_glo(ig,ib) = irrig_deficit_glo(ig,ib) - adduction … … 2476 2780 pond_diag(:) = zero 2477 2781 irrigation(:) = zero 2782 irrigdeficit(:) = zero ! 2783 irrigadduct(:) = zero ! 2478 2784 ! 2479 2785 ! … … 2482 2788 DO ig=1,nbpt 2483 2789 IF (hydrodiag(ig,ib) > 0 ) THEN 2484 hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + & 2485 & stream_flow(ig,ib) 2790 hydrographs(ig) = hydrographs(ig) + fast_flow(ig,ib) + slow_flow(ig,ib) + & 2791 & stream_flow(ig,ib) 2486 2792 slowflow_diag(ig) = slowflow_diag(ig) + slow_flow(ig,ib) 2487 2793 ENDIF … … 2491 2797 flood_diag(ig) = flood_diag(ig) + flood_reservoir(ig,ib) 2492 2798 irrigation (ig) = irrigation (ig) + irrig_actual(ig,ib) + irrig_adduct(ig,ib) 2799 irrigdeficit(ig) = irrigdeficit(ig) + irrig_deficit(ig,ib)! 2800 irrigadduct (ig) = irrigadduct (ig) + irrig_adduct(ig,ib)! 2493 2801 ENDDO 2494 2802 ENDDO … … 2502 2810 ! 2503 2811 irrigation(ig) = irrigation(ig)/totarea(ig) 2504 ! 2505 ! The three output types for the routing : endoheric basins,, rivers and 2812 irrigdeficit(ig) = irrigdeficit(ig)/totarea(ig)! 2813 irrigadduct(ig) = irrigadduct(ig)/totarea(ig)! 2814 irrig_gw_source(ig) = irrig_gw_source(ig)/totarea(ig)! 2815 irrig_fast_source(ig) = irrig_fast_source(ig)/totarea(ig)! 2816 irrig_str_source(ig) = irrig_str_source(ig)/totarea(ig)! 2817 2818 ! 2819 ! The three output types for the routing : endoheric basins,, rivers and 2506 2820 ! diffuse coastal flow. 2507 2821 ! … … 2513 2827 ! 2514 2828 flood_res = flood_diag + pond_diag 2515 2516 2517 !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it 2829 2830 2831 !! Remove water from lake reservoir if it exceeds the maximum limit and distribute it 2518 2832 !! uniformly over all possible the coastflow gridcells 2519 2833 2520 2834 ! Calculate lake_overflow and remove it from lake_reservoir 2521 2835 DO ig=1,nbpt … … 2541 2855 ! Transform from kg/grid-cell/dt_routing into m^3/grid-cell/s to match output unit of coastalflow 2542 2856 CALL xios_orchidee_send_field("lake_overflow_coast",lake_overflow_coast/mille/dt_routing) 2543 2857 2858 ! Counter of reservoir failure to assure irrigation 2859 CALL xios_orchidee_send_field("Count_failure_slow", Count_failure_slow) 2860 CALL xios_orchidee_send_field("Count_failure_fast", Count_failure_fast) 2861 CALL xios_orchidee_send_field("Count_failure_stre", Count_failure_stre) 2862 2544 2863 2545 2864 END SUBROUTINE routing_flow … … 2590 2909 ! 2591 2910 lake_reservoir(ig) = lake_reservoir(ig) + lakeinflow(ig) 2592 2911 2593 2912 IF ( doswamps ) THEN 2594 2913 ! Calculate a return flow that will be extracted from the lake reservoir and reinserted in the soil in hydrol … … 2635 2954 ! 2636 2955 IMPLICIT NONE 2637 2956 2638 2957 !! INPUT VARIABLES 2639 2958 INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) … … 2652 2971 2653 2972 !_ ================================================================================================================================ 2654 routing_area => routing_area_glo 2973 routing_area => routing_area_glo 2655 2974 topo_resid => topo_resid_glo 2656 2975 route_togrid => route_togrid_glo … … 2660 2979 hydrodiag=>hydrodiag_glo 2661 2980 hydroupbasin=>hydroupbasin_glo 2662 2981 2663 2982 IF (is_root_prc) CALL routing_diagnostic(nbp_glo, index_g, lalo_g, resolution_g, contfrac_g, nbrivers_g,basinmap_g) 2664 2983 2665 routing_area => routing_area_loc 2984 routing_area => routing_area_loc 2666 2985 topo_resid => topo_resid_loc 2667 2986 route_togrid => route_togrid_loc … … 2671 2990 hydrodiag=>hydrodiag_loc 2672 2991 hydroupbasin=>hydroupbasin_loc 2673 2992 2674 2993 CALL scatter(nbrivers_g,nbrivers) 2675 2994 CALL scatter(basinmap_g,basinmap) 2676 2995 CALL scatter(hydrodiag_glo,hydrodiag_loc) 2677 2996 CALL scatter(hydroupbasin_glo,hydroupbasin_loc) 2678 2997 2679 2998 CALL xios_orchidee_send_field("basinmap",basinmap) 2680 2999 CALL xios_orchidee_send_field("nbrivers",nbrivers) … … 2692 3011 ENDIF 2693 3012 ENDIF 2694 2695 3013 3014 2696 3015 END SUBROUTINE routing_diagnostic_p 2697 3016 … … 2699 3018 !! SUBROUTINE : routing_diagnostic 2700 3019 !! 2701 !>\BRIEF This non-parallelized subroutine gives a diagnostic of the basins used. This produces some information 3020 !>\BRIEF This non-parallelized subroutine gives a diagnostic of the basins used. This produces some information 2702 3021 !! on the rivers which are being diagnosed. 2703 3022 !! 2704 3023 !! DESCRIPTION (definitions, functional, design, flags) : As not all rivers can be monitored in the model, we will only 2705 3024 !! archive num_largest rivers. In this routine we will diagnose the num_largest largest rivers and print to the standard 2706 !! output the names of these basins and their area. The list of names of these largest rivers are taken from a list coded in the 2707 !! routine routing_names. As this standard output is not sufficient, we will also write it to a netCDF file with the routine 3025 !! output the names of these basins and their area. The list of names of these largest rivers are taken from a list coded in the 3026 !! routine routing_names. As this standard output is not sufficient, we will also write it to a netCDF file with the routine 2708 3027 !! routing_diagncfile. It is important to keep for diagnostic the fraction of the largest basins in each grid box and keep information 2709 3028 !! how they are linked one to the other. … … 2726 3045 INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) 2727 3046 INTEGER(i_std), INTENT(in) :: l_index(nbpt) !! Indices of the points on the map (unitless) 2728 REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) 3047 REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) 2729 3048 REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) 2730 3049 REAL(r_std), INTENT(in) :: contfrac(nbpt) !! Fraction of land in each grid box (unitless;0-1) … … 2741 3060 INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: nb_pts !! Number of points in the basin (unitless) 2742 3061 REAL(r_std), ALLOCATABLE, DIMENSION(:) :: totarea !! Total area of basin (m^2) 2743 REAL(r_std), ALLOCATABLE, DIMENSION(:) :: tmparea !! 3062 REAL(r_std), ALLOCATABLE, DIMENSION(:) :: tmparea !! 2744 3063 INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: topids !! The IDs of the first num_largest basins (unitless) 2745 3064 CHARACTER(LEN=25), ALLOCATABLE, DIMENSION(:) :: basin_names !! Names of the rivers (unitless) … … 2818 3137 ALLOCATE(streams_resid(num_largest), stat=ier) 2819 3138 IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for streams_resid','','') 2820 3139 2821 3140 ALLOCATE(lbasin_area(num_largest,nbpt), lbasin_uparea(num_largest,nbpt), lrivercode(num_largest,nbpt), stat=ier) 2822 3141 IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for lbasin_area','','') 2823 3142 2824 3143 IF ( .NOT. is_root_prc) THEN 2825 3144 WRITE(numout,*) "routing_diagnostic is not suitable for running in parallel" … … 2828 3147 CALL ipslerr_p(3,'routing_diagnostic','This routine is not suitable for running in parallel','','') 2829 3148 ENDIF 2830 2831 3149 3150 2832 3151 !Config Key = RIVER_DESC 2833 3152 !Config Desc = Writes out a description of the rivers 2834 3153 !Config If = RIVER_ROUTING 2835 3154 !Config Def = n 2836 !Config Help = This flag allows to write out a file containing the list of 3155 !Config Help = This flag allows to write out a file containing the list of 2837 3156 !Config rivers which are beeing simulated. It provides location of outflow 2838 3157 !Config drainage area, name and ID. … … 2880 3199 ENDDO 2881 3200 ENDDO 2882 3201 2883 3202 nb_small = MIN(nb_small, 349) 2884 3203 2885 3204 ALLOCATE(basin_names(nb_small), stat=ier) 2886 3205 IF (ier /= 0) CALL ipslerr_p(3,'routing_diagnostic','Pb in allocate for basins_names','','') … … 2922 3241 longest_river = MAX(longest_river, ic) 2923 3242 ! 2924 ! Now that we have an outflow check if it is one of the num_largest rivers. 3243 ! Now that we have an outflow check if it is one of the num_largest rivers. 2925 3244 ! In this case we keeps the location so we diagnose it. 2926 3245 ! … … 3123 3442 tuparea(1:ii) = lbasin_uparea(idbas,1:ii) 3124 3443 ! 3125 CALL routing_diagcode(ii, tpts, tptbas, tuparea, tslen, MAXVAL(tslen), allstreams, upstreamchange, tcode) 3444 CALL routing_diagcode(ii, tpts, tptbas, tuparea, tslen, MAXVAL(tslen), allstreams, upstreamchange, tcode) 3126 3445 ! 3127 3446 lrivercode(idbas,:) = 0 … … 3178 3497 ENDIF 3179 3498 ENDDO 3180 3499 3181 3500 IF ( name_found > 1 ) THEN 3182 3501 DO ic=num_largest,1,-1 … … 3199 3518 ENDDO 3200 3519 ENDIF 3201 3520 3202 3521 ENDDO 3203 3522 ! … … 3281 3600 WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. zero) 3282 3601 END IF 3283 3602 3284 3603 DEALLOCATE(pts) 3285 3604 DEALLOCATE(outpt) … … 3301 3620 WRITE(numout,*) 'Minimum topographic index :', MINVAL(topo_resid, MASK=topo_resid .GT. 0.) 3302 3621 END IF 3303 3622 3304 3623 END SUBROUTINE routing_diagnostic 3305 3624 ! … … 3308 3627 !! 3309 3628 !>\BRIEF This subroutine determines the code in the Pfafstetter system for all points 3310 !! within the given catchment. 3629 !! within the given catchment. 3311 3630 !! 3312 3631 !! DESCRIPTION (definitions, functional, design, flags) : None … … 3322 3641 !_ ================================================================================================================================ 3323 3642 3324 SUBROUTINE routing_diagcode(ip, tpts, tpbas, tuparea, tslen, ls, allstreams, upstreamchange, streamcode) 3643 SUBROUTINE routing_diagcode(ip, tpts, tpbas, tuparea, tslen, ls, allstreams, upstreamchange, streamcode) 3325 3644 ! 3326 3645 IMPLICIT NONE … … 3388 3707 IF ( streamcode(ic) == indsubbas(ib) ) THEN 3389 3708 it =it+1 3390 iw(it)=ic 3709 iw(it)=ic 3391 3710 ENDIF 3392 3711 ENDDO … … 3407 3726 ENDDO 3408 3727 ELSE 3409 ! 3728 ! 3410 3729 ! Else do the Pfafstetter numbering 3411 ! 3730 ! 3412 3731 ! 3413 3732 ! Where do we have the 4 largest change in upstream area on this stream. … … 3427 3746 ENDDO 3428 3747 ! 3429 ! Find all streams which are identical up to that junction and increase their code accordingly 3748 ! Find all streams which are identical up to that junction and increase their code accordingly 3430 3749 ! 3431 3750 DO i=1,it … … 3474 3793 !! 3475 3794 !>\BRIEF This subroutine creates a netCDF file containing all the informations 3476 !! on the largest rivers which can be used for a refined analysis. 3795 !! on the largest rivers which can be used for a refined analysis. 3477 3796 !! 3478 3797 !! DESCRIPTION (definitions, functional, design, flags) : None … … 3571 3890 & TRIM(river_file_name),'(Solution ?)') 3572 3891 ENDIF 3573 3892 3574 3893 iret = NF90_DEF_DIM(fid, 'y', jjm_g, dims(2)) 3575 3894 IF (iret /= NF90_NOERR) THEN … … 3624 3943 ! 1.3 Add attributes to the coordinate variables 3625 3944 ! 3626 iret = NF90_PUT_ATT(fid, nlonid, 'units', "degrees_east") 3945 iret = NF90_PUT_ATT(fid, nlonid, 'units', "degrees_east") 3627 3946 IF (iret /= NF90_NOERR) THEN 3628 3947 CALL ipslerr_p (3,'routing_diagncfile', 'Could not add attribut to variable '//lon_name//' for the file :', & … … 3754 4073 ELSE 3755 4074 basinfrac(i,j) = MIN(basinfrac(i,j), un) 3756 ENDIF 4075 ENDIF 3757 4076 ENDDO 3758 4077 ENDDO 3759 4078 ! 3760 ! 4079 ! 3761 4080 ! 2.2 Define the variables in the netCDF file 3762 4081 ! … … 3834 4153 ENDIF 3835 4154 ! 3836 ! Longitude of outflow point 4155 ! Longitude of outflow point 3837 4156 att_str='Longitude_of_outflow_point' 3838 4157 iret = NF90_PUT_ATT(fid, varid, att_str, lalo(outpt(ib,1),2)) … … 3877 4196 ENDIF 3878 4197 ! 3879 ! Average number of hops to go to the ocean for any stream 4198 ! Average number of hops to go to the ocean for any stream 3880 4199 att_str='Average_number_of_hops_to_ocean_for_any_stream' 3881 4200 iret = NF90_PUT_ATT(fid, varid, att_str, streams_avehops(ib)) … … 3885 4204 ENDIF 3886 4205 ! 3887 ! Maximum number of hops to go to the ocean for any stream 4206 ! Maximum number of hops to go to the ocean for any stream 3888 4207 att_str='Maximum_number_of_hops_to_ocean_for_any_stream' 3889 4208 iret = NF90_PUT_ATT(fid, varid, att_str, streams_maxhops(ib)) … … 4027 4346 !! SUBROUTINE : routing_basins_p 4028 4347 !! 4029 !>\BRIEF This parallelized subroutine computes the routing map if needed. 4348 !>\BRIEF This parallelized subroutine computes the routing map if needed. 4030 4349 !! 4031 4350 !! DESCRIPTION (definitions, functional, design, flags) : None … … 4056 4375 ! INTEGER(i_std) :: neighbours_tmp(nbpt,8) 4057 4376 ! INTEGER(i_std) :: i,j 4058 4377 4059 4378 ! DO i=1,nbp_loc 4060 4379 ! DO j=1,NbNeighb … … 4063 4382 ! ELSE 4064 4383 ! neighbours_tmp(i,j)=neighbours(i,j)+nbp_para_begin(mpi_rank)-1 4065 ! ENDIF 4384 ! ENDIF 4066 4385 ! ENDDO 4067 4386 ! ENDDO 4068 4387 4069 routing_area => routing_area_glo 4388 routing_area => routing_area_glo 4070 4389 topo_resid => topo_resid_glo 4071 4390 route_togrid => route_togrid_glo … … 4073 4392 route_nbintobas => route_nbintobas_glo 4074 4393 global_basinid => global_basinid_glo 4075 4394 4076 4395 IF (is_root_prc) CALL routing_basins(nbp_glo,lalo_g, neighbours_g, resolution_g, contfrac_g) 4077 4396 4078 routing_area => routing_area_loc 4397 routing_area => routing_area_loc 4079 4398 topo_resid => topo_resid_loc 4080 4399 route_togrid => route_togrid_loc … … 4089 4408 CALL scatter(route_nbintobas_glo,route_nbintobas_loc) 4090 4409 CALL scatter(global_basinid_glo,global_basinid_loc) 4091 4410 4092 4411 END SUBROUTINE routing_basins_p 4093 ! 4094 4412 ! 4413 4095 4414 !! ================================================================================================================================ 4096 4415 !! SUBROUTINE : routing_basins 4097 4416 !! 4098 4417 !>\BRIEF This non-parallelized subroutine reads in the map of basins and flow direction to construct 4099 !! the catchments of each grid box. 4418 !! the catchments of each grid box. 4100 4419 !! 4101 4420 !! DESCRIPTION (definitions, functional, design, flags) : 4102 !! The work is done in a number of steps which are performed locally on the 4421 !! The work is done in a number of steps which are performed locally on the 4103 4422 !! GCM grid: 4104 4423 !! 1) First we find the grid-points of the high resolution routing grid which are … … 4130 4449 INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless) 4131 4450 REAL(r_std), INTENT(in) :: lalo(nbpt,2) !! Vector of latitude and longitudes (beware of the order !) 4132 INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 4451 INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 4133 4452 !! (1=North and then cloxkwise) 4134 4453 REAL(r_std), INTENT(in) :: resolution(nbpt,2) !! The size of each grid box in X and Y (m) … … 4212 4531 ! 4213 4532 IF (debug) THEN 4214 IF (ANY(diagbox .GT. nbpt)) THEN 4533 IF (ANY(diagbox .GT. nbpt)) THEN 4215 4534 WRITE(numout,*) "Debug diganostics : nbpt, diagbox", nbpt, diagbox 4216 4535 call ipslerr_p(3,'routing_basin', & 4217 & 'Problem with diagbox in debug mode.', & 4536 & 'Problem with diagbox in debug mode.', & 4218 4537 & 'diagbox values can''t be greater than land points number.', & 4219 4538 & '(decrease diagbox wrong value)') … … 4230 4549 !Config Def = routing.nc 4231 4550 !Config Help = The file provided here should alow the routing module to 4232 !Config read the high resolution grid of basins and the flow direction 4551 !Config read the high resolution grid of basins and the flow direction 4233 4552 !Config from one mesh to the other. 4234 4553 !Config Units = [FILE] … … 4280 4599 !! Basins : Provides a uniqe ID for each basin. These IDs are also used to get 4281 4600 !! the name of the basin from the table in routine routing_names. 4282 !! 4601 !! 4283 4602 !! Topoind : is the topographic index for the retention time of the water in the 4284 !! grid box. It has been computed with the following formula : 1000 x sqrt(d^3/Dz) 4603 !! grid box. It has been computed with the following formula : 1000 x sqrt(d^3/Dz) 4285 4604 !! where d is the distance of the river from the current grid box to the next one 4286 4605 !! as indicated by the variable trip. … … 4333 4652 ! Resolution in longitude 4334 4653 ! 4335 coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos ) 4654 coslat = MAX( COS( lat_rel(ip,jp) * pi/180. ), mincos ) 4336 4655 IF ( ip .EQ. 1 ) THEN 4337 4656 resol_lu(ip,jp,1) = ABS( lon_rel(ip+1,jp) - lon_rel(ip,jp) ) * pi/180. * R_Earth * coslat … … 4361 4680 callsign = "routing_basins" 4362 4681 ok_interpol = .FALSE. 4363 ! 4682 ! 4364 4683 nix=INT(MAXVAL(resolution_g(:,1))/MAXVAL(resol_lu(:,:,1)))+2 4365 4684 njx=INT(MAXVAL(resolution_g(:,2))/MAXVAL(resol_lu(:,:,2)))+2 … … 4454 4773 ALLOCATE (nbcoastal(nbpt), coastal_basin(nbpt,nwbas), stat=ALLOC_ERR) 4455 4774 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_basins','Pb in allocate for nbcoastal','','') 4456 4775 4457 4776 ! Order all sub points in each grid_box and find the sub basins 4458 4777 ! … … 4527 4846 & nb_basin, basin_inbxid, basin_sz, basin_pts, basin_bxout, coast_pts, nwbas, basin_count,& 4528 4847 & basin_area, basin_hierarchy, basin_topoind, basin_id, basin_flowdir, outflow_grid,& 4529 & nbcoastal, coastal_basin) 4848 & nbcoastal, coastal_basin) 4530 4849 ! 4531 4850 ! … … 4548 4867 & basin_hierarchy, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, & 4549 4868 & nbcoastal, coastal_basin, invented_basins) 4550 ! 4869 ! 4551 4870 ! 4552 4871 IF (printlev>=1) WRITE(numout,*) 'The maximum number of basins in any grid :', MAXVAL(basin_count) … … 4780 5099 ! We simplify the outflow. We only need the direction normal to the 4781 5100 ! box boundary and the 4 corners. 4782 ! 5101 ! 4783 5102 ! Northern border 4784 5103 IF ( trip_bx(1,1) .EQ. 102 ) trip_bx(1,1) = 101 … … 4804 5123 DO jp=2,nbj-1 4805 5124 IF ( trip_bx(1,jp) .EQ. 106 .OR. trip_bx(1,jp) .EQ. 108 ) trip_bx(1,jp) = 107 4806 ENDDO 5125 ENDDO 4807 5126 ! 4808 5127 ! … … 4854 5173 IF ( coords(ipos+1) /= undef_sechiba) THEN 4855 5174 IF ( COUNT(coords(ipos:nb_out) == coords(ipos)) > 1 ) THEN 4856 coords(ipos:nb_out-1) = coords(ipos+1:nb_out) 5175 coords(ipos:nb_out-1) = coords(ipos+1:nb_out) 4857 5176 coords(nb_out:nb_in) = undef_sechiba 4858 5177 nb_out = nb_out - 1 … … 4890 5209 DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1) 4891 5210 ll = MINLOC(coords_tmp(:), coords_tmp /= undef_sechiba) 4892 ind(ipos) = ll(1) 5211 ind(ipos) = ll(1) 4893 5212 coords_tmp(ll(1)) = undef_sechiba 4894 5213 ipos = ipos + 1 … … 4897 5216 DO WHILE (COUNT(ABS(coords_tmp(:)-undef_sechiba) > EPSILON(undef_sechiba)*10.) >= 1) 4898 5217 ll = MAXLOC(coords_tmp(:), coords_tmp /= undef_sechiba) 4899 ind(ipos) = ll(1) 5218 ind(ipos) = ll(1) 4900 5219 coords_tmp(ll(1)) = undef_sechiba 4901 5220 ipos = ipos + 1 … … 5066 5385 coast_pts(sz(ipb)) = bname(trans(ip)) 5067 5386 sz(trans(ip)) = 0 5068 pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1) 5069 pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2) 5387 pts(ipb, sz(ipb), 1) = pts(trans(ip), 1, 1) 5388 pts(ipb, sz(ipb), 2) = pts(trans(ip), 1, 2) 5070 5389 ENDDO 5071 5390 ENDIF … … 5116 5435 ENDIF 5117 5436 ! 5118 ! Third issue : we have more than one outflow from the boxes. This could be 5437 ! Third issue : we have more than one outflow from the boxes. This could be 5119 5438 ! - flow into 2 or more neighboring GCM grids 5120 5439 ! - flow into a neighboring GCM grids and into the ocean or be a return flow (=97. =98, =99) … … 5181 5500 ! 5182 5501 ! 5183 ! Sort the output by size of the various basins. 5502 ! Sort the output by size of the various basins. 5184 5503 ! 5185 5504 nb_basin = COUNT(sz(1:nbb) .GT. 0) … … 5218 5537 CALL ipslerr_p(3,'routing_findbasins','Probleme less outflow points than basins','','') 5219 5538 ENDIF 5220 5539 5221 5540 END SUBROUTINE routing_findbasins 5222 5541 ! … … 5224 5543 !! SUBROUTINE : routing_simplify 5225 5544 !! 5226 !>\BRIEF This subroutine symplifies the routing out of each basin by taking 5545 !>\BRIEF This subroutine symplifies the routing out of each basin by taking 5227 5546 !! out redundancies at the borders of the GCM box. 5228 5547 !! The aim is to have only one outflow point per basin and grid box. 5229 !! But here we will not change the direction of the outflow. 5548 !! But here we will not change the direction of the outflow. 5230 5549 !! 5231 5550 !! DESCRIPTION (definitions, functional, design, flags) : None … … 5353 5672 IF ( trip_flow(ip,jp,1) .EQ. outflow(ismall,1) .AND.& 5354 5673 & trip_flow(ip,jp,2) .EQ. outflow(ismall,2) ) THEN 5355 ! Now that we have found a point of the smallest sub-basin we 5674 ! Now that we have found a point of the smallest sub-basin we 5356 5675 ! look around for another sub-basin 5357 5676 ib = 1 5358 5677 not_found = .TRUE. 5359 DO WHILE ( not_found .AND. ib .LE. itodo ) 5678 DO WHILE ( not_found .AND. ib .LE. itodo ) 5360 5679 IF ( ib .NE. ill(1) ) THEN 5361 5680 ibas = todopt(ib) … … 5415 5734 !! 5416 5735 !>\BRIEF This subroutine cuts the original basin which has more than one outflow 5417 !! into as many subbasins as outflow directions. 5736 !! into as many subbasins as outflow directions. 5418 5737 !! 5419 5738 !! DESCRIPTION (definitions, functional, design, flags) : None … … 5499 5818 ! 5500 5819 CALL routing_findrout(nbi, nbj, trip_tmp, basin_sz, basin_inbxid, nbout, outflow, trip_flow, outsz) 5501 ! 5820 ! 5502 5821 ! IF ( debug ) THEN 5503 5822 ! DO ib = nb_in+1,nb … … 5516 5835 ! Take out the small sub-basins. That is those which have only one grid box 5517 5836 ! This is only done if we need to save space in the number of basins. Else we 5518 ! can take it easy and keep diverging sub-basins for the moment. 5837 ! can take it easy and keep diverging sub-basins for the moment. 5519 5838 ! 5520 5839 IF ( nbbasins .GE. nbasmax ) THEN … … 5543 5862 DO ibb=1,nbout 5544 5863 IF ( iip .EQ. outflow(ibb,1) .AND. jjp .EQ. outflow(ibb,2) ) THEN 5545 outsz(ibb) = outsz(ibb)+1 5864 outsz(ibb) = outsz(ibb)+1 5546 5865 trip_flow(ip,jp,1) = outflow(ibb,1) 5547 5866 trip_flow(ip,jp,2) = outflow(ibb,2) … … 5584 5903 ENDIF 5585 5904 ENDDO 5586 ! A short verification 5905 ! A short verification 5587 5906 IF ( SUM(sz(nb_in+1:nb)) .NE. basin_sz) THEN 5588 5907 WRITE(numout,*) 'Lost some points while spliting the basin' … … 5598 5917 CALL ipslerr_p(3,'routing_cutbasin','Lost some points while spliting the basin','','') 5599 5918 ENDIF 5600 5919 5601 5920 IF ( debug ) THEN 5602 5921 DO ib = nb_in+1,nb … … 5623 5942 !! 5624 5943 !>\BRIEF This subroutine finds, for each point, the distance to the outflow 5625 !! point along the flowlines of the basin. 5944 !! point along the flowlines of the basin. 5626 5945 !! 5627 5946 !! DESCRIPTION (definitions, functional, design, flags) : None … … 5691 6010 ENDIF 5692 6011 ! 5693 DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml) 6012 DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. iml*jml) 5694 6013 cnt = cnt + 1 5695 6014 ntripi = ntripi + inc(trp,1) … … 5707 6026 trp = NINT(trip(ntripi, ntripj)) 5708 6027 ENDDO 5709 6028 5710 6029 IF ( cnt .EQ. iml*jml) THEN 5711 6030 WRITE(numout,*) 'We could not route point (routing_findrout) :', ip, jp … … 5717 6036 CALL ipslerr_p(3,'routing_hierarchy','We could not route point','','') 5718 6037 ENDIF 5719 6038 5720 6039 hierarchy(ip,jp) = topohier 5721 6040 5722 6041 ENDIF 5723 6042 ENDDO … … 5731 6050 !! 5732 6051 !>\BRIEF This subroutine simply computes the route to each outflow point 5733 !! and returns the outflow point for each point in the basin. 6052 !! and returns the outflow point for each point in the basin. 5734 6053 !! 5735 6054 !! DESCRIPTION (definitions, functional, design, flags) : None … … 5759 6078 INTEGER(i_std), DIMENSION(nbvmax,2), INTENT(out) :: outflow !! 5760 6079 INTEGER(i_std), DIMENSION(nbvmax,nbvmax,2), INTENT(out) :: trip_flow !! 5761 INTEGER(i_std), INTENT(out) :: nbout !! 6080 INTEGER(i_std), INTENT(out) :: nbout !! 5762 6081 INTEGER(i_std), DIMENSION(nbvmax), INTENT(out) :: outsz !! 5763 6082 ! … … 5790 6109 ! 5791 6110 ! 5792 ! Get the outflows and determine for each point to which outflow point it belong 6111 ! Get the outflows and determine for each point to which outflow point it belong 5793 6112 ! 5794 6113 nbout = 0 … … 5800 6119 outflow(nbout,1) = ip 5801 6120 outflow(nbout,2) = jp 5802 ENDIF 6121 ENDIF 5803 6122 IF ( trip(ip,jp) .GT. 0) THEN 5804 6123 trip_flow(ip,jp,1) = ip … … 5815 6134 trp = trip(trip_flow(ip,jp,1), trip_flow(ip,jp,2)) 5816 6135 cnt = 0 5817 DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj) 6136 DO WHILE ( trp .GT. 0 .AND. trp .LT. 9 .AND. cnt .LT. nbi*nbj) 5818 6137 cnt = cnt + 1 5819 6138 trip_flow(ip,jp,1) = trip_flow(ip,jp,1) + inc(trp,1) … … 5862 6181 !! 5863 6182 !>\BRIEF This subroutine puts the basins found for grid box in the global map. 5864 !! Connection can only be made later when all information is together. 6183 !! Connection can only be made later when all information is together. 5865 6184 !! 5866 6185 !! DESCRIPTION (definitions, functional, design, flags) : None … … 5888 6207 INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) 5889 6208 INTEGER(i_std), INTENT (in) :: ib !! Current basin (unitless) 5890 INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 6209 INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point 5891 6210 !! (1=North and then clockwise) 5892 !! LOCAL VARIABLES 6211 !! LOCAL VARIABLES 5893 6212 REAL(r_std), DIMENSION(nbvmax,nbvmax) :: area_bx !! Area of each small box in the grid box (m^2) 5894 6213 INTEGER(i_std), DIMENSION(nbvmax,nbvmax) :: trip_bx !! The trip field for each of the smaller boxes (unitless) … … 5924 6243 ! 5925 6244 basin_count(ib) = basin_count(ib)+1 5926 if (basin_count(ib) > nwbas) then 6245 if (basin_count(ib) > nwbas) then 5927 6246 WRITE(numout,*) 'ib=',ib 5928 6247 call ipslerr_p(3,'routing_globalize', & 5929 & 'Problem with basin_count : ', & 6248 & 'Problem with basin_count : ', & 5930 6249 & 'It is greater than number of allocated basin nwbas.', & 5931 6250 & '(stop to count basins)') 5932 endif 6251 endif 5933 6252 basin_id(ib,basin_count(ib)) = basin_inbxid(ij) 5934 6253 ! … … 6025 6344 !! SUBROUTINE : routing_linkup 6026 6345 !! 6027 !>\BRIEF This subroutine makes the connections between the basins and ensure global coherence. 6346 !>\BRIEF This subroutine makes the connections between the basins and ensure global coherence. 6028 6347 !! 6029 6348 !! DESCRIPTION (definitions, functional, design, flags) : … … 6078 6397 INTEGER(i_std) :: ff(1) !! 6079 6398 ! 6080 ! ERRORS 6399 ! ERRORS 6081 6400 LOGICAL :: error1, error2, error3, error4, error5 !! (true/false) 6082 6401 ! … … 6084 6403 LOGICAL, PARAMETER :: check = .TRUE. !! (true/false) 6085 6404 6086 !_ ================================================================================================================================ 6405 !_ ================================================================================================================================ 6087 6406 error1=.FALSE. 6088 6407 error2=.FALSE. … … 6133 6452 ! 6134 6453 ENDDO 6135 ! 6454 ! 6136 6455 ! If we have more than one basin we will take the one which is lowest 6137 6456 ! in the hierarchy. … … 6146 6465 bop = sbl 6147 6466 ELSE 6148 ! The same hierarchy is allowed if both grids flow in about 6467 ! The same hierarchy is allowed if both grids flow in about 6149 6468 ! the same direction : 6150 6469 IF ( ( MOD(basin_flowdir(inp,sbl)+1-1, 8)+1 .EQ. basin_flowdir(sp,sb)).OR. & … … 6179 6498 ! 6180 6499 ! In case the outflow point was ocean or we did not find the correct basin we start to look 6181 ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding 6500 ! around. We find two options for the outflow direction (dp1 & dm1) and the corresponding 6182 6501 ! basin index (bp1 & bm1). 6183 6502 ! … … 6202 6521 bp1 = sbl 6203 6522 ELSE 6204 ! The same hierarchy is allowed if both grids flow in about 6523 ! The same hierarchy is allowed if both grids flow in about 6205 6524 ! the same direction : 6206 6525 angle=MODULO(basin_flowdir(dp1,sbl)-basin_flowdir(sp,sb)+8,8) … … 6223 6542 IF ( basin_hierarchy(sp,sb) .GT. basin_hierarchy(dm1,sbl) ) THEN 6224 6543 bm1 = sbl 6225 ELSE 6226 ! The same hierarchy is allowed if both grids flow in about 6544 ELSE 6545 ! The same hierarchy is allowed if both grids flow in about 6227 6546 ! the same direction : 6228 6547 angle=MODULO(basin_flowdir(dm1,sbl)-basin_flowdir(sp,sb)+8,8) … … 6240 6559 ! First deal with the case on land. 6241 6560 ! 6242 ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1 6561 ! For that we need to check if the water will be able to flow out of the grid dp1 or dm1 6243 6562 ! and not return to our current grid. If it is the current grid 6244 6563 ! then we can not do anything with that neighbour. Thus we set the … … 6250 6569 IF (basin_flowdir(dp1,bp1) .GT. 0) THEN 6251 6570 outdp1 = neighbours(dp1,basin_flowdir(dp1,bp1)) 6252 IF ( outdp1 .EQ. sp ) outdp1 = undef_int 6571 IF ( outdp1 .EQ. sp ) outdp1 = undef_int 6253 6572 ELSE 6254 6573 outdp1 = nbpt + 1 … … 6270 6589 bop = undef_int 6271 6590 IF ( outdp1 .LT. undef_int .AND. outdm1 .LT. undef_int) THEN 6272 ! 6591 ! 6273 6592 ! In this case we let the current basin flow into the smaller one 6274 6593 ! … … 6292 6611 ELSE 6293 6612 ! 6294 ! Now we are at the point where none of the neighboring points is suitable 6613 ! Now we are at the point where none of the neighboring points is suitable 6295 6614 ! or we have a coastal point. 6296 6615 ! … … 6301 6620 ELSE 6302 6621 ! 6303 ! If we are on a land point with only land neighbors but no one suitable to let the 6622 ! If we are on a land point with only land neighbors but no one suitable to let the 6304 6623 ! water flow into we have to look for a solution in the current grid box. 6305 ! 6624 ! 6306 6625 ! 6307 6626 IF ( bp1 .LT. 0 .AND. bm1 .LT. 0 ) THEN … … 6402 6721 & .AND. basin_flowdir(sp,sb) .GT. 0) THEN 6403 6722 ! 6404 6723 6405 6724 IF (check) & 6406 6725 WRITE(numout,*) 'There is no reason to here, this part of the code should be dead :', sp,sb … … 6451 6770 IF (error2) THEN 6452 6771 CALL ipslerr_p(3,'routing_linkup', & 6453 & 'In the routine which make connections between the basins and ensure global coherence,', & 6772 & 'In the routine which make connections between the basins and ensure global coherence,', & 6454 6773 & 'there is a problem with outflow linkup without any valid direction. Try with check=.TRUE.', & 6455 6774 & '(Perhaps there is a problem with the grid.)') … … 6461 6780 ENDIF 6462 6781 IF (error4) THEN 6463 WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", & 6464 & " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, & 6782 WRITE(numout,*) " routing_linkup : (sbl .NE. sb) .AND. (basin_id(sp,sbl) .EQ. bid) ", & 6783 & " .AND. (basin_hierarchy(sp,sb) .GT. basin_hierarchy(sp,sbl))",sbl,sb,basin_id(sp,sbl),bid, & 6465 6784 & basin_hierarchy(sp,sb),basin_hierarchy(sp,sbl) 6466 6785 CALL ipslerr_p(3,'routing_linkup', & … … 6471 6790 WRITE(numout,*) 'We could not find the basin into which we need to flow' 6472 6791 WRITE(numout,*) 'Grid point ', sp, ' and basin ', sb 6473 WRITE(numout,*) 'Explored neighbours :', dm1, dp1 6792 WRITE(numout,*) 'Explored neighbours :', dm1, dp1 6474 6793 WRITE(numout,*) 'Outflow direction :', basin_flowdir(sp,sb) 6475 6794 WRITE(numout,*) 'Outlfow grid :', outflow_grid(sp,sb) … … 6507 6826 !! will know how much area is upstream. It will help decide how to procede with the 6508 6827 !! the truncation later and allow to set correctly in outflow_grid the distinction 6509 !! between coastal and river flow. 6828 !! between coastal and river flow. 6510 6829 !! 6511 6830 !! DESCRIPTION (definitions, functional, design, flags) : None … … 6613 6932 WRITE(numout,*) 'The largest FETCH :', MAXVAL(fetch_basin) 6614 6933 ! 6615 ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow 6616 ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow 6934 ! Now we set for the 'num_largest' largest basins the outflow condition as stream flow 6935 ! (i.e. outflow_grid = -1) and all other outflow basins are set to coastal flow 6617 6936 ! (i.e. outflow_grid = -2). The return flow is not touched. 6618 6937 ! … … 6654 6973 !>\BRIEF This subroutine reduces the number of basins per grid to the value chosen by the user. 6655 6974 !! It also computes the final field which will be used to route the water at the 6656 !! requested truncation. 6657 !! 6658 !! DESCRIPTION (definitions, functional, design, flags) : 6975 !! requested truncation. 6976 !! 6977 !! DESCRIPTION (definitions, functional, design, flags) : 6659 6978 !! Truncate if needed and find the path closest to the high resolution data. 6660 6979 !! 6661 6980 !! The algorithm : 6662 !! 6981 !! 6663 6982 !! We only go through this procedure only as many times as there are basins to take out at most. 6664 6983 !! This is important as it allows the simplifications to spread from one grid to the other. 6665 !! The for each step of the iteration and at each grid point we check the following options for 6984 !! The for each step of the iteration and at each grid point we check the following options for 6666 6985 !! simplifying the pathways of water : 6667 6986 !! 1) If the basin of a grid flows into another basin of the same grid. Kill the one which only 6668 6987 !! served as a transition 6669 !! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow. 6670 !! We kill the smallest one and put into the largest basin. There is no need to manage many 6671 !! basins going into the ocean as coastal flows. 6672 !! 3) If we have streams run in parallel from one gird box to the others (that is these are 6988 !! 2) If in one grid box we have a number of basins which flow into the ocean as coastal flow. 6989 !! We kill the smallest one and put into the largest basin. There is no need to manage many 6990 !! basins going into the ocean as coastal flows. 6991 !! 3) If we have streams run in parallel from one gird box to the others (that is these are 6673 6992 !! different basins) we will put the smaller one in the larger one. This may hapen at any 6674 6993 !! level of the flow but in theory it should propagate downstream. 6675 6994 !! 4) If we have two basins with the same ID but flow into different grid boxes we sacrifice 6676 !! the smallest one and route it through the largest. 6995 !! the smallest one and route it through the largest. 6677 6996 !! 6678 6997 !! Obviously if any of the options find something then we skip the rest and take out the basin.:\n … … 6758 7077 ENDDO 6759 7078 ! 6760 ! Go through the basins which need to be truncated. 7079 ! Go through the basins which need to be truncated. 6761 7080 ! 6762 7081 DO ibt=1,nbtruncate … … 6782 7101 ENDDO 6783 7102 ! 6784 ! 2) Merge two basins which flow into the ocean as coastal or return flow 6785 ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and 7103 ! 2) Merge two basins which flow into the ocean as coastal or return flow 7104 ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and 6786 7105 ! have not found anything yet! 6787 7106 ! … … 6818 7137 ENDIF 6819 7138 ! 6820 ! 3) If we have basins flowing into the same grid but different basins then we put them 7139 ! 3) If we have basins flowing into the same grid but different basins then we put them 6821 7140 ! together. Obviously we first work with the grid which has most streams running into it 6822 7141 ! and putting the smallest in the largests catchments. … … 6947 7266 ENDDO 6948 7267 ! 6949 ! 7268 ! 6950 7269 ENDDO 6951 7270 ! … … 7103 7422 WRITE(numout,*) 'We should not be here as the basin flows into the pampa' 7104 7423 WRITE(numout,*) 'Last correct point :', pold, bold 7105 WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold) 7106 WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold) 7424 WRITE(numout,*) 'It pointed to in the new variables :', route_togrid(pold, bold),route_tobasin(pold, bold) 7425 WRITE(numout,*) 'The old variables gave :', outflow_grid(pold, bold), outflow_basin(pold, bold) 7107 7426 WRITE(numout,*) 'Where we ended up :', igrif,ibasf 7108 7427 CALL ipslerr_p(3,'routing_truncate','Problem with routing..','','') … … 7139 7458 DO ib=1, pickmax 7140 7459 ff = MAXLOC(floflo) 7141 ! tdo - To take into account rivers that do not flow to the oceans 7460 ! tdo - To take into account rivers that do not flow to the oceans 7142 7461 IF ( route_tobasin(ff(1), ff(2)) .GT. nbasmax ) THEN 7143 7462 ! IF ( route_tobasin(ff(1), ff(2)) .EQ. nbasmax + 2) THEN … … 7169 7488 !! 7170 7489 !>\BRIEF The aim of this subroutine is to kill a basin (that is put into another larger one). 7171 !! When we do this we need to be careful and change all associated variables. 7490 !! When we do this we need to be careful and change all associated variables. 7172 7491 !! 7173 7492 !! DESCRIPTION (definitions, functional, design, flags) : None … … 7232 7551 inf = 0 7233 7552 DO WHILE (igrif .GT. 0) 7234 fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill) 7553 fetch_basin(igrif,ibasf) = fetch_basin(igrif,ibasf) + fetch_basin(ib, tokill) 7235 7554 it = outflow_grid(igrif, ibasf) 7236 7555 ibasf = outflow_basin(igrif, ibasf) … … 7249 7568 ibasf = outflow_basin(igrif, ibasf) 7250 7569 igrif = it 7251 ENDDO 7570 ENDDO 7252 7571 ! 7253 7572 ! Redirect the flows which went into the basin to be killed before we change everything … … 7291 7610 ENDIF 7292 7611 inflow_number(ing, inb) = inflow_number(ing, inb) - 1 7293 7612 7294 7613 ENDIF 7295 7614 ! … … 7340 7659 basin_count(ib) = basin_count(ib) - 1 7341 7660 ! 7342 END SUBROUTINE routing_killbas 7661 END SUBROUTINE routing_killbas 7343 7662 ! 7344 7663 !! ================================================================================================================================ … … 7741 8060 ! 7742 8061 !! ================================================================================================================================ 7743 !! SUBROUTINE : routing_ irrigmap8062 !! SUBROUTINE : routing_floodmap 7744 8063 !! 7745 8064 !>\BRIEF This subroutine interpolates the 0.5x0.5 degree based map of irrigated areas to the resolution of the model. … … 7747 8066 !! DESCRIPTION (definitions, functional, design, flags) : None 7748 8067 !! 7749 !! RECENT CHANGE(S): None8068 !! RECENT CHANGE(S): Irrigated is interpolated from slowproc as irrigated_next 7750 8069 !! 7751 8070 !! MAIN OUTPUT VARIABLE(S): … … 7757 8076 !_ ================================================================================================================================ 7758 8077 7759 SUBROUTINE routing_ irrigmap (nbpt, index, lalo, neighbours, resolution, contfrac, &7760 & init_irrig, irrigated,init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id)8078 SUBROUTINE routing_floodmap (nbpt, index, lalo, neighbours, resolution, contfrac, & 8079 & init_flood, floodplains, init_swamp, swamp, hist_id, hist2_id) 7761 8080 ! 7762 8081 IMPLICIT NONE … … 7780 8099 INTEGER(i_std), INTENT(in) :: hist_id !! Access to history file (unitless) 7781 8100 INTEGER(i_std), INTENT(in) :: hist2_id !! Access to history file 2 (unitless) 7782 LOGICAL, INTENT(in) :: init_irrig !! Logical to initialize the irrigation (true/false) 8101 7783 8102 LOGICAL, INTENT(in) :: init_flood !! Logical to initialize the floodplains (true/false) 7784 8103 LOGICAL, INTENT(in) :: init_swamp !! Logical to initialize the swamps (true/false) 7785 8104 ! 7786 8105 !! OUTPUT VARIABLES 7787 REAL(r_std), INTENT(out) :: irrigated(:) !! Irrigated surface in each grid box (m^2) 8106 7788 8107 REAL(r_std), INTENT(out) :: floodplains(:) !! Surface which can be inundated in each grid box (m^2) 7789 8108 REAL(r_std), INTENT(out) :: swamp(:) !! Surface which can be swamp in each grid box (m^2) … … 7791 8110 !! LOCAL VARIABLES 7792 8111 ! Interpolation variables 7793 ! 8112 ! 7794 8113 INTEGER(i_std) :: nbpmax, nix, njx, fopt !! 7795 8114 CHARACTER(LEN=30) :: callsign !! … … 7807 8126 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: latrel !! Latitude 7808 8127 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lonrel !! Longitude 7809 REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: irrigated_frac !! Irrigated fraction of the grid box (unitless;0-1) 8128 7810 8129 REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: flood_fracmax !! Maximal flooded fraction of the grid box (unitless;0-1) 7811 REAL(r_std) :: area_irrig !! Irrigated surface in the grid box (m^2) 8130 7812 8131 REAL(r_std) :: area_flood(ntype) !! Flooded surface in the grid box (m^2) 7813 8132 REAL(r_std) :: resolution_1 !! temporary variable … … 7826 8145 !Config with the area in m^2 of the area irrigated within each 7827 8146 !Config 0.5 0.5 deg grid box. The map currently used is the one 7828 !Config developed by the Center for Environmental Systems Research 8147 !Config developed by the Center for Environmental Systems Research 7829 8148 !Config in Kassel (1995). 7830 8149 !Config Units = [FILE] … … 7851 8170 ! 7852 8171 ALLOCATE (latrel(iml,jml), STAT=ALLOC_ERR) 7853 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_ irrigmap','Pb in allocate for latrel','','')8172 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for latrel','','') 7854 8173 7855 8174 ALLOCATE (lonrel(iml,jml), STAT=ALLOC_ERR) 7856 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for lonrel','','') 7857 7858 ALLOCATE (irrigated_frac(iml,jml), STAT=ALLOC_ERR) 7859 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_irrigmap','Pb in allocate for irrigated_frac','','') 8175 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for lonrel','','') 7860 8176 7861 8177 ALLOCATE (flood_fracmax(iml,jml,ntype), STAT=ALLOC_ERR) 7862 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_ irrigmap','Pb in allocate for flood_fracmax','','')8178 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for flood_fracmax','','') 7863 8179 7864 8180 IF (is_root_prc) CALL flinopen(filename, .FALSE., iml, jml, lml, lonrel, latrel, lev, tml, itau, date, dt, fid) … … 7867 8183 CALL bcast(latrel) 7868 8184 ! 7869 IF (is_root_prc) CALL flinget(fid, 'irrig', iml, jml, lml, tml, 1, 1, irrigated_frac)7870 CALL bcast(irrigated_frac)7871 8185 IF (is_root_prc) CALL flinget(fid, 'lake', iml, jml, lml, tml, 1, 1, flood_fracmax(:,:,ilake)) 7872 8186 IF (is_root_prc) CALL flinget(fid, 'dam', iml, jml, lml, tml, 1, 1, flood_fracmax(:,:,idam)) … … 7883 8197 DO ip=1,iml 7884 8198 DO jp=1,jml 7885 ! 7886 IF ( irrigated_frac(ip,jp) .LT. undef_sechiba-un) THEN 7887 irrigated_frac(ip,jp) = irrigated_frac(ip,jp)/100. 7888 IF ( irrigated_frac(ip,jp) < 0.005 ) irrigated_frac(ip,jp) = zero 7889 ENDIF 7890 ! 8199 7891 8200 DO itype=1,ntype 7892 8201 IF ( flood_fracmax(ip,jp,itype) .LT. undef_sechiba-1.) THEN … … 7898 8207 ENDDO 7899 8208 ENDDO 7900 8209 7901 8210 IF (printlev>=2) THEN 7902 8211 WRITE(numout,*) 'lonrel : ', MAXVAL(lonrel), MINVAL(lonrel) 7903 8212 WRITE(numout,*) 'latrel : ', MAXVAL(latrel), MINVAL(latrel) 7904 WRITE(numout,*) 'irrigated_frac : ', MINVAL(irrigated_frac, MASK=irrigated_frac .GT. 0), & 7905 MAXVAL(irrigated_frac, MASK=irrigated_frac .LT. undef_sechiba) 8213 7906 8214 WRITE(numout,*) 'flood_fracmax : ', MINVAL(flood_fracmax, MASK=flood_fracmax .GT. 0), & 7907 8215 MAXVAL(flood_fracmax, MASK=flood_fracmax .LT. undef_sechiba) … … 7911 8219 ! 7912 8220 ALLOCATE(resol_lu(iml,jml,2), STAT=ALLOC_ERR) 7913 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_ irrigmap','Pb in allocate for resol_lu','','')8221 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for resol_lu','','') 7914 8222 7915 8223 ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR) 7916 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_ irrigmap','Pb in allocate for mask','','')8224 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for mask','','') 7917 8225 mask(:,:) = 0 7918 8226 … … 7930 8238 ! Resolution in longitude 7931 8239 ! 7932 coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos ) 8240 coslat = MAX( COS( latrel(ip,jp) * pi/180. ), mincos ) 7933 8241 IF ( ip .EQ. 1 ) THEN 7934 8242 resol_lu(ip,jp,1) = ABS( lonrel(ip+1,jp) - lonrel(ip,jp) ) * pi/180. * R_Earth * coslat … … 7955 8263 ! Some lmargin is taken. 7956 8264 ! 7957 callsign = ' Irrigationmap'8265 callsign = 'Flood map' 7958 8266 ok_interpol = .FALSE. 7959 8267 IF (is_root_prc) THEN … … 7969 8277 7970 8278 ALLOCATE(irrsub_index(nbpt, nbpmax, 2), STAT=ALLOC_ERR) 7971 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_ irrigmap','Pb in allocate for irrsub_index','','')8279 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for irrsub_index','','') 7972 8280 irrsub_index(:,:,:)=0 7973 8281 7974 8282 ALLOCATE(irrsub_area(nbpt, nbpmax), STAT=ALLOC_ERR) 7975 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_ irrigmap','Pb in allocate for irrsub_area','','')8283 IF (ALLOC_ERR /= 0) CALL ipslerr_p(3,'routing_floodmap','Pb in allocate for irrsub_area','','') 7976 8284 irrsub_area(:,:)=zero 7977 8285 … … 7982 8290 ! 7983 8291 WHERE (irrsub_area < 0) irrsub_area=zero 7984 ! 8292 ! 7985 8293 ! Test here if not all sub_area are larger than 0 if so, then we need to increase nbpmax 7986 8294 ! 7987 8295 DO ib=1,nbpt 7988 8296 ! 7989 area_irrig = 0.07990 8297 area_flood = 0.0 7991 8298 ! … … 7994 8301 ip = irrsub_index(ib, fopt, 1) 7995 8302 jp = irrsub_index(ib, fopt, 2) 7996 !7997 IF (irrigated_frac(ip,jp) .LT. undef_sechiba-1.) THEN7998 area_irrig = area_irrig + irrsub_area(ib,fopt)*irrigated_frac(ip,jp)7999 ENDIF8000 8303 ! 8001 8304 DO itype=1,ntype … … 8006 8309 ENDDO 8007 8310 ! 8008 ! Put the total irrigated and flooded areas in the output variables 8009 ! 8010 IF ( init_irrig ) THEN 8311 ! Put the total flooded areas in the output variables 8312 ! 8313 ! 8314 IF ( init_flood ) THEN 8011 8315 ! if we are at the poles resolution(ib,1) = 0 8012 8316 IF (resolution(ib,1) == 0) THEN 8013 ! use pi*resolution(ib,2) to get the disc area8014 resolution_1 = pi*resolution(ib,2)8015 ELSE8016 resolution_1 = resolution(ib,1)8017 END IF8018 irrigated(ib) = MIN(area_irrig, resolution_1*resolution(ib,2)*contfrac(ib))8019 IF ( irrigated(ib) < 0 ) THEN8020 WRITE(numout,*) 'We have a problem here : ', irrigated(ib)8021 WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2)8022 WRITE(numout,*) area_irrig8023 CALL ipslerr_p(3,'routing_irrigmap','Problem with irrigated...','','')8024 ENDIF8025 !!$ ! Compute a diagnostic of the map.8026 !!$ IF(contfrac(ib).GT.zero) THEN8027 !!$ irrigmap (ib) = irrigated(ib) / ( resolution(ib,1)*resolution(ib,2)*contfrac(ib) )8028 !!$ ELSE8029 !!$ irrigmap (ib) = zero8030 !!$ ENDIF8031 !8032 ENDIF8033 !8034 IF ( init_flood ) THEN8035 ! if we are at the poles resolution(ib,1) = 08036 IF (resolution(ib,1) == 0) THEN8037 8317 ! use pi*resolution(ib,2) to get the disc area 8038 8318 resolution_1 = pi*resolution(ib,2) … … 8043 8323 & resolution_1*resolution(ib,2)*contfrac(ib)) 8044 8324 IF ( floodplains(ib) < 0 ) THEN 8045 WRITE(numout,*) 'We have a problem here : ', floodplains(ib) 8325 WRITE(numout,*) 'We have a problem here : ', floodplains(ib) 8046 8326 WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 8047 8327 WRITE(numout,*) area_flood 8048 CALL ipslerr_p(3,'routing_ irrigmap','Problem with floodplains..','','')8328 CALL ipslerr_p(3,'routing_floodmap','Problem with floodplains..','','') 8049 8329 ENDIF 8050 8330 !!$ ! Compute a diagnostic of the map. … … 8058 8338 IF ( init_swamp ) THEN 8059 8339 ! if we are at the poles resolution(ib,1) = 0 8060 IF (resolution(ib,1) == 0) THEN 8061 ! use pi*resolution(ib,2) to get the disc area 8340 IF (resolution(ib,1) == 0) THEN 8341 ! use pi*resolution(ib,2) to get the disc area 8062 8342 resolution_1 = pi*resolution(ib,2) 8063 8343 ELSE … … 8066 8346 swamp(ib) = MIN(area_flood(iswamp), resolution_1*resolution(ib,2)*contfrac(ib)) 8067 8347 IF ( swamp(ib) < 0 ) THEN 8068 WRITE(numout,*) 'We have a problem here : ', swamp(ib) 8348 WRITE(numout,*) 'We have a problem here : ', swamp(ib) 8069 8349 WRITE(numout,*) 'resolution :', resolution_1, resolution(ib,2) 8070 8350 WRITE(numout,*) area_flood 8071 CALL ipslerr_p(3,'routing_ irrigmap','Problem with swamp...','','')8351 CALL ipslerr_p(3,'routing_floodmap','Problem with swamp...','','') 8072 8352 ENDIF 8073 8353 !!$ ! Compute a diagnostic of the map. … … 8083 8363 ! 8084 8364 ! 8085 8365 8086 8366 IF (printlev>=1) THEN 8087 IF ( init_irrig ) WRITE(numout,*) "Diagnostics irrigated :", MINVAL(irrigated), MAXVAL(irrigated)8088 8367 IF ( init_flood ) WRITE(numout,*) "Diagnostics floodplains :", MINVAL(floodplains), MAXVAL(floodplains) 8089 8368 IF ( init_swamp ) WRITE(numout,*) "Diagnostics swamp :", MINVAL(swamp), MAXVAL(swamp) … … 8112 8391 DEALLOCATE (latrel) 8113 8392 ! 8114 END SUBROUTINE routing_ irrigmap8393 END SUBROUTINE routing_floodmap 8115 8394 ! 8116 8395 !! ================================================================================================================================ … … 8167 8446 REAL(r_std), SAVE :: totw_out !! Sum of the water flow out to the routing scheme 8168 8447 !$OMP THREADPRIVATE(totw_out) 8169 REAL(r_std), SAVE :: totw_return !! 8448 REAL(r_std), SAVE :: totw_return !! 8170 8449 !$OMP THREADPRIVATE(totw_return) 8171 REAL(r_std), SAVE :: totw_irrig !! 8450 REAL(r_std), SAVE :: totw_irrig !! 8172 8451 !$OMP THREADPRIVATE(totw_irrig) 8173 REAL(r_std), SAVE :: totw_river !! 8452 REAL(r_std), SAVE :: totw_river !! 8174 8453 !$OMP THREADPRIVATE(totw_river) 8175 REAL(r_std), SAVE :: totw_coastal !! 8454 REAL(r_std), SAVE :: totw_coastal !! 8176 8455 !$OMP THREADPRIVATE(totw_coastal) 8177 8456 REAL(r_std) :: totarea !! Total area of basin (m^2) 8178 8457 REAL(r_std) :: area !! Total area of routing (m^2) 8179 INTEGER(i_std) :: ig !! 8458 INTEGER(i_std) :: ig !! 8180 8459 ! 8181 8460 ! Just to make sure we do not get too large numbers ! … … 8194 8473 totw_slow = zero 8195 8474 totw_lake = zero 8196 totw_pond = zero 8475 totw_pond = zero 8197 8476 totw_in = zero 8198 8477 !
Note: See TracChangeset
for help on using the changeset viewer.