Small corrections Linkup + improvement routing_reg_divbas

F90subroutines/routing_reg.f90

@@ -1089,9 +1089,9 @@ CONTAINS
     INTEGER(i_std), DIMENSION(nbi,nbj) :: basintmp !!
     REAL(r_std), DIMENSION(nbi,nbj) :: lontmp, lattmp !!
     !
-    INTEGER(i_std) :: il, jl, ill, jll !!
+    INTEGER(i_std) :: il, jl, ill, jll, jp !!
     INTEGER(i_std) :: ie, je !!
-    INTEGER(i_std) :: p, cnt, k, l, ic, ik !!
+    INTEGER(i_std) :: p, cnt, k, l, ic, ik, cut !!
     INTEGER(i_std) :: ip, isz, checksz !!
     INTEGER(i_std) :: ff(1) !!
     !
@@ -1123,6 +1123,8 @@ CONTAINS
     !
     INTEGER(i_std) :: main_loc(nbv,2), tri_loc(nbne,nbv,2)
     REAL(r_std) :: main_fac(nbv), tri_fac(nbne,nbv)
+    REAL(r_std) :: diff_main_fac(nbv)
+    REAL(r_std) :: half_output_fac
     REAL(r_std) :: tmptri_fac(nbne)  ! Flow accumulation of all neighbour points
     INTEGER(i_std) :: sortedtrifac(nbne) ! And sort of tmptri_fac(nbne)
     REAL(r_std) :: alltri_fac(nbv) ! Sort flow accumulation of all tributaries
@@ -1191,6 +1193,24 @@ CONTAINS
        lattmp(il,jl) = lat(il,jl)
     ENDDO
     !
+    ! Calculation of local flow accumulation
+    factmp(:,:) = 0
+    DO isz = 1, tsz
+       ip = tpts(ibas, isz, 1)
+       jp = tpts(ibas, isz, 2)
+       IF ( trip(ip,jp) .GT. 0 .AND. trip(ip,jp) .LT. 109 ) THEN
+         p = trip(ip,jp)
+         il = ip ; jl = jp
+         DO WHILE ( p .GT. 0 .AND. p .LT. 9 .AND. cnt .LT. nbi*nbj)
+           ill = il + inc(p,1)
+           jll = jl + inc(p,2)
+           il = ill ; jl = jll
+           factmp(il, jl) = factmp(il, jl) + 1
+           p = trip(il,jl)
+         ENDDO
+       ENDIF
+    END DO
+    !
     ! Print out information of grid box
     !
     IF ( checkgrid ) THEN
@@ -1329,183 +1349,42 @@ CONTAINS
           okpoint = .FALSE.
        ENDIF
        !
-       IF ( cnt .EQ. nbi*nbj) THEN
-          IF ( debug ) THEN
-             WRITE(numout,*) "Error: cnt .EQ. nbi*nbj "
-             WRITE(numout,*) "Point: ", il, jl
-             WRITE(numout,*) "Per: ", nbi, nbj
-             WRITE(fmt,"('(',I3,'I6)')") nbi
-             WRITE(fmtr,"('(',I3,'F8.1)')") nbi
-             !
-             WRITE(numout,*) '+++++++++++++++++++ TRIP IN DIVBAS ++++++++++++++++++++'
-             DO je=1,nbj
-                WRITE(numout,fmt) triptemp(1:nbi,je)
-             ENDDO
-             !
-             WRITE(numout,*) '+++++++++++++++++++ BASIN IDs IN DIVBAS ++++++++++++++++++++'
-             DO je=1,nbj
-                WRITE(numout,fmt) basintmp(1:nbi,je)
-             ENDDO
-             !
-             WRITE(numout,*) '+++++++++++++++++++ LONGITUDE IN DIVBAS ++++++++++++++++++++'
-             DO je=1,nbj
-                WRITE(numout,fmtr) lontmp(1:nbi,je)
-             ENDDO
-             !
-             WRITE(numout,*) '+++++++++++++++++++ LATITUDE IN DIVBAS ++++++++++++++++++++'
-             DO je=1,nbj
-                WRITE(numout,fmtr) lattmp(1:nbi,je)
-             ENDDO
-             !
-             WRITE(numout,*) '+++++++++++++++++++ FAC IN DIVBAS ++++++++++++++++++++'
-             DO je=1,nbj
-                WRITE(numout,fmtr) factmp(1:nbi,je)
-             ENDDO
-             CALL FLUSH(numout)
-          ENDIF
-          CALL ipslerr_p(3,'routing_reg_divbas','We could not route point','','')
-       ENDIF
-       !
     ENDDO
     IF ( debug ) WRITE (numout,*) 'Length of mainstream: ', cnt
     !
     ! 3.0 Sort the flow accumulation of the tributaries and find the
     ! divided location.
     !
-    alltri_fac(:) = flag
-    alltri_loc1(:) = 0
-    alltri_loc2(:) = 0
-    divloc(:,:) = 0
-    !
-    ip = 0
-    DO k = 1, cnt
-       ! Actually, we should only DO l = 1, nbne-2
-       ! Because, there are usually 2 neighbour points belong to main stream
-       ! and maximum (nbne-2) neighbour points can be tributaries.
-       ! DO l = 1, nbne-2
-       DO l = 1, nbne
-          IF ( tri_fac(l,k) .NE. flag ) THEN
-             ip = ip + 1
-             alltri_fac(ip) = tri_fac(l,k)
-             ! Attention: here we store location of each tributary in main stream
-             ! ( k <= cnt ) and its direction compared with main stream ( l <= nbne )
-             alltri_loc1(ip) = k
-             alltri_loc2(ip) = l
-          ENDIF
-       ENDDO
-    ENDDO
-    IF ( debug ) WRITE (numout,*) 'Number of tributaries: ', ip
-    ! If we have at least one tributary
-    IF ( ip .GT. 0 ) THEN
-       ! Original output
-       oic = 0
-       !
-       sortedallfac(:) = 0
-       DO k = 1, ip
-          ff = MAXLOC(alltri_fac)
-          sortedallfac(k) = ff(1)
-          alltri_fac(ff(1)) = flag
-       ENDDO
-       ! Get the points in the mainstream
-       tmpmain_fac(:) = flag
-       tmpmain_loc(:,:) = 0
-       numtri(:) = 0
-       usetri_loc(:,:,:) = 0
-       ! Using l to count the actual number of points in main stream need to use
-       ! for dividing (not always equal 4 as ideal case).
-       l = 0
-       ! According to Pfafstetter coding, we only look for 4 greatest
-       ! tributaries:
-       DO k = 1, 4
-          !
-          IF ( sortedallfac(k) .NE. 0 ) THEN
-             !
-             ik = alltri_loc1(sortedallfac(k))
-             il = alltri_loc2(sortedallfac(k))
-             ! Starting store location and flow accumulation of points in main stream:
-             IF ( COUNT ( tmpmain_fac(:) .EQ. main_fac(ik) ) .EQ. 0 ) THEN
-                l = l + 1
-                tmpmain_loc(l,1) = main_loc(ik,1)
-                tmpmain_loc(l,2) = main_loc(ik,2)
-                tmpmain_fac(l) = main_fac(ik)
-                !
-                numtri(l) = 1
-                usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1)
-                usetri_loc(numtri(l),l,2) = tri_loc(il,ik,2)
-                !
-             ELSE
-                ! If this point in main stream have more than 1 tributary
-                ! belongs to top 4 greatest tributaries
-                numtri(l) = numtri(l) + 1
-                usetri_loc(numtri(l),l,1) = tri_loc(il,ik,1)
-                usetri_loc(numtri(l),l,2) = tri_loc(il,ik,2)
-                !
-             ENDIF
-          ENDIF
-       ENDDO
-       ! Sort the mainstream points
-       sortedmainfac(:) = 0
-       DO k = 1, l
-          ff = MAXLOC(tmpmain_fac)
-          sortedmainfac(k) = ff(1)
-          tmpmain_fac(ff(1)) = flag
-       ENDDO
-       !
-       ic = 0
-       DO k = l, 1, -1
-          IF ( tmpmain_loc(sortedmainfac(k),1) .NE. 0 ) THEN
-             ic = ic + 1
-             divloc(ic,1) = tmpmain_loc(sortedmainfac(k),1)
-             divloc(ic,2) = tmpmain_loc(sortedmainfac(k),2)
-             !
-             new_outbas(ic) = ic + 1
-             !
-          ENDIF
-       ENDDO
-       ! If the last divide point is not the original outlet
-       IF ( divloc(ic,1) .NE. iloc .OR. divloc(ic,2) .NE. jloc ) THEN
-          ic = ic + 1
-          divloc(ic,1) = iloc
-          divloc(ic,2) = jloc
-       ENDIF
-       !
-       new_outbas(ic) = tout
-       ! Number of mainstream sub-basin
-       oic = ic
-       !
-       DO k = l, 1, -1
-          !
-          !DO il = 1, nbne-2
-          ! Maximum 4 tributaries:
-          DO il = 1, 4
-             IF ( usetri_loc(il,sortedmainfac(k),1) .NE. 0 ) THEN
-                ic = ic + 1
-                divloc(ic,1) = usetri_loc(il,sortedmainfac(k),1)
-                divloc(ic,2) = usetri_loc(il,sortedmainfac(k),2)
-                !
-                !new_outbas(ic) = new_outbas(l-k+1)
-                new_outbas(ic) = l-k+1
-                !
-             ENDIF
-          ENDDO
-       ENDDO
-       ! Save number of dividing points
-       l = ic
-       !
-    ! If we don't have any tributary
+    ! Here the process is simplified we just cut the basin in 2
+    ! over the main river (defined by local fac) 
+    ! Where the local flow accumulation is half the local fac at the output
+    !
+    oic = 2
+    l = 2
+    ! Find the main_fac(cnt)
+    diff_main_fac(:)=9999
+    ! half_output_fac is half the local fac value for the output 
+    half_output_fac = main_fac(1)/2
+    ! Get the difference between local fac and half_output_fac for each point
+    ! of the main river
+    DO ik=1,cnt
+        diff_main_fac(ik) = abs(main_fac(ik)-half_output_fac)
+    END DO
+    ff = MINLOC(diff_main_fac(:))
+    IF ((ff(1) .EQ. 1) .OR. (ff(1) .EQ. cnt)) THEN
+        ! In order to avoid limit cases
+        cut = cnt/2 + 1
     ELSE
-       oic = 2
-       l = 2
-       cnt = cnt/2 + 1
-       divloc(1,1) = main_loc(cnt,1)
-       divloc(1,2) = main_loc(cnt,2)
-       divloc(2,1) = iloc
-       divloc(2,2) = jloc
-       !
-       new_outbas(1) = 2
-       new_outbas(2) = tout
-       !
-    ENDIF
+        cut = ff(1)
+    END IF
+    divloc(1,1) = main_loc(cut,1)
+    divloc(1,2) = main_loc(cut,2)
+    divloc(2,1) = iloc
+    divloc(2,2) = jloc
+    !
+    new_outbas(1) = 2
+    new_outbas(2) = tout
+    !
     IF ( debug ) WRITE (numout,*) 'Number of new sub-basin: ', l, oic
     !
     ! Now, cut it ! Cut the sub-basin ! Release The Kraken !
@@ -1513,16 +1392,7 @@ CONTAINS
     new_nb = l
     !
     DO ik = 1, l
-       ! A small trick here
-       IF ( oic .EQ. 2 .AND. l .EQ. 2 ) THEN
-          k = ik
-       ELSE
-          IF ( ik .LE. (l - oic) ) THEN
-             k = ik + oic
-          ELSE
-             k = ik - (l - oic)
-          ENDIF
-       ENDIF
+       k = ik
        ! I'm so stupid careful here with this IF
        IF ( divloc(k,1) .NE. 0 .AND. divloc(k,2) .NE. 0 ) THEN
           !
@@ -1771,7 +1641,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        !
        basin_area(ib,ij) = zero
        basin_cg(ib,ij,:) = zero
-       basin_hierarchy(ib,ij) = zero
+       basin_hierarchy(ib,ij) = 0
        basin_fac(ib,ij) = zero
        basin_orog_mean(ib,ij) = zero
        basin_orog_min(ib,ij) = 99999
@@ -1807,6 +1677,8 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
           ! - Take the minimum value within the basin
           ! - Take the value at the outflow point
           ! Probably taking the value of the outflow point is the best solution.
+          !
+          
           !
           SELECT CASE (hierar_method)
              !
@@ -1818,10 +1690,13 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
                      & fac_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
              CASE("MINI")
                 ! The smallest value of the basin
+                IF (fac_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .GT. basin_fac(ib,ij)) THEN
+                   basin_fac(ib,ij) = fac_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
+                END IF
                 IF ( hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) .LT. basin_hierarchy(ib,ij)) THEN
                    basin_hierarchy(ib,ij) = hierarchy_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
                    ! Trung: We should take the "fac" value at the same point
-                   basin_fac(ib,ij) = fac_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
+                   !basin_fac(ib,ij) = fac_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
                 ENDIF
              CASE("OUTP")
                 ! Value at the outflow point
@@ -2078,6 +1953,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
           found = 0
           IF (  outflow_grid(sp,sb) == 0 ) THEN
               found = 1
+              IF (debug) WRITE(numout,*) sp, sb, "Linkup 1.0 -- Flow out of Halo zone"
           ELSE
               found = 0
           END IF
@@ -2091,6 +1967,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                 IF ( outflow_basin(sp,sb) == bop ) THEN
                    found = 1
                    solved(sp,1) = solved(sp,1) + 1
+                   IF (debug) WRITE(numout,*) sp, sb, "flows in the same grid !"
                 ELSE
                    WRITE(numout,*) sp, sb, "flows in the same grid but has an error !"
                 ENDIF
@@ -2100,16 +1977,19 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                 ! Nothing to do but just remember it is done.
                 found = 1
                 solved(sp,1) = solved(sp,1) + 1
+                IF (debug) WRITE(numout,*) sp,sb,"is a return flow"
              ELSE IF ( outflow_grid(sp,sb) .EQ. -2 ) THEN
                 ! Coastal flow
                 ! Nothing to do but just remember it is done.
                 found = 1
                 solved(sp,1) = solved(sp,1) + 1
+                IF (debug) WRITE(numout,*) sp,sb,"is a coastal flow"
              ELSE IF ( outflow_grid(sp,sb) .EQ. -1 ) THEN
                 ! River flow
                 ! Nothing to do but just remember it is done.
                 found = 1
                 solved(sp,1) = solved(sp,1) + 1
+                IF (debug) WRITE(numout,*) sp,sb,"is a river outflow"
              ENDIF
           END IF
           IF ( found .EQ. 0 ) THEN
@@ -2123,18 +2003,20 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                   &                    nbpt, nwbas, inp, basin_count, basin_id, basin_hierarchy, basin_fac, &
                   &                    basin_flowdir, nbcoastal, coastal_basin, bop, bopqual)
              !
-             IF ( bop .LT. undef_int .AND. bop .NE. sb) THEN
+             IF ( bop .LT. undef_int) THEN
                 !
                 CALL routing_updateflow(sp, sb, inp, bop, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
                      &                  inflow_number, inflow_grid, inflow_basin)
                 IF ( outflow_basin(sp,sb) == bop ) THEN
                    solved(sp,1) = solved(sp,1) + 1
                    found = 1
+                   IF (debug) WRITE (numout,*) sp, sb, "Flows in the 1st approx grid point:", inp, bop
                 ENDIF
              !
              ENDIF
              !
              IF ( found == 0 ) THEN
+                IF (debug) WRITE (numout,*) "Establishing the list of neighbours"
                 ! Organize the location of the neighbours to visit by order of priority
                 ! first the outflow grid then 2 by 2 till the opposite side (by +1/-1 - +2/-2 ...)
                 ! if NbNeighb is odd we have to had the opposite neighbour
@@ -2247,6 +2129,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                          IF ( outflow_basin(sp,sb) == bop ) THEN
                             solved(sp,2) = solved(sp,2) + 1
                             found = 1
+                            IF (debug) WRITE (numout,*) sp, sb,"->Sol. in neighbours,output found at level:",nb,"dop,bop=",dop,bop
                          ELSE
                             nb = nb+1
                          ENDIF
@@ -2262,10 +2145,35 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
              !   ! Do like first one
              !  ! loc(NbNeighb)
              ENDIF
+             ! Opposite point
+             IF (found .EQ. 0) THEN
+               IF ( modulo(NbNeighb,2) == 0 ) THEN
+                 dop = neighbours_g(sp,order_ref(NbNeighb))
+                 IF (dop .GT. 0) THEN
+                   CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb),basin_hierarchy(sp,sb), basin_fac(sp,sb), &
+                                &                    basin_flowdir(sp,sb),invented_basins, &
+                                &                    nbpt, nwbas, dop ,basin_count, basin_id, basin_hierarchy, basin_fac, &
+                                &                    basin_flowdir, nbcoastal,coastal_basin, bop, bopqual)
+                   IF ( dop > 0 .AND. dop < undef_int .AND. bop < undef_int ) THEN
+                           IF (debug) WRITE (numout,*) sp, sb,"Opposite neighbours WRITE", dop
+                           !
+                           CALL routing_updateflow(sp, sb, dop, bop, nbpt, nwbas,inflowmax, outflow_grid, outflow_basin, &
+                                &                  inflow_number, inflow_grid,inflow_basin)
+                           !
+                           IF ( outflow_basin(sp,sb) == bop ) THEN
+                              found = 1
+                              IF (debug) WRITE (numout,*) sp, sb,"Sol. in opposite neighbours dop,bop=",dop,bop
+                           ENDIF
+                           !
+                   ENDIF
+                 ENDIF
+               END IF
+             END IF
              !
              ! Looking for a solution in the grid -> HTU with a similar hierarchy or lowest hierarchy
              !
              IF ( found == 0 ) THEN
+                IF (debug) WRITE (numout,*) "Looking for a solution in the grid"
                 sbint = undef_int
                 DO sba=1,basin_count(sp)
                    IF ( sba .NE. sb ) THEN
@@ -2315,6 +2223,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                         &                  inflow_number, inflow_grid,inflow_basin)
                    IF (outflow_basin(sp,sb) == sbint) THEN
                       found = 1
+                      IF (debug) WRITE (numout,*) sp, sb, "Lowest basin hierarchy in the grid file"
                    END IF
                 ENDIF
              ENDIF
@@ -2372,9 +2281,44 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                 !        But I need the model works now !!! => so, come back
                 !        here later !
                 !
-                ! Coastal flow or river flow is both ok here
+                ! If debug -> More detailled evaluation of situation
+                IF (debug) THEN
+                   WRITE (numout,*) "initial outflowgrid", outflow_grid(sp,sb)
+                   WRITE (numout,*) "initial hierarchy",basin_hierarchy(sp,sb) 
+                   WRITE (numout,*) "initial fac", basin_fac(sp,sb)
+                   !
+                   WRITE (numout,*) "Neighbours with basin_id = ", basin_id(sp,sb)
+                   DO nb=1,NbNeighb
+                     IF (neighbours_g(sp,nb) .GT. 0) THEN
+                       dop = neighbours_g(sp,nb)
+                       DO sba=1,basin_count(dop)
+                         IF (basin_id(dop,sba) ==basin_id(sp,sb)) THEN
+                           IF (basin_hierarchy(dop,sba) .LT. basin_hierarchy(sp,sb)) THEN
+                             IF (basin_fac(dop,sba) .GT. basin_fac(sp,sb)) THEN
+                               val =basin_hierarchy(dop,sba)
+                               dm1 = basin_fac(dop,sba)
+                               WRITE (numout,*) "NBPT, minimum hierarchy for basin(npt,hierarch,fac)", dop,val,dm1 
+                             END IF
+                           END IF
+                         END IF
+                       END DO
+                     END IF
+                   END DO
+                   !
+                   val = 999999998
+                   DO sba=1,basin_count(sp)
+                     IF (basin_id(sp,sba) ==basin_id(sp,sb)) THEN
+                       IF (basin_hierarchy(sp,sba) .LT. basin_hierarchy(sp,sb)) THEN
+                          val = basin_hierarchy(sp,sba)
+                          dm1 = basin_fac(dop,sba)
+                          WRITE (numout,*) "SAME PT, minimum hierarchy for basin",sp,val,dm1
+                       END IF
+                     END IF
+                   END DO
+                END IF
                 outflow_grid(sp,sb) = -2
                 basin_hierarchy(sp,sb) = 0.00
+                WRITE (numout,*) "Linkup : Made a NEW OUTLET at sp & sb:",sp,sb 
              ENDIF
              IF ( outflow_basin(sp,sb) < undef_int ) solved(sp,4) = solved(sp,4)+1
           ENDIF