diff --git a/main/HYDRO/MOD_Catch_BasinNetwork.F90 b/main/HYDRO/MOD_Catch_BasinNetwork.F90 index 2a97d386..2119e718 100644 --- a/main/HYDRO/MOD_Catch_BasinNetwork.F90 +++ b/main/HYDRO/MOD_Catch_BasinNetwork.F90 @@ -14,6 +14,9 @@ MODULE MOD_Catch_BasinNetwork ! -- instances -- integer :: numbasin integer, allocatable :: basinindex(:) + + integer :: numrivmth + integer, allocatable :: rivermouth(:) integer :: numbsnhru type(subset_type) :: basin_hru @@ -28,6 +31,7 @@ MODULE MOD_Catch_BasinNetwork integer, allocatable :: paddr (:) integer, allocatable :: ndata (:) integer, allocatable :: ipush (:) + integer, allocatable :: isdrv (:) CONTAINS final :: basin_pushdata_free_mem END type basin_pushdata_type @@ -231,6 +235,20 @@ SUBROUTINE build_basin_network () ENDIF ENDDO + allocate (rivermouth (totalnumbasin)) + ithis = 0 + DO i = totalnumbasin, 1, -1 + j = basindown(b_up2down(i)) + IF (j <= 0) THEN + ithis = ithis + 1 + rivermouth(b_up2down(i)) = ithis + ELSE + rivermouth(b_up2down(i)) = rivermouth(j) + ENDIF + ENDDO + + numrivmth = ithis + deallocate (b_up2down) deallocate (nups_all ) deallocate (orderbsn ) @@ -238,6 +256,8 @@ SUBROUTINE build_basin_network () ENDIF + CALL mpi_bcast (numrivmth, 1, MPI_INTEGER, p_address_master, p_comm_glb, p_err) + ! 3-3: send basin index to workers IF (p_is_master) THEN @@ -266,6 +286,10 @@ SUBROUTINE build_basin_network () CALL mpi_send (icache, nbasin, MPI_INTEGER, p_address_worker(iworker), & mpi_tag_data, p_comm_glb, p_err) + icache = rivermouth(bindex) + CALL mpi_send (icache, nbasin, MPI_INTEGER, p_address_worker(iworker), & + mpi_tag_data, p_comm_glb, p_err) + nhru_in_bsn = nhru_all(bindex) CALL mpi_send (nhru_in_bsn, nbasin, MPI_INTEGER, p_address_worker(iworker), & mpi_tag_data, p_comm_glb, p_err) @@ -297,6 +321,10 @@ SUBROUTINE build_basin_network () CALL mpi_recv (basindown, numbasin, MPI_INTEGER, p_address_master, & mpi_tag_data, p_comm_glb, p_stat, p_err) + allocate (rivermouth (numbasin)) + CALL mpi_recv (rivermouth, numbasin, MPI_INTEGER, p_address_master, & + mpi_tag_data, p_comm_glb, p_stat, p_err) + allocate (nhru_in_bsn (numbasin)) CALL mpi_recv (nhru_in_bsn, numbasin, MPI_INTEGER, p_address_master, & mpi_tag_data, p_comm_glb, p_stat, p_err) @@ -592,7 +620,7 @@ END SUBROUTINE build_basin_network ! ---------- - SUBROUTINE worker_push_data_real8 (send_pointer, recv_pointer, accum, vec_send, vec_recv) + SUBROUTINE worker_push_data_real8 (send_pointer, recv_pointer, vec_send, vec_recv) USE MOD_Precision USE MOD_SPMD_Task @@ -600,7 +628,6 @@ SUBROUTINE worker_push_data_real8 (send_pointer, recv_pointer, accum, vec_send, type(basin_pushdata_type) :: send_pointer type(basin_pushdata_type) :: recv_pointer - logical, intent(in) :: accum real(r8), intent(in) :: vec_send(:) real(r8), intent(inout) :: vec_recv(:) @@ -619,14 +646,7 @@ SUBROUTINE worker_push_data_real8 (send_pointer, recv_pointer, accum, vec_send, IF (p_is_worker) THEN IF (send_pointer%nself > 0) THEN - IF (.not. accum) THEN - vec_recv(recv_pointer%iself) = vec_send(send_pointer%iself) - ELSE - DO i = 1, send_pointer%nself - vec_recv(recv_pointer%iself(i)) = & - vec_recv(recv_pointer%iself(i)) + vec_send(send_pointer%iself(i)) - ENDDO - ENDIF + vec_recv(recv_pointer%iself) = vec_send(send_pointer%iself) ENDIF #ifdef USEMPI @@ -675,17 +695,8 @@ SUBROUTINE worker_push_data_real8 (send_pointer, recv_pointer, accum, vec_send, ENDIF IF (recv_pointer%nproc > 0) THEN - CALL mpi_waitall(recv_pointer%nproc, req_recv, MPI_STATUSES_IGNORE, p_err) - - IF (accum) THEN - DO i = 1, ndatarecv - vec_recv(recv_pointer%ipush(i)) = & - vec_recv(recv_pointer%ipush(i)) + recvcache(i) - ENDDO - ELSE - vec_recv(recv_pointer%ipush) = recvcache - ENDIF + vec_recv(recv_pointer%ipush) = recvcache ENDIF IF (send_pointer%nproc > 0) THEN @@ -705,7 +716,7 @@ SUBROUTINE worker_push_data_real8 (send_pointer, recv_pointer, accum, vec_send, END SUBROUTINE worker_push_data_real8 ! ---------- - SUBROUTINE worker_push_data_int32 (send_pointer, recv_pointer, accum, vec_send, vec_recv) + SUBROUTINE worker_push_data_int32 (send_pointer, recv_pointer, vec_send, vec_recv) USE MOD_Precision USE MOD_SPMD_Task @@ -713,7 +724,6 @@ SUBROUTINE worker_push_data_int32 (send_pointer, recv_pointer, accum, vec_send, type(basin_pushdata_type) :: send_pointer type(basin_pushdata_type) :: recv_pointer - logical, intent(in) :: accum integer, intent(in) :: vec_send(:) integer, intent(inout) :: vec_recv(:) @@ -732,14 +742,7 @@ SUBROUTINE worker_push_data_int32 (send_pointer, recv_pointer, accum, vec_send, IF (p_is_worker) THEN IF (send_pointer%nself > 0) THEN - IF (.not. accum) THEN - vec_recv(recv_pointer%iself) = vec_send(send_pointer%iself) - ELSE - DO i = 1, send_pointer%nself - vec_recv(recv_pointer%iself(i)) = & - vec_recv(recv_pointer%iself(i)) + vec_send(send_pointer%iself(i)) - ENDDO - ENDIF + vec_recv(recv_pointer%iself) = vec_send(send_pointer%iself) ENDIF #ifdef USEMPI @@ -788,17 +791,8 @@ SUBROUTINE worker_push_data_int32 (send_pointer, recv_pointer, accum, vec_send, ENDIF IF (recv_pointer%nproc > 0) THEN - CALL mpi_waitall(recv_pointer%nproc, req_recv, MPI_STATUSES_IGNORE, p_err) - - IF (accum) THEN - DO i = 1, ndatarecv - vec_recv(recv_pointer%ipush(i)) = & - vec_recv(recv_pointer%ipush(i)) + recvcache(i) - ENDDO - ELSE - vec_recv(recv_pointer%ipush) = recvcache - ENDIF + vec_recv(recv_pointer%ipush) = recvcache ENDIF IF (send_pointer%nproc > 0) THEN @@ -1003,6 +997,16 @@ SUBROUTINE worker_push_subset_data (send_pointer, recv_pointer, & END SUBROUTINE worker_push_subset_data + ! --------- + SUBROUTINE basin_network_final () + + IMPLICIT NONE + + IF (allocated(basinindex)) deallocate(basinindex) + IF (allocated(rivermouth)) deallocate(rivermouth) + + END SUBROUTINE basin_network_final + ! --------- SUBROUTINE basin_pushdata_free_mem (this) @@ -1013,6 +1017,7 @@ SUBROUTINE basin_pushdata_free_mem (this) IF (allocated(this%paddr)) deallocate(this%paddr) IF (allocated(this%ndata)) deallocate(this%ndata) IF (allocated(this%ipush)) deallocate(this%ipush) + IF (allocated(this%isdrv)) deallocate(this%isdrv) END SUBROUTINE basin_pushdata_free_mem diff --git a/main/HYDRO/MOD_Catch_Hist.F90 b/main/HYDRO/MOD_Catch_Hist.F90 index a62e16c4..2853d1aa 100644 --- a/main/HYDRO/MOD_Catch_Hist.F90 +++ b/main/HYDRO/MOD_Catch_Hist.F90 @@ -164,7 +164,7 @@ SUBROUTINE hist_basin_out (file_hist, idate) END WHERE ENDIF - CALL worker_push_data (iam_bsn, iam_elm, .false., a_wdsrf_bsn, a_wdsrf_elm) + CALL worker_push_data (iam_bsn, iam_elm, a_wdsrf_bsn, a_wdsrf_elm) CALL vector_write_basin (& file_hist_basin, a_wdsrf_elm, numelm, totalnumelm, 'wdsrf_bsn', 'basin', elm_data_address, & @@ -177,7 +177,7 @@ SUBROUTINE hist_basin_out (file_hist, idate) END WHERE ENDIF - CALL worker_push_data (iam_bsn, iam_elm, .false., a_veloc_riv, a_veloc_elm) + CALL worker_push_data (iam_bsn, iam_elm, a_veloc_riv, a_veloc_elm) CALL vector_write_basin (& file_hist_basin, a_veloc_elm, numelm, totalnumelm, 'veloc_riv', 'basin', elm_data_address, & @@ -190,7 +190,7 @@ SUBROUTINE hist_basin_out (file_hist, idate) END WHERE ENDIF - CALL worker_push_data (iam_bsn, iam_elm, .false., a_discharge, a_dschg_elm) + CALL worker_push_data (iam_bsn, iam_elm, a_discharge, a_dschg_elm) CALL vector_write_basin (& file_hist_basin, a_dschg_elm, numelm, totalnumelm, 'discharge', 'basin', elm_data_address, & diff --git a/main/HYDRO/MOD_Catch_LateralFlow.F90 b/main/HYDRO/MOD_Catch_LateralFlow.F90 index 1ca274d8..e1418297 100644 --- a/main/HYDRO/MOD_Catch_LateralFlow.F90 +++ b/main/HYDRO/MOD_Catch_LateralFlow.F90 @@ -234,7 +234,7 @@ SUBROUTINE lateral_flow (deltime) #ifdef RangeCheck IF (p_is_worker .and. (p_iam_worker == 0)) THEN write(*,'(/,A)') 'Checking Lateral Flow Variables ...' - write(*,'(A,F12.5,A)') 'River Lake Flow average timestep: ', & + write(*,'(A,F12.5,A)') 'River Lake Flow minimum average timestep: ', & dt_average/nsubstep, ' seconds' ENDIF @@ -265,7 +265,7 @@ SUBROUTINE lateral_flow (deltime) CALL mpi_allreduce (MPI_IN_PLACE, toldis, 1, MPI_REAL8, MPI_SUM, p_comm_worker, p_err) #endif IF (p_iam_worker == 0) THEN - write(*,'(A,F10.2,A,ES10.3,A,ES10.3,A)') 'Total surface water error: ', dtolw, & + write(*,'(A,F12.2,A,ES8.1,A,ES10.3,A)') 'Total surface water error: ', dtolw, & '(m^3) in area ', landarea, '(m^2), discharge ', toldis, '(m^3)' ENDIF @@ -275,7 +275,7 @@ SUBROUTINE lateral_flow (deltime) CALL mpi_allreduce (MPI_IN_PLACE, dtolw, 1, MPI_REAL8, MPI_SUM, p_comm_worker, p_err) #endif IF (p_iam_worker == 0) THEN - write(*,'(A,F10.2,A,ES10.3,A)') 'Total ground water error: ', dtolw, & + write(*,'(A,F12.2,A,ES8.1,A)') 'Total ground water error: ', dtolw, & '(m^3) in area ', landarea, '(m^2)' ENDIF ENDIF @@ -291,6 +291,7 @@ SUBROUTINE lateral_flow_final () CALL river_lake_network_final () CALL subsurface_network_final () + CALL basin_network_final () #ifdef CoLMDEBUG IF (allocated(patcharea)) deallocate(patcharea) diff --git a/main/HYDRO/MOD_Catch_RiverLakeFlow.F90 b/main/HYDRO/MOD_Catch_RiverLakeFlow.F90 index bc4c97ad..e4c9439d 100644 --- a/main/HYDRO/MOD_Catch_RiverLakeFlow.F90 +++ b/main/HYDRO/MOD_Catch_RiverLakeFlow.F90 @@ -36,6 +36,7 @@ MODULE MOD_Catch_RiverLakeFlow SUBROUTINE river_lake_flow (dt) USE MOD_SPMD_Task + USE MOD_Utils USE MOD_Catch_BasinNetwork USE MOD_Catch_HillslopeNetwork USE MOD_Catch_RiverLakeNetwork @@ -47,7 +48,9 @@ SUBROUTINE river_lake_flow (dt) real(r8), intent(in) :: dt ! Local Variables - integer :: hs, he, i, j + integer :: hs, he, i, j, nsync + real(r8) :: dt_this + logical :: allfinished real(r8), allocatable :: wdsrf_bsn_ds(:) real(r8), allocatable :: veloc_riv_ds(:) @@ -65,12 +68,12 @@ SUBROUTINE river_lake_flow (dt) real(r8) :: bedelv_fc, height_up, height_dn real(r8) :: vwave_up, vwave_dn, hflux_up, hflux_dn, mflux_up, mflux_dn real(r8) :: totalvolume, loss, friction, dvol, nextl, nexta, nextv, ddep - real(r8) :: dt_res, dt_this - logical, allocatable :: mask(:) + real(r8), allocatable :: dt_res(:), dt_all(:), dt_tmp(:) + logical, allocatable :: hmask(:), bsnfilter(:), linkfilter(:), syncfilter(:) IF (p_is_worker) THEN - + ! update water depth in basin by aggregating water depths in patches DO i = 1, numbasin hs = basin_hru%substt(i) @@ -115,23 +118,48 @@ SUBROUTINE river_lake_flow (dt) allocate (sum_hflux_riv (numbasin)) allocate (sum_mflux_riv (numbasin)) allocate (sum_zgrad_riv (numbasin)) + allocate (bsnfilter (numbasin)) ENDIF + IF (numbsnlink > 0) allocate (linkfilter (numbsnlink)) + + allocate (dt_res (numrivmth)) + allocate (dt_all (numrivmth)) + allocate (syncfilter (numrivmth)) + ntimestep_riverlake = 0 - dt_res = dt - DO WHILE (dt_res > 0) + + WHERE (riv_multi_proc) + dt_res = dt + ELSEWHERE + dt_res = 0 + ENDWHERE + + DO i = 1, numbasin + dt_res(rivermouth(i)) = dt + ENDDO + + allfinished = .false. + DO WHILE (.not. allfinished) ntimestep_riverlake = ntimestep_riverlake + 1 - + DO i = 1, numbasin - sum_hflux_riv(i) = 0. - sum_mflux_riv(i) = 0. - sum_zgrad_riv(i) = 0. + bsnfilter(i) = dt_res(rivermouth(i)) > 0 + IF (bsnfilter(i)) THEN + sum_hflux_riv(i) = 0. + sum_mflux_riv(i) = 0. + sum_zgrad_riv(i) = 0. + ENDIF ENDDO - CALL worker_push_data (river_iam_dn, river_iam_up, .false., wdsrf_bsn, wdsrf_bsn_ds) - CALL worker_push_data (river_iam_dn, river_iam_up, .false., veloc_riv, veloc_riv_ds) - CALL worker_push_data (river_iam_dn, river_iam_up, .false., momen_riv, momen_riv_ds) + DO i = 1, numbsnlink + linkfilter(i) = dt_res(linkrivmth(i)) > 0 + ENDDO + + CALL pull_from_downstream (wdsrf_bsn, wdsrf_bsn_ds, bsnfilter, linkfilter) + CALL pull_from_downstream (veloc_riv, veloc_riv_ds, bsnfilter, linkfilter) + CALL pull_from_downstream (momen_riv, momen_riv_ds, bsnfilter, linkfilter) ! velocity in ocean or inland depression is assumed to be 0. IF (numbasin > 0) THEN @@ -140,9 +168,12 @@ SUBROUTINE river_lake_flow (dt) END WHERE ENDIF - dt_this = dt_res + dt_all(:) = dt_res(:) DO i = 1, numbasin + + IF (.not. bsnfilter(i)) CYCLE + IF (riverdown(i) >= 0) THEN IF (riverdown(i) > 0) THEN @@ -260,6 +291,7 @@ SUBROUTINE river_lake_flow (dt) ENDIF IF ((lake_id(i) < 0) .and. (hflux_fc(i) < 0)) THEN + dt_this = dt_all(rivermouth(i)) hflux_fc(i) = & max(hflux_fc(i), (height_up-height_dn) / dt_this * sum(hillslope_basin(i)%area, & mask = hillslope_basin(i)%hand <= wdsrf_bsn(i) + handmin(i))) @@ -273,16 +305,21 @@ SUBROUTINE river_lake_flow (dt) IF (numbasin > 0) THEN hflux_fc = - hflux_fc; mflux_fc = - mflux_fc; zgrad_dn = - zgrad_dn ENDIF + + CALL push_to_downstream (hflux_fc, sum_hflux_riv, bsnfilter, linkfilter) + CALL push_to_downstream (mflux_fc, sum_mflux_riv, bsnfilter, linkfilter) + CALL push_to_downstream (zgrad_dn, sum_zgrad_riv, bsnfilter, linkfilter) - CALL worker_push_data (river_iam_up, river_iam_dn, .true., hflux_fc, sum_hflux_riv) - CALL worker_push_data (river_iam_up, river_iam_dn, .true., mflux_fc, sum_mflux_riv) - CALL worker_push_data (river_iam_up, river_iam_dn, .true., zgrad_dn, sum_zgrad_riv) - IF (numbasin > 0) THEN hflux_fc = - hflux_fc; mflux_fc = - mflux_fc; zgrad_dn = - zgrad_dn ENDIF DO i = 1, numbasin + + IF (.not. bsnfilter(i)) CYCLE + + dt_this = dt_all(rivermouth(i)) + ! constraint 1: CFL condition (only for rivers) IF (lake_id(i) == 0) THEN IF ((veloc_riv(i) /= 0.) .or. (wdsrf_bsn(i) > 0.)) THEN @@ -314,25 +351,38 @@ SUBROUTINE river_lake_flow (dt) abs(momen_riv(i) * riverarea(i) / (sum_mflux_riv(i)-sum_zgrad_riv(i)))) ENDIF ENDIF + + dt_all(rivermouth(i)) = min(dt_this, dt_all(rivermouth(i))) + ENDDO #ifdef USEMPI - CALL mpi_allreduce (MPI_IN_PLACE, dt_this, 1, MPI_REAL8, MPI_MIN, p_comm_worker, p_err) + syncfilter = riv_multi_proc .and. (dt_res > 0) + nsync = count(syncfilter) + IF (nsync > 0) THEN + allocate (dt_tmp (nsync)) + dt_tmp = pack(dt_all, mask = syncfilter) + CALL mpi_allreduce (MPI_IN_PLACE, dt_tmp, nsync, MPI_REAL8, MPI_MIN, p_comm_worker, p_err) + CALL unpack_inplace(dt_tmp, syncfilter, dt_all) + deallocate(dt_tmp) + ENDIF #endif DO i = 1, numbasin + IF (.not. bsnfilter(i)) CYCLE + IF (lake_id(i) <= 0) THEN ! rivers or lake catchments hs = basin_hru%substt(i) he = basin_hru%subend(i) - allocate (mask (hillslope_basin(i)%nhru)) + allocate (hmask (hillslope_basin(i)%nhru)) totalvolume = sum((wdsrf_bsn(i) + handmin(i) - hillslope_basin(i)%hand) & * hillslope_basin(i)%area, & mask = wdsrf_bsn(i) + handmin(i) >= hillslope_basin(i)%hand) - totalvolume = totalvolume - sum_hflux_riv(i) * dt_this + totalvolume = totalvolume - sum_hflux_riv(i) * dt_all(rivermouth(i)) IF (totalvolume < VOLUMEMIN) THEN DO j = 1, hillslope_basin(i)%nhru @@ -344,12 +394,12 @@ SUBROUTINE river_lake_flow (dt) wdsrf_bsn(i) = 0 ELSE - dvol = sum_hflux_riv(i) * dt_this + dvol = sum_hflux_riv(i) * dt_all(rivermouth(i)) IF (dvol > VOLUMEMIN) THEN DO WHILE (dvol > VOLUMEMIN) - mask = hillslope_basin(i)%hand < wdsrf_bsn(i) + handmin(i) - nextl = maxval(hillslope_basin(i)%hand, mask = mask) - nexta = sum (hillslope_basin(i)%area, mask = mask) + hmask = hillslope_basin(i)%hand < wdsrf_bsn(i) + handmin(i) + nextl = maxval(hillslope_basin(i)%hand, mask = hmask) + nexta = sum (hillslope_basin(i)%area, mask = hmask) nextv = nexta * (wdsrf_bsn(i)+handmin(i)-nextl) IF (nextv > dvol) THEN ddep = dvol/nexta @@ -362,22 +412,22 @@ SUBROUTINE river_lake_flow (dt) wdsrf_bsn(i) = wdsrf_bsn(i) - ddep DO j = 1, hillslope_basin(i)%nhru - IF (mask(j)) THEN + IF (hmask(j)) THEN wdsrf_bsnhru(j+hs-1) = wdsrf_bsnhru(j+hs-1) - ddep ENDIF ENDDO ENDDO ELSEIF (dvol < -VOLUMEMIN) THEN - mask = .true. + hmask = .true. nexta = 0. DO WHILE (dvol < -VOLUMEMIN) - IF (any(mask)) THEN - j = minloc(hillslope_basin(i)%hand + wdsrf_bsnhru(hs:he), 1, mask = mask) + IF (any(hmask)) THEN + j = minloc(hillslope_basin(i)%hand + wdsrf_bsnhru(hs:he), 1, mask = hmask) nexta = nexta + hillslope_basin(i)%area(j) - mask(j) = .false. + hmask(j) = .false. ENDIF - IF (any(mask)) THEN - nextl = minval(hillslope_basin(i)%hand + wdsrf_bsnhru(hs:he), mask = mask) + IF (any(hmask)) THEN + nextl = minval(hillslope_basin(i)%hand + wdsrf_bsnhru(hs:he), mask = hmask) nextv = nexta*(nextl-(wdsrf_bsn(i)+handmin(i))) IF ((-dvol) > nextv) THEN ddep = nextl - (wdsrf_bsn(i)+handmin(i)) @@ -394,7 +444,7 @@ SUBROUTINE river_lake_flow (dt) wdsrf_bsn(i) = wdsrf_bsn(i) + ddep DO j = 1, hillslope_basin(i)%nhru - IF (.not. mask(j)) THEN + IF (.not. hmask(j)) THEN wdsrf_bsnhru(j+hs-1) = wdsrf_bsnhru(j+hs-1) + ddep ENDIF ENDDO @@ -402,11 +452,11 @@ SUBROUTINE river_lake_flow (dt) ENDIF ENDIF - deallocate(mask) + deallocate(hmask) ELSE totalvolume = lakeinfo(i)%volume(wdsrf_bsn(i)) - totalvolume = totalvolume - sum_hflux_riv(i) * dt_this + totalvolume = totalvolume - sum_hflux_riv(i) * dt_all(rivermouth(i)) wdsrf_bsn(i) = lakeinfo(i)%surface(totalvolume) ENDIF @@ -416,8 +466,8 @@ SUBROUTINE river_lake_flow (dt) ELSE friction = grav * nmanning_riv**2 / wdsrf_bsn(i)**(7.0/3.0) * abs(momen_riv(i)) momen_riv(i) = (momen_riv(i) & - - (sum_mflux_riv(i) - sum_zgrad_riv(i)) / riverarea(i) * dt_this) & - / (1 + friction * dt_this) + - (sum_mflux_riv(i) - sum_zgrad_riv(i)) / riverarea(i) * dt_all(rivermouth(i))) & + / (1 + friction * dt_all(rivermouth(i))) veloc_riv(i) = momen_riv(i) / wdsrf_bsn(i) ENDIF @@ -432,29 +482,37 @@ SUBROUTINE river_lake_flow (dt) ENDDO - IF (numbasin > 0) THEN - wdsrf_bsn_ta (:) = wdsrf_bsn_ta (:) + wdsrf_bsn(:) * dt_this - momen_riv_ta (:) = momen_riv_ta (:) + momen_riv(:) * dt_this - discharge_ta (:) = discharge_ta (:) + hflux_fc (:) * dt_this - ENDIF - DO i = 1, numbasin - IF (lake_id(i) > 0) THEN ! for lakes - hs = basin_hru%substt(i) - he = basin_hru%subend(i) - DO j = hs, he - wdsrf_bsnhru(j) = max(wdsrf_bsn(i) - (lakeinfo(i)%depth(1) - lakeinfo(i)%depth0(j-hs+1)), 0.) - wdsrf_bsnhru_ta(j) = wdsrf_bsnhru_ta(j) + wdsrf_bsnhru(j) * dt_this - ENDDO + IF (bsnfilter(i)) THEN + + wdsrf_bsn_ta (i) = wdsrf_bsn_ta (i) + wdsrf_bsn(i) * dt_all(rivermouth(i)) + momen_riv_ta (i) = momen_riv_ta (i) + momen_riv(i) * dt_all(rivermouth(i)) + discharge_ta (i) = discharge_ta (i) + hflux_fc (i) * dt_all(rivermouth(i)) + + IF (lake_id(i) > 0) THEN ! for lakes + hs = basin_hru%substt(i) + he = basin_hru%subend(i) + DO j = hs, he + wdsrf_bsnhru(j) = max(wdsrf_bsn(i) - (lakeinfo(i)%depth(1) - lakeinfo(i)%depth0(j-hs+1)), 0.) + wdsrf_bsnhru_ta(j) = wdsrf_bsnhru_ta(j) + wdsrf_bsnhru(j) * dt_all(rivermouth(i)) + ENDDO + ENDIF + ENDIF ENDDO - dt_res = dt_res - dt_this - + dt_res = dt_res - dt_all + + allfinished = all(dt_res == 0) +#ifdef USEMPI + CALL mpi_allreduce (MPI_IN_PLACE, allfinished, 1, MPI_LOGICAL, MPI_LAND, p_comm_worker, p_err) +#endif ENDDO IF (numbasin > 0) wdsrf_bsn_prev(:) = wdsrf_bsn(:) + CALL mpi_allreduce (MPI_IN_PLACE, ntimestep_riverlake, 1, MPI_INTEGER, MPI_MAX, p_comm_worker, p_err) + IF (allocated(wdsrf_bsn_ds )) deallocate(wdsrf_bsn_ds ) IF (allocated(veloc_riv_ds )) deallocate(veloc_riv_ds ) IF (allocated(momen_riv_ds )) deallocate(momen_riv_ds ) @@ -464,6 +522,8 @@ SUBROUTINE river_lake_flow (dt) IF (allocated(sum_hflux_riv)) deallocate(sum_hflux_riv) IF (allocated(sum_mflux_riv)) deallocate(sum_mflux_riv) IF (allocated(sum_zgrad_riv)) deallocate(sum_zgrad_riv) + IF (allocated(bsnfilter )) deallocate(bsnfilter ) + IF (allocated(syncfilter )) deallocate(syncfilter ) ENDIF diff --git a/main/HYDRO/MOD_Catch_RiverLakeNetwork.F90 b/main/HYDRO/MOD_Catch_RiverLakeNetwork.F90 index 07ea5470..e23ab0ca 100644 --- a/main/HYDRO/MOD_Catch_RiverLakeNetwork.F90 +++ b/main/HYDRO/MOD_Catch_RiverLakeNetwork.F90 @@ -34,6 +34,7 @@ MODULE MOD_Catch_RiverLakeNetwork ! index of downstream river ! > 0 : other catchment; 0 : river mouth; -1 : inland depression integer, allocatable :: riverdown (:) + integer, allocatable :: ilocdown (:) logical, allocatable :: to_lake (:) real(r8), allocatable :: riverlen_ds (:) @@ -43,6 +44,18 @@ MODULE MOD_Catch_RiverLakeNetwork real(r8), allocatable :: outletwth (:) + integer :: numbsnlink +#ifdef USEMPI + integer, allocatable :: linkbindex (:) + integer, allocatable :: linkrivmth (:) + + integer :: nlink_me + integer, allocatable :: linkpush (:) + integer, allocatable :: linkpull (:) + + logical, allocatable :: riv_multi_proc (:) +#endif + ! -- lake data type -- type :: lake_info_type integer :: nsub @@ -66,8 +79,14 @@ MODULE MOD_Catch_RiverLakeNetwork ! -- information of HRU in basin -- type(hillslope_network_type), pointer :: hillslope_basin (:) - type(basin_pushdata_type), target :: river_iam_dn - type(basin_pushdata_type), target :: river_iam_up + + ! ----- subroutines and functions ----- + PUBLIC :: river_lake_network_init + PUBLIC :: pull_from_downstream + PUBLIC :: push_to_downstream + PUBLIC :: calc_riverdepth_from_runoff + PUBLIC :: retrieve_lake_surface_from_volume + PUBLIC :: river_lake_network_final CONTAINS @@ -93,23 +112,20 @@ SUBROUTINE river_lake_network_init () character(len=256) :: river_file, rivdpt_file logical :: use_calc_rivdpt - integer :: totalnumbasin, ibasin, inb - integer :: iworker, mesg(4), isrc, idest, iproc - integer :: nrecv, irecv, ifrom, ito, iup, idn - integer :: ndata, idata, ip, nup, ndn - integer :: iloc, iloc1, iloc2 - integer :: ielm, i, j, ithis, nave + integer :: totalnumbasin, ibasin, inb, iloc, ielm, i, j, ilink + integer :: iworker, mesg(2), isrc, idest, nrecv + logical , allocatable :: is_link (:) + logical , allocatable :: link_on_me (:) integer , allocatable :: icache (:) real(r8), allocatable :: rcache (:) logical , allocatable :: lcache (:) - type(pointer_int32_2d), allocatable :: datapush_w (:) - integer, allocatable :: datapush(:,:) - integer, allocatable :: bindex(:), addrbasin(:), addrdown(:), nelm_wrk(:), paddr(:), ndata_w(:) + integer, allocatable :: bindex(:), addrbasin(:), addrdown(:) integer, allocatable :: basin_sorted(:), basin_order(:), order (:) - integer, allocatable :: river_up_ups(:), river_up_paddr(:), river_dn_ups(:), river_dn_paddr(:) + + logical, allocatable :: bsnfilter(:), linkfilter(:) ! for lakes integer :: ps, pe, nsublake, hs, he, ihru, ipxl @@ -301,75 +317,51 @@ SUBROUTINE river_lake_network_init () ENDIF - CALL mpi_barrier (p_comm_glb, p_err) -#endif -#ifdef USEMPI IF (p_is_master) THEN - allocate (addrdown (totalnumbasin)) - allocate (ndata_w (0:p_np_worker-1)) + allocate (addrdown (totalnumbasin)) + allocate (is_link (totalnumbasin)); is_link = .false. - ndata_w(:) = 0 DO ibasin = 1, totalnumbasin IF (riverdown(ibasin) >= 1) THEN addrdown(ibasin) = addrbasin(riverdown(ibasin)) - IF (addrbasin(ibasin) /= addrdown(ibasin)) THEN - ifrom = p_itis_worker(addrbasin(ibasin)) - ito = p_itis_worker(addrdown(ibasin)) - ndata_w(ifrom) = ndata_w(ifrom) + 1 - ndata_w(ito) = ndata_w(ito) + 1 + IF (addrdown(ibasin) /= addrbasin(ibasin)) THEN + is_link(riverdown(ibasin)) = .true. ENDIF + ELSE + addrdown(ibasin) = addrbasin(ibasin) ENDIF ENDDO + + numbsnlink = count(is_link) - allocate (datapush_w (0:p_np_worker-1)) - DO iworker = 0, p_np_worker-1 - IF (ndata_w(iworker) > 0) THEN - allocate (datapush_w(iworker)%val (4,ndata_w(iworker))) - ENDIF - ENDDO + IF (numbsnlink > 0) THEN + allocate (linkbindex (numbsnlink)) + allocate (linkrivmth (numbsnlink)) + linkbindex = pack((/(ibasin, ibasin = 1, totalnumbasin)/), mask = is_link) + linkrivmth = pack(rivermouth, mask = is_link) + ENDIF - ndata_w(:) = 0 - DO ibasin = 1, totalnumbasin - IF ((riverdown(ibasin) >= 1) .and. (addrbasin(ibasin) /= addrdown(ibasin))) THEN - ifrom = p_itis_worker(addrbasin(ibasin)) - ito = p_itis_worker(addrdown(ibasin)) - ndata_w(ifrom) = ndata_w(ifrom) + 1 - ndata_w(ito) = ndata_w(ito) + 1 - - datapush_w(ifrom)%val(:,ndata_w(ifrom)) = & - (/addrbasin(ibasin), ibasin, addrdown(ibasin), riverdown(ibasin)/) - datapush_w(ito)%val(:,ndata_w(ito)) = & - (/addrbasin(ibasin), ibasin, addrdown(ibasin), riverdown(ibasin)/) - ENDIF - ENDDO + deallocate (addrdown) + deallocate (is_link ) - DO iworker = 0, p_np_worker-1 - CALL mpi_send (ndata_w(iworker), 1, MPI_INTEGER, & - p_address_worker(iworker), mpi_tag_size, p_comm_glb, p_err) - IF (ndata_w(iworker) > 0) THEN - CALL mpi_send (datapush_w(iworker)%val, 4*ndata_w(iworker), MPI_INTEGER, & - p_address_worker(iworker), mpi_tag_data, p_comm_glb, p_err) - ENDIF - ENDDO + write(*,'(/,A,I5,A,/)') 'There are ', numbsnlink, ' links between processors.' - deallocate (addrdown ) - deallocate (ndata_w ) - deallocate (datapush_w) + ENDIF + CALL mpi_bcast (numbsnlink, 1, MPI_INTEGER, p_address_master, p_comm_glb, p_err) + IF (numbsnlink > 0) THEN + IF (.not. allocated(linkbindex)) allocate (linkbindex (numbsnlink)) + IF (.not. allocated(linkrivmth)) allocate (linkrivmth (numbsnlink)) + CALL mpi_bcast (linkbindex, numbsnlink, MPI_INTEGER, p_address_master, p_comm_glb, p_err) + CALL mpi_bcast (linkrivmth, numbsnlink, MPI_INTEGER, p_address_master, p_comm_glb, p_err) ENDIF +#else + numbsnlink = 0 #endif IF (p_is_worker) THEN -#ifdef USEMPI - CALL mpi_recv (ndata, 1, MPI_INTEGER, p_address_master, mpi_tag_size, p_comm_glb, p_stat, p_err) - IF (ndata > 0) THEN - allocate (datapush(4,ndata)) - CALL mpi_recv (datapush, 4*ndata, MPI_INTEGER, & - p_address_master, mpi_tag_data, p_comm_glb, p_stat, p_err) - ENDIF -#endif IF (numbasin > 0) THEN allocate (basin_sorted (numbasin)) @@ -380,164 +372,61 @@ SUBROUTINE river_lake_network_init () CALL quicksort (numbasin, basin_sorted, basin_order) ENDIF - river_iam_up%nself = 0 - river_iam_dn%nself = 0 - river_iam_up%nproc = 0 - river_iam_dn%nproc = 0 - IF (numbasin > 0) THEN + + allocate (ilocdown (numbasin)); ilocdown(:) = 0 DO ibasin = 1, numbasin IF (riverdown(ibasin) > 0) THEN iloc = find_in_sorted_list1 (riverdown(ibasin), numbasin, basin_sorted) IF (iloc > 0) THEN - river_iam_up%nself = river_iam_up%nself + 1 + ilocdown(ibasin) = basin_order(iloc) +#ifdef USEMPI + ELSE + iloc = find_in_sorted_list1 (riverdown(ibasin), numbsnlink, linkbindex) + ilocdown(ibasin) = - iloc +#endif ENDIF ENDIF ENDDO - - IF (river_iam_up%nself > 0) THEN - - allocate (river_iam_up%iself (river_iam_up%nself)) - - river_iam_dn%nself = river_iam_up%nself - allocate (river_iam_dn%iself (river_iam_dn%nself)) - - idata = 0 - DO ibasin = 1, numbasin - IF (riverdown(ibasin) > 0) THEN - iloc = find_in_sorted_list1 (riverdown(ibasin), numbasin, basin_sorted) - IF (iloc > 0) THEN - idata = idata + 1 - river_iam_up%iself(idata) = ibasin - river_iam_dn%iself(idata) = basin_order(iloc) - ENDIF - ENDIF - ENDDO - - ENDIF + ENDIF #ifdef USEMPI - IF (ndata > 0) THEN - - nup = count(datapush(3,:) == p_iam_glb) - ndn = count(datapush(1,:) == p_iam_glb) - - IF (nup > 0) allocate (river_up_paddr (nup)) - IF (nup > 0) allocate (river_up_ups (nup)) - IF (ndn > 0) allocate (river_dn_paddr (ndn)) - IF (ndn > 0) allocate (river_dn_ups (ndn)) - - iup = 0 - idn = 0 - DO idata = 1, ndata - IF (datapush(3,idata) == p_iam_glb) THEN - CALL insert_into_sorted_list2 (datapush(2,idata), datapush(1,idata), & - iup, river_up_ups, river_up_paddr, iloc) - ELSEIF (datapush(1,idata) == p_iam_glb) THEN - CALL insert_into_sorted_list2 (datapush(2,idata), datapush(3,idata), & - idn, river_dn_ups, river_dn_paddr, iloc) - ENDIF - ENDDO - - IF (nup > 0) allocate (river_iam_dn%ipush (nup)) - IF (ndn > 0) allocate (river_iam_up%ipush (ndn)) - - DO idata = 1, ndata - IF (datapush(3,idata) == p_iam_glb) THEN - - iloc1 = find_in_sorted_list2 (datapush(2,idata), datapush(1,idata), & - nup, river_up_ups, river_up_paddr) - iloc2 = find_in_sorted_list1 (datapush(4,idata), numbasin, basin_sorted) - - river_iam_dn%ipush(iloc1) = basin_order(iloc2) - - ELSEIF (datapush(1,idata) == p_iam_glb) THEN - - iloc1 = find_in_sorted_list2 (datapush(2,idata), datapush(3,idata), & - ndn, river_dn_ups, river_dn_paddr) - iloc2 = find_in_sorted_list1 (datapush(2,idata), numbasin, basin_sorted) - - river_iam_up%ipush(iloc1) = basin_order(iloc2) - - ENDIF - ENDDO - + IF ((numbsnlink > 0) .and. (numbasin > 0)) THEN - IF (nup > 0) THEN - - DO iup = 1, nup - IF (iup == 1) THEN - river_iam_dn%nproc = 1 - ELSEIF (river_up_paddr(iup) /= river_up_paddr(iup-1)) THEN - river_iam_dn%nproc = river_iam_dn%nproc + 1 - ENDIF - ENDDO - - allocate (river_iam_dn%paddr(river_iam_dn%nproc)) - allocate (river_iam_dn%ndata(river_iam_dn%nproc)) - - DO iup = 1, nup - IF (iup == 1) THEN - ip = 1 - river_iam_dn%paddr(ip) = river_up_paddr(iup) - river_iam_dn%ndata(ip) = 1 - ELSEIF (river_up_paddr(iup) /= river_up_paddr(iup-1)) THEN - ip = ip + 1 - river_iam_dn%paddr(ip) = river_up_paddr(iup) - river_iam_dn%ndata(ip) = 1 - ELSE - river_iam_dn%ndata(ip) = river_iam_dn%ndata(ip) + 1 - ENDIF - ENDDO + allocate (link_on_me (numbsnlink)); link_on_me = .false. + DO ibasin = 1, numbsnlink + iloc = find_in_sorted_list1 (linkbindex(ibasin), numbasin, basin_sorted) + IF (iloc > 0) THEN + link_on_me(ibasin) = .true. ENDIF - + ENDDO - IF (ndn > 0) THEN + nlink_me = count(link_on_me) + + IF (nlink_me > 0) THEN + allocate (linkpush (nlink_me)) + allocate (linkpull (nlink_me)) + ilink = 0 + DO ibasin = 1, numbsnlink + IF (link_on_me(ibasin)) THEN + ilink = ilink + 1 + iloc = find_in_sorted_list1 (linkbindex(ibasin), numbasin, basin_sorted) + linkpush(ilink) = basin_order(iloc) + linkpull(ilink) = ibasin + ENDIF + ENDDO + ENDIF - DO idn = 1, ndn - IF (idn == 1) THEN - river_iam_up%nproc = 1 - ELSEIF (river_dn_paddr(idn) /= river_dn_paddr(idn-1)) THEN - river_iam_up%nproc = river_iam_up%nproc + 1 - ENDIF - ENDDO - - allocate (river_iam_up%paddr(river_iam_up%nproc)) - allocate (river_iam_up%ndata(river_iam_up%nproc)) - - DO idn = 1, ndn - IF (idn == 1) THEN - ip = 1 - river_iam_up%paddr(ip) = river_dn_paddr(idn) - river_iam_up%ndata(ip) = 1 - ELSEIF (river_dn_paddr(idn) /= river_dn_paddr(idn-1)) THEN - ip = ip + 1 - river_iam_up%paddr(ip) = river_dn_paddr(idn) - river_iam_up%ndata(ip) = 1 - ELSE - river_iam_up%ndata(ip) = river_iam_up%ndata(ip) + 1 - ENDIF - ENDDO + deallocate (link_on_me) - ENDIF - - IF (nup > 0) deallocate (river_up_paddr) - IF (nup > 0) deallocate (river_up_ups ) - IF (ndn > 0) deallocate (river_dn_paddr) - IF (ndn > 0) deallocate (river_dn_ups ) - - deallocate(datapush) - - ENDIF -#endif + ELSE + nlink_me = 0 ENDIF +#endif ENDIF - IF (allocated(basin_sorted )) deallocate(basin_sorted ) - IF (allocated(basin_order )) deallocate(basin_order ) - CALL hillslope_network_init (numbasin, basinindex, hillslope_basin) @@ -569,8 +458,8 @@ SUBROUTINE river_lake_network_init () ENDIF - CALL worker_push_data (iam_bsn, iam_elm, .false., lake_id, lake_id_elm) - CALL worker_push_data (iam_bsn, iam_elm, .false., lakedown_id_bsn, lakedown_id_elm) + CALL worker_push_data (iam_bsn, iam_elm, lake_id, lake_id_elm) + CALL worker_push_data (iam_bsn, iam_elm, lakedown_id_bsn, lakedown_id_elm) IF (p_is_worker) THEN @@ -605,7 +494,7 @@ SUBROUTINE river_lake_network_init () ENDDO ENDIF - CALL worker_push_data (iam_elm, iam_bsn, .false., lakeoutlet_elm, lakeoutlet_bsn) + CALL worker_push_data (iam_elm, iam_bsn, lakeoutlet_elm, lakeoutlet_bsn) CALL worker_push_subset_data (iam_elm, iam_bsn, elm_hru, basin_hru, lakedepth_hru, lakedepth_bsnhru) CALL worker_push_subset_data (iam_elm, iam_bsn, elm_hru, basin_hru, lakearea_hru, lakearea_bsnhru ) @@ -656,7 +545,7 @@ SUBROUTINE river_lake_network_init () wtsrfelv(ibasin) = basinelv(ibasin) - bedelv(ibasin) = basinelv(ibasin) - minval(lakedepth_bsnhru(hs:he)) + bedelv(ibasin) = basinelv(ibasin) - maxval(lakedepth_bsnhru(hs:he)) nsublake = he - hs + 1 lakeinfo(ibasin)%nsub = nsublake @@ -709,10 +598,16 @@ SUBROUTINE river_lake_network_init () ENDDO ENDIF - CALL worker_push_data (river_iam_dn, river_iam_up, .false., riverlen, riverlen_ds) - CALL worker_push_data (river_iam_dn, river_iam_up, .false., wtsrfelv, wtsrfelv_ds) - CALL worker_push_data (river_iam_dn, river_iam_up, .false., riverwth, riverwth_ds) - CALL worker_push_data (river_iam_dn, river_iam_up, .false., bedelv , bedelv_ds ) + IF (numbasin > 0) allocate(bsnfilter (numbasin )); bsnfilter(:) = .true. + IF (numbsnlink > 0) allocate(linkfilter(numbsnlink)); linkfilter(:) = .true. + + CALL pull_from_downstream (riverlen, riverlen_ds, bsnfilter, linkfilter) + CALL pull_from_downstream (wtsrfelv, wtsrfelv_ds, bsnfilter, linkfilter) + CALL pull_from_downstream (riverwth, riverwth_ds, bsnfilter, linkfilter) + CALL pull_from_downstream (bedelv , bedelv_ds , bsnfilter, linkfilter) + + IF (allocated(bsnfilter )) deallocate(bsnfilter ) + IF (allocated(linkfilter)) deallocate(linkfilter) DO ibasin = 1, numbasin IF (lake_id(ibasin) < 0) THEN @@ -754,6 +649,31 @@ SUBROUTINE river_lake_network_init () IF (allocated (lakeoutlet_bsn )) deallocate (lakeoutlet_bsn ) IF (allocated (lakearea_bsnhru )) deallocate (lakearea_bsnhru ) + + IF (p_is_worker) THEN + + allocate (riv_multi_proc (numrivmth)) + + riv_multi_proc(:) = .false. + + DO ibasin = 1, numbasin + j = riverdown(ibasin) + IF (j > 0) THEN + iloc = find_in_sorted_list1 (j, numbasin, basin_sorted) + IF (iloc <= 0) THEN + riv_multi_proc(rivermouth(ibasin)) = .true. + ENDIF + ENDIF + ENDDO + + CALL mpi_allreduce (MPI_IN_PLACE, riv_multi_proc, numrivmth, & + MPI_LOGICAL, MPI_LOR, p_comm_worker, p_err) + + ENDIF + + IF (allocated(basin_sorted )) deallocate(basin_sorted ) + IF (allocated(basin_order )) deallocate(basin_order ) + #ifdef USEMPI CALL mpi_barrier (p_comm_glb, p_err) IF (p_is_master) write(*,'(A)') 'Building river network information done.' @@ -762,6 +682,146 @@ SUBROUTINE river_lake_network_init () END SUBROUTINE river_lake_network_init + ! ----- pull data from downstream basin ----- + SUBROUTINE pull_from_downstream (datain, dataout, bsnfilter, linkfilter) + + USE MOD_SPMD_Task + + IMPLICIT NONE + + real(r8), intent(in) :: datain (:) + real(r8), intent(inout) :: dataout(:) + + logical, intent(in) :: bsnfilter (:) + logical, intent(in) :: linkfilter (:) + + ! local variables + integer :: ndata, i, idata, ibasin + real(r8), allocatable :: xdata (:) + real(r8), allocatable :: datalink(:), fill(:) + + IF (p_is_worker) THEN + +#ifdef USEMPI + IF (numbsnlink > 0) THEN + IF (any(linkfilter)) THEN + + ndata = count(linkfilter) + allocate (xdata (ndata)); xdata(:) = 0. + DO i = 1, nlink_me + IF (linkfilter(linkpull(i))) THEN + idata = count(linkfilter(1:linkpull(i))) + xdata(idata) = datain(linkpush(i)) + ENDIF + ENDDO + + CALL mpi_allreduce (MPI_IN_PLACE, xdata, ndata, MPI_REAL8, MPI_SUM, p_comm_worker, p_err) + + allocate (datalink (numbsnlink)) + allocate (fill (numbsnlink)); fill(:) = 0. + datalink = unpack(xdata, linkfilter, fill) + + deallocate (xdata) + deallocate (fill ) + + ENDIF + ENDIF +#endif + + DO ibasin = 1, numbasin + IF (bsnfilter(ibasin)) THEN + IF (riverdown(ibasin) > 0) THEN + IF (ilocdown(ibasin) > 0) THEN + dataout(ibasin) = datain(ilocdown(ibasin)) +#ifdef USEMPI + ELSEIF (ilocdown(ibasin) < 0) THEN + dataout(ibasin) = datalink(-ilocdown(ibasin)) +#endif + ENDIF + ENDIF + ENDIF + ENDDO + +#ifdef USEMPI + if (allocated(datalink)) deallocate (datalink) +#endif + ENDIF + + END SUBROUTINE pull_from_downstream + + ! ----- push data to downstream basin ----- + SUBROUTINE push_to_downstream (datain, dataout, bsnfilter, linkfilter) + + USE MOD_SPMD_Task + + IMPLICIT NONE + + real(r8), intent(in) :: datain (:) + real(r8), intent(inout) :: dataout(:) + + logical, intent(in) :: bsnfilter (:) + logical, intent(in) :: linkfilter (:) + + ! local variables + integer :: ndata, i, idata, ibasin + real(r8), allocatable :: xdata (:) + real(r8), allocatable :: datalink(:), fill(:) + + + IF (p_is_worker) THEN + +#ifdef USEMPI + IF (numbsnlink > 0) THEN + allocate (datalink (numbsnlink)) + datalink(:) = 0. + ENDIF +#endif + + DO ibasin = 1, numbasin + IF (bsnfilter(ibasin)) THEN + IF (riverdown(ibasin) > 0) THEN + IF (ilocdown(ibasin) > 0) THEN + dataout(ilocdown(ibasin)) = dataout(ilocdown(ibasin)) + datain(ibasin) +#ifdef USEMPI + ELSEIF (ilocdown(ibasin) < 0) THEN + datalink(-ilocdown(ibasin)) = datalink(-ilocdown(ibasin)) + datain(ibasin) +#endif + ENDIF + ENDIF + ENDIF + ENDDO + +#ifdef USEMPI + IF (numbsnlink > 0) THEN + IF (any(linkfilter)) THEN + + ndata = count(linkfilter) + allocate (xdata (ndata)) + xdata(:) = pack(datalink, linkfilter) + + CALL mpi_allreduce (MPI_IN_PLACE, xdata, ndata, MPI_REAL8, MPI_SUM, p_comm_worker, p_err) + + allocate (fill (numbsnlink)); fill(:) = 0. + datalink = unpack(xdata, linkfilter, fill) + + DO i = 1, nlink_me + IF (linkfilter(linkpull(i))) THEN + idata = count(linkfilter(1:linkpull(i))) + dataout(linkpush(i)) = dataout(linkpush(i)) + xdata(idata) + ENDIF + ENDDO + + deallocate (datalink) + deallocate (xdata ) + deallocate (fill ) + + ENDIF + ENDIF +#endif + ENDIF + + END SUBROUTINE push_to_downstream + ! ----- retrieve river depth from runoff ----- SUBROUTINE calc_riverdepth_from_runoff () @@ -1024,28 +1084,33 @@ SUBROUTINE river_lake_network_final () ! Local Variables integer :: ilake - IF (allocated(lake_id )) deallocate(lake_id ) - IF (allocated(riverlen )) deallocate(riverlen ) - IF (allocated(riverelv )) deallocate(riverelv ) - IF (allocated(riverarea)) deallocate(riverarea) - IF (allocated(riverwth )) deallocate(riverwth ) - IF (allocated(riverdpth)) deallocate(riverdpth) - IF (allocated(basinelv )) deallocate(basinelv ) - IF (allocated(bedelv )) deallocate(bedelv ) - IF (allocated(handmin )) deallocate(handmin ) - IF (allocated(wtsrfelv )) deallocate(wtsrfelv ) - IF (allocated(riverdown)) deallocate(riverdown) - IF (allocated(to_lake )) deallocate(to_lake ) - - IF (allocated(riverlen_ds)) deallocate(riverlen_ds) - IF (allocated(wtsrfelv_ds)) deallocate(wtsrfelv_ds) - IF (allocated(riverwth_ds)) deallocate(riverwth_ds) - IF (allocated(bedelv_ds )) deallocate(bedelv_ds ) - IF (allocated(outletwth )) deallocate(outletwth ) - - IF (allocated(lakeinfo)) deallocate(lakeinfo) - - IF (associated(hillslope_basin)) deallocate(hillslope_basin) + IF (allocated (riverlen )) deallocate(riverlen ) + IF (allocated (riverelv )) deallocate(riverelv ) + IF (allocated (riverwth )) deallocate(riverwth ) + IF (allocated (riverarea )) deallocate(riverarea ) + IF (allocated (riverdpth )) deallocate(riverdpth ) + IF (allocated (basinelv )) deallocate(basinelv ) + IF (allocated (bedelv )) deallocate(bedelv ) + IF (allocated (handmin )) deallocate(handmin ) + IF (allocated (wtsrfelv )) deallocate(wtsrfelv ) + IF (allocated (riverdown )) deallocate(riverdown ) + IF (allocated (ilocdown )) deallocate(ilocdown ) + IF (allocated (to_lake )) deallocate(to_lake ) + IF (allocated (riverlen_ds )) deallocate(riverlen_ds ) + IF (allocated (wtsrfelv_ds )) deallocate(wtsrfelv_ds ) + IF (allocated (riverwth_ds )) deallocate(riverwth_ds ) + IF (allocated (bedelv_ds )) deallocate(bedelv_ds ) + IF (allocated (outletwth )) deallocate(outletwth ) +#ifdef USEMPI + IF (allocated (linkbindex )) deallocate(linkbindex ) + IF (allocated (linkrivmth )) deallocate(linkrivmth ) + IF (allocated (linkpush )) deallocate(linkpush ) + IF (allocated (linkpull )) deallocate(linkpull ) + IF (allocated (riv_multi_proc )) deallocate(riv_multi_proc ) +#endif + IF (allocated (lake_id )) deallocate(lake_id ) + IF (allocated (lakeinfo )) deallocate(lakeinfo ) + IF (associated(hillslope_basin)) deallocate(hillslope_basin) END SUBROUTINE river_lake_network_final diff --git a/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 b/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 index 1d20905a..613017e6 100644 --- a/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 +++ b/main/HYDRO/MOD_Catch_SubsurfaceFlow.F90 @@ -94,8 +94,8 @@ SUBROUTINE subsurface_network_init () IF (numelm > 0) allocate (lakedepth_elm(numelm)) IF (numelm > 0) allocate (wdsrf_elm (numelm)) - CALL worker_push_data (iam_bsn, iam_elm, .false., lake_id, lake_id_elm ) - CALL worker_push_data (iam_bsn, iam_elm, .false., riverdpth, riverdpth_elm) + CALL worker_push_data (iam_bsn, iam_elm, lake_id, lake_id_elm ) + CALL worker_push_data (iam_bsn, iam_elm, riverdpth, riverdpth_elm) DO ielm = 1, numelm IF (lake_id_elm(ielm) <= 0) THEN diff --git a/main/HYDRO/MOD_Catch_Vars_TimeVariables.F90 b/main/HYDRO/MOD_Catch_Vars_TimeVariables.F90 index 2a7e20f1..0d0943f5 100644 --- a/main/HYDRO/MOD_Catch_Vars_TimeVariables.F90 +++ b/main/HYDRO/MOD_Catch_Vars_TimeVariables.F90 @@ -88,10 +88,10 @@ SUBROUTINE READ_CatchTimeVariables (file_restart) character(len=*), intent(in) :: file_restart CALL vector_read_basin (file_restart, veloc_elm, numelm, 'veloc_riv', elm_data_address) - CALL worker_push_data (iam_elm, iam_bsn, .false., veloc_elm, veloc_riv) + CALL worker_push_data (iam_elm, iam_bsn, veloc_elm, veloc_riv) CALL vector_read_basin (file_restart, wdsrf_elm_prev, numelm, 'wdsrf_bsn_prev', elm_data_address) - CALL worker_push_data (iam_elm, iam_bsn, .false., wdsrf_elm_prev, wdsrf_bsn_prev) + CALL worker_push_data (iam_elm, iam_bsn, wdsrf_elm_prev, wdsrf_bsn_prev) CALL vector_read_basin (file_restart, veloc_hru, numhru, 'veloc_hru', hru_data_address) CALL worker_push_subset_data (iam_elm, iam_bsn, elm_hru, basin_hru, veloc_hru, veloc_bsnhru) @@ -132,11 +132,11 @@ SUBROUTINE WRITE_CatchTimeVariables (file_restart) 'long_name', 'index of hydrological units inside basin') ENDIF - CALL worker_push_data (iam_bsn, iam_elm, .false., veloc_riv, veloc_elm) + CALL worker_push_data (iam_bsn, iam_elm, veloc_riv, veloc_elm) CALL vector_write_basin (& file_restart, veloc_elm, numelm, totalnumelm, 'veloc_riv', 'basin', elm_data_address) - CALL worker_push_data (iam_bsn, iam_elm, .false., wdsrf_bsn_prev, wdsrf_elm_prev) + CALL worker_push_data (iam_bsn, iam_elm, wdsrf_bsn_prev, wdsrf_elm_prev) CALL vector_write_basin (& file_restart, wdsrf_elm_prev, numelm, totalnumelm, 'wdsrf_bsn_prev', 'basin', elm_data_address) diff --git a/main/MOD_LAIReadin.F90 b/main/MOD_LAIReadin.F90 index 7e394966..77cb80dd 100644 --- a/main/MOD_LAIReadin.F90 +++ b/main/MOD_LAIReadin.F90 @@ -197,6 +197,13 @@ SUBROUTINE LAI_readin (year, time, dir_landdata) green(npatch) = 1. fveg (npatch) = fveg0(m) + IF (m == WATERBODY) THEN + fveg(npatch) = 0. + tlai(npatch) = 0. + tsai(npatch) = 0. + green(npatch) = 0. + ENDIF + ENDDO ENDIF ENDIF diff --git a/main/MOD_Lake.F90 b/main/MOD_Lake.F90 index b589825d..c1a7142a 100644 --- a/main/MOD_Lake.F90 +++ b/main/MOD_Lake.F90 @@ -2103,6 +2103,8 @@ SUBROUTINE adjust_lake_layer (nl_lake, dz_lake, t_lake, lake_icefrac) dz_lake_new(nl_lake) = dzlak(nl_lake)*depthratio - (dz_lake_new(1) - dzlak(1)*depthratio) ELSEIF(wdsrfm > 0. .and. wdsrfm <= 1.)THEN dz_lake_new(:) = wdsrfm / nl_lake + ELSE + write(*,*) 'Warning: lake depth is over 2000 meters!' ENDIF j = 1 diff --git a/mkinidata/MOD_Initialize.F90 b/mkinidata/MOD_Initialize.F90 index 48c69ddb..b1b4f4ee 100644 --- a/mkinidata/MOD_Initialize.F90 +++ b/mkinidata/MOD_Initialize.F90 @@ -1407,7 +1407,7 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & hs = basin_hru%substt(i) he = basin_hru%subend(i) IF (lake_id(i) <= 0) THEN - wdsrf_bsn(i) = minval(hillslope_basin(i)%hand + wdsrf_bsnhru(hs:he)) + wdsrf_bsn(i) = minval(hillslope_basin(i)%hand + wdsrf_bsnhru(hs:he)) - handmin(i) ELSE ! lake totalvolume = sum(wdsrf_bsnhru(hs:he) * lakeinfo(i)%area0) @@ -1442,6 +1442,9 @@ SUBROUTINE initialize (casename, dir_landdata, dir_restart, & dz_lake(:,i) = 0. ENDIF ENDDO + + CALL check_vector_data ('Basin Water Depth [m] ', wdsrf_bsn) + CALL check_vector_data ('HRU Water Depth [m] ', wdsrf_bsnhru) ENDIF #endif diff --git a/mksrfdata/Aggregation_LAI.F90 b/mksrfdata/Aggregation_LAI.F90 index 3ae2cd2f..d448d2e7 100644 --- a/mksrfdata/Aggregation_LAI.F90 +++ b/mksrfdata/Aggregation_LAI.F90 @@ -428,6 +428,10 @@ SUBROUTINE Aggregation_LAI (gridlai, dir_rawdata, dir_model_landdata, lc_year) write(cyear,'(i4.4)') iy suffix = 'MOD'//trim(cyear) CALL system('mkdir -p ' // trim(landdir) // trim(cyear)) + + IF (p_is_master) THEN + write(*,'(A,I4)') 'Aggregate LAI : ', iy + ENDIF IF (p_is_io) THEN CALL read_5x5_data_pft (dir_5x5, suffix, gridlai, 'PCT_PFT', pftPCT) diff --git a/share/MOD_Mesh.F90 b/share/MOD_Mesh.F90 index 14c4f1ff..54ea8704 100644 --- a/share/MOD_Mesh.F90 +++ b/share/MOD_Mesh.F90 @@ -695,7 +695,7 @@ SUBROUTINE mesh_build () idest = gblock%pio (meshtmp(ie)%xblk, meshtmp(ie)%yblk) ! send(09) - elmtag = meshtmp(ie)%indx + elmtag = mod(meshtmp(ie)%indx, 30000) smesg(1:5) = (/p_iam_glb, elmtag, meshtmp(ie)%xblk, meshtmp(ie)%yblk, meshtmp(ie)%npxl/) CALL mpi_send (smesg(1:5), 5, MPI_INTEGER, idest, mpi_tag_mesg, p_comm_glb, p_err) diff --git a/share/MOD_NetCDFSerial.F90 b/share/MOD_NetCDFSerial.F90 index c2f4a762..0e0e14db 100644 --- a/share/MOD_NetCDFSerial.F90 +++ b/share/MOD_NetCDFSerial.F90 @@ -590,6 +590,7 @@ SUBROUTINE ncio_read_serial_int32_2d (filename, dataname, rdata) ! Local variables integer :: ncid, varid integer, allocatable :: varsize(:) + integer :: dsp, nread CALL check_ncfile_exist (filename) @@ -598,7 +599,19 @@ SUBROUTINE ncio_read_serial_int32_2d (filename, dataname, rdata) CALL nccheck( nf90_open(trim(filename), NF90_NOWRITE, ncid) ) CALL nccheck( nf90_inq_varid(ncid, trim(dataname), varid) ) - CALL nccheck( nf90_get_var(ncid, varid, rdata) ) + + IF ((varsize(1) > 1000) .and. (varsize(2) > 100000)) THEN + dsp = 0 + DO WHILE (dsp < varsize(2)) + nread = min(100000,varsize(2)-dsp) + CALL nccheck (nf90_get_var(ncid, varid, & + rdata(1:varsize(1),dsp+1:dsp+nread), (/1,dsp+1/), (/varsize(1),nread/))) + dsp = dsp + nread + ENDDO + ELSE + CALL nccheck( nf90_get_var(ncid, varid, rdata) ) + ENDIF + CALL nccheck( nf90_close(ncid) ) deallocate (varsize) @@ -648,6 +661,7 @@ SUBROUTINE ncio_read_serial_real8_2d (filename, dataname, rdata) ! Local variables integer :: ncid, varid integer, allocatable :: varsize(:) + integer :: dsp, nread CALL check_ncfile_exist (filename) @@ -656,7 +670,19 @@ SUBROUTINE ncio_read_serial_real8_2d (filename, dataname, rdata) CALL nccheck( nf90_open(trim(filename), NF90_NOWRITE, ncid) ) CALL nccheck( nf90_inq_varid(ncid, trim(dataname), varid) ) - CALL nccheck( nf90_get_var(ncid, varid, rdata) ) + + IF ((varsize(1) > 1000) .and. (varsize(2) > 100000)) THEN + dsp = 0 + DO WHILE (dsp < varsize(2)) + nread = min(100000,varsize(2)-dsp) + CALL nccheck (nf90_get_var(ncid, varid, & + rdata(1:varsize(1),dsp+1:dsp+nread), (/1,dsp+1/), (/varsize(1),nread/))) + dsp = dsp + nread + ENDDO + ELSE + CALL nccheck( nf90_get_var(ncid, varid, rdata) ) + ENDIF + CALL nccheck( nf90_close(ncid) ) deallocate (varsize)