!******************************************************************************* !Subroutine - rapid_routing !******************************************************************************* subroutine rapid_routing(ZV_C1,ZV_C2,ZV_C3,ZV_Qext, & ZV_QoutinitR,ZV_VinitR, & ZV_QoutR,ZV_QoutbarR,ZV_VR,ZV_VbarR) !Purpose: !Performs flow calculation in each reach of a river network using the Muskingum !method (McCarthy 1938). Also calculates the volume of each reach using a !simple first order approximation !Author: !Cedric H. David, 2008-2015. !******************************************************************************* !Declaration of variables !******************************************************************************* use netcdf use rapid_var, only : & ZS_dtR,IS_R,JS_R, & ZM_Net,ZM_TC1, & ZV_b,ZV_babsmax,ZV_bhat, & ZV_QoutprevR,ZV_VprevR,ZV_QoutRabsmin,ZV_QoutRabsmax, & ZV_QoutRhat, & ZV_VoutR,ZV_Vext, & ierr,ksp, & ZS_one,IS_ksp_iter,IS_ksp_iter_max, & vecscat,ZV_SeqZero,ZV_pointer,rank, & IS_nc_status,IS_nc_id_fil_Qout,IS_nc_id_var_Qout, & IV_nc_start,IV_nc_count2, & IS_riv_bas,JS_riv_bas,IM_index_up, & IS_opt_routing,IV_nbup,IV_riv_index, & BS_opt_influence implicit none !******************************************************************************* !Includes !******************************************************************************* #include "finclude/petscsys.h" !base PETSc routines #include "finclude/petscvec.h" #include "finclude/petscvec.h90" !vectors, and vectors in Fortran90 #include "finclude/petscmat.h" !matrices #include "finclude/petscksp.h" !Krylov subspace methods #include "finclude/petscpc.h" !preconditioners #include "finclude/petscviewer.h" !viewers (allows writing results in file for example) !******************************************************************************* !Intent (in/out), and local variables !******************************************************************************* Vec, intent(in) :: ZV_C1,ZV_C2,ZV_C3,ZV_Qext, & ZV_QoutinitR,ZV_VinitR Vec, intent(out) :: ZV_QoutR,ZV_QoutbarR Vec :: ZV_VR,ZV_VbarR PetscInt :: IS_localsize,JS_localsize PetscScalar, pointer :: ZV_QoutR_p(:),ZV_QoutprevR_p(:),ZV_QoutinitR_p(:), & ZV_QoutbarR_p(:),ZV_Qext_p(:),ZV_C1_p(:),ZV_C2_p(:), & ZV_C3_p(:),ZV_b_p(:), & ZV_babsmax_p(:),ZV_QoutRabsmin_p(:),ZV_QoutRabsmax_p(:) !******************************************************************************* !Get local sizes for vectors !******************************************************************************* call VecGetLocalSize(ZV_QoutR,IS_localsize,ierr) !******************************************************************************* !Set mean values to zero initialize QoutprevR with QoutinitR !******************************************************************************* call VecSet(ZV_QoutbarR,0*ZS_one,ierr) !Qoutbar=0 !call VecSet(ZV_VbarR,0*ZS_one,ierr) !Vbar=0 !set the means to zero at beginning of iterations over routing time step call VecCopy(ZV_QoutinitR,ZV_QoutprevR,ierr) !QoutprevR=QoutinitR !call VecCopy(ZV_VinitR,ZV_VprevR,ierr) !VprevR=VinitR !set the previous value to the initial value given as input to subroutine !******************************************************************************* !Temporal loop !******************************************************************************* call VecGetArrayF90(ZV_C1,ZV_C1_p,ierr) call VecGetArrayF90(ZV_C2,ZV_C2_p,ierr) call VecGetArrayF90(ZV_C3,ZV_C3_p,ierr) call VecGetArrayF90(ZV_Qext,ZV_Qext_p,ierr) do JS_R=1,IS_R !------------------------------------------------------------------------------- !Update mean !------------------------------------------------------------------------------- call VecAXPY(ZV_QoutbarR,ZS_one/IS_R,ZV_QoutprevR,ierr) !Qoutbar=Qoutbar+Qoutprev/IS_R !call VecAXPY(ZV_VbarR,ZS_one/IS_R,ZV_VprevR,ierr) !Vbar=Vbar+Vprev/IS_R !------------------------------------------------------------------------------- !Calculation of the right hand size, b !------------------------------------------------------------------------------- call MatMult(ZM_Net,ZV_QoutprevR,ZV_b,ierr) !b2=Net*Qoutprev call VecGetArrayF90(ZV_b,ZV_b_p,ierr) call VecGetArrayF90(ZV_QoutprevR,ZV_QoutprevR_p,ierr) do JS_localsize=1,IS_localsize ZV_b_p(JS_localsize)=ZV_b_p(JS_localsize)*ZV_C2_p(JS_localsize) & +(ZV_C1_p(JS_localsize)+ZV_C2_p(JS_localsize)) & *ZV_Qext_p(JS_localsize) & +ZV_C3_p(JS_localsize)*ZV_QoutprevR_p(JS_localsize) end do call VecRestoreArrayF90(ZV_QoutprevR,ZV_QoutprevR_p,ierr) call VecRestoreArrayF90(ZV_b,ZV_b_p,ierr) !------------------------------------------------------------------------------- !Routing with PETSc using a matrix method !------------------------------------------------------------------------------- if (IS_opt_routing==1) then call KSPSolve(ksp,ZV_b,ZV_QoutR,ierr) !solves A*Qout=b call KSPGetIterationNumber(ksp,IS_ksp_iter,ierr) if (IS_ksp_iter>IS_ksp_iter_max) IS_ksp_iter_max=IS_ksp_iter end if !------------------------------------------------------------------------------- !Routing with Fortran using the traditional Muskingum method !------------------------------------------------------------------------------- if (IS_opt_routing==2) then call VecGetArrayF90(ZV_QoutR,ZV_QoutR_p,ierr) call VecGetArrayF90(ZV_QoutprevR,ZV_QoutprevR_p,ierr) call VecGetArrayF90(ZV_b,ZV_b_p,ierr) do JS_riv_bas=1,IS_riv_bas ZV_QoutR_p(JS_riv_bas)=ZV_b_p(JS_riv_bas) & +sum(ZV_C1_p(JS_riv_bas) & *ZV_QoutR_p(IM_index_up(JS_riv_bas,1: & IV_nbup(IV_riv_index(JS_riv_bas))))) end do !Taking into account the knowledge of how many upstream locations exist. !Similar to exact preallocation of network matrix call VecRestoreArrayF90(ZV_QoutR,ZV_QoutR_p,ierr) call VecRestoreArrayF90(ZV_QoutprevR,ZV_QoutprevR_p,ierr) call VecRestoreArrayF90(ZV_b,ZV_b_p,ierr) end if !------------------------------------------------------------------------------- !Routing with PETSc using a matrix method with transboundary matrix !------------------------------------------------------------------------------- if (IS_opt_routing==3) then call KSPSolve(ksp,ZV_b,ZV_QoutRhat,ierr) !solves A*Qouthat=b call KSPGetIterationNumber(ksp,IS_ksp_iter,ierr) if (IS_ksp_iter>IS_ksp_iter_max) IS_ksp_iter_max=IS_ksp_iter call MatMult(ZM_TC1,ZV_QoutRhat,ZV_bhat,ierr) call VecAYPX(ZV_bhat,ZS_one,ZV_b,ierr) call KSPSolve(ksp,ZV_bhat,ZV_QoutR,ierr) !solves A*Qout=bhat call KSPGetIterationNumber(ksp,IS_ksp_iter,ierr) if (IS_ksp_iter>IS_ksp_iter_max) IS_ksp_iter_max=IS_ksp_iter end if !------------------------------------------------------------------------------- !Calculation of babsmax, QoutRabsmin and QoutRabsmax !------------------------------------------------------------------------------- if (BS_opt_influence) then call VecGetArrayF90(ZV_b,ZV_b_p,ierr) call VecGetArrayF90(ZV_babsmax,ZV_babsmax_p,ierr) do JS_localsize=1,IS_localsize if (ZV_babsmax_p(JS_localsize)<=abs(ZV_b_p(JS_localsize))) then ZV_babsmax_p(JS_localsize) =abs(ZV_b_p(JS_localsize)) end if end do call VecRestoreArrayF90(ZV_b,ZV_b_p,ierr) call VecRestoreArrayF90(ZV_babsmax,ZV_babsmax_p,ierr) call VecGetArrayF90(ZV_QoutR,ZV_QoutR_p,ierr) call VecGetArrayF90(ZV_QoutRabsmin,ZV_QoutRabsmin_p,ierr) call VecGetArrayF90(ZV_QoutRabsmax,ZV_QoutRabsmax_p,ierr) do JS_localsize=1,IS_localsize if (ZV_QoutRabsmin_p(JS_localsize)>=abs(ZV_QoutR_p(JS_localsize))) then ZV_QoutRabsmin_p(JS_localsize) =abs(ZV_QoutR_p(JS_localsize)) end if if (ZV_QoutRabsmax_p(JS_localsize)<=abs(ZV_QoutR_p(JS_localsize))) then ZV_QoutRabsmax_p(JS_localsize) =abs(ZV_QoutR_p(JS_localsize)) end if end do call VecRestoreArrayF90(ZV_QoutR,ZV_QoutR_p,ierr) call VecRestoreArrayF90(ZV_QoutRabsmin,ZV_QoutRabsmin_p,ierr) call VecRestoreArrayF90(ZV_QoutRabsmax,ZV_QoutRabsmax_p,ierr) end if !------------------------------------------------------------------------------- !Calculation of V (this part can be commented to accelerate parameter !estimation in calibration mode) !------------------------------------------------------------------------------- !call VecCopy(ZV_QoutR,ZV_VoutR,ierr) !Vout=Qout !call VecScale(ZV_VoutR,ZS_dtR,ierr) !Vout=Vout*dt !!result Vout=Qout*dt ! !call VecCopy(ZV_Qext,ZV_Vext,ierr) !Vext=Qext !call VecScale(ZV_Vext,ZS_dtR,ierr) !Vext=Vext*dt !!result Vext=Qext*dt ! !call MatMult(ZM_Net,ZV_VoutR,ZV_VR,ierr) !V=Net*Vout !call VecAXPY(ZV_VR,ZS_one,ZV_Vext,ierr) !V=V+Vext !call VecAXPY(ZV_VR,-ZS_one,ZV_VoutR,ierr) !V=V-Vout !call VecAXPY(ZV_VR,ZS_one,ZV_VprevR,ierr) !V=V+Vprev !!result V=Vprev+(Net*Vout+Vext)-Vout !------------------------------------------------------------------------------- !Reset previous !------------------------------------------------------------------------------- call VecCopy(ZV_QoutR,ZV_QoutprevR,ierr) !Qoutprev=Qout !call VecCopy(ZV_VR,ZV_VprevR,ierr) !Vprev=V !reset previous !------------------------------------------------------------------------------- !optional write outputs !------------------------------------------------------------------------------- !call VecScatterBegin(vecscat,ZV_QoutR,ZV_SeqZero, & ! INSERT_VALUES,SCATTER_FORWARD,ierr) !call VecScatterEnd(vecscat,ZV_QoutR,ZV_SeqZero, & ! INSERT_VALUES,SCATTER_FORWARD,ierr) !call VecGetArrayF90(ZV_SeqZero,ZV_pointer,ierr) !!if (rank==0) write (99,'(10e10.3)') ZV_pointer !if (rank==0) IS_nc_status=NF90_PUT_VAR(IS_nc_id_fil_Qout,IS_nc_id_var_Qout, & ! ZV_pointer, & ! [IV_nc_start(1),(IV_nc_start(2)-1)*IS_R+JS_R],IV_nc_count2) !call VecRestoreArrayF90(ZV_SeqZero,ZV_pointer,ierr) !------------------------------------------------------------------------------- !End temporal loop !------------------------------------------------------------------------------- end do call VecRestoreArrayF90(ZV_C1,ZV_C1_p,ierr) call VecRestoreArrayF90(ZV_C2,ZV_C2_p,ierr) call VecRestoreArrayF90(ZV_C3,ZV_C3_p,ierr) call VecRestoreArrayF90(ZV_Qext,ZV_Qext_p,ierr) !******************************************************************************* !End !******************************************************************************* end subroutine rapid_routing