module RvicMod use shr_kind_mod , only : r8 => shr_kind_r8, SHR_KIND_CL use shr_file_mod , only : shr_file_getUnit, shr_file_freeunit !use shr_const_mod , only : SHR_CONST_RHOFW use RvicVar use RvicRestFile use netcdf use rvic_cpl_indices, only : nt_rvic, nt_uh, rvic_tracers use RvicHistFlds , only : RvicHistFldsInit, RvicHistFldsSet use RvicHistFile , only : RvicHistUpdateHbuf, RvicHistHtapesWrapup, & RvicHistHtapesBuild, rvichist_ndens, & rvichist_mfilt, rvichist_nhtfrq, & rvichist_avgflag_pertape, rvichist_fincl1, & rvichist_fincl2, rvichist_fincl3, & rvichist_fexcl1, rvichist_fexcl2, & rvichist_fexcl3, max_tapes, max_namlen use RvicSpmd , only : masterproc, npes, iam, mpicom_rof, & MPI_REAL8,MPI_INTEGER,MPI_CHARACTER,MPI_LOGICAL,MPI_MAX use RvicTimeManager , only : timemgr_init, get_nstep, get_curr_date, RvicTimeMode use perf_mod use RvicIO use mct_mod implicit none SAVE private ! PUBLIC MEMBER FUNCTIONS: public :: rvic_init public :: rvic_run ! PUBLIC DATA MEMBERS: ! None ! PRIVATE MEMBER FUNCTIONS: private :: RvicParamFileRead ! rvic read parameters routine private :: RvicGridFileRead ! rvic read grid file routine private :: RvicConvolve private :: RvicUpdate private :: handle_err ! PRIVATE DATA MEMBERS: ! None ! REVISION HISTORY: ! 07/21/2013 ! Author: Joe Hamman, University of Washington ! Comment: Initial RVIC model ! !=============================================================== contains !=============================================================== !----------------------------------------------------------------------- subroutine rvic_init(NLFilename, rvic_active) !---------------------------------------------------------------------------! ! DESCRIPTION: ! Initialize RVIC variables ! Intended to be called from rof_comp_mct.F90 ! ! ARGUMENTS: implicit none character(len=*), intent(in) :: NLFilename logical, intent(out) :: rvic_active ! LOCAL VARIABLES: logical :: rvar_set ! status of runoff variable logical :: do_rvic, ice_runoff integer :: nunit ! unit number integer :: ier ! error code integer :: i, j, n ! Timestep Info integer :: rvic_tstep ! model timestep ! Input file names character(len=SHR_KIND_CL) :: gridfile ! input grid file character(len=SHR_KIND_CL) :: paramfile ! input parameter file character(len=256):: finidat_rvic ! input state file character(len=256):: pnamer ! full pathname of netcdf restart file character(len=256):: fnamer ! name of netcdf restart file ! Subroutine descriptions character(len=*), parameter :: subname = 'rvic_init' character(len=*), parameter :: format = "('("//trim(subname)//") :',A)" !----- define namelist -----! namelist / rvic_inparm / rvic_tstep, gridfile, paramfile, & finidat_rvic, nrevsn_rvic, rvichist_ndens, rvichist_mfilt, & rvichist_nhtfrq, rvichist_fincl1, rvichist_fincl2, rvichist_fincl3, & rvichist_fexcl1, rvichist_fexcl2, rvichist_fexcl3, & rvichist_avgflag_pertape !--------------------------------------------------------- ! Read rvic_inparam namelist !--------------------------------------------------------- rvic_tstep = 0 gridfile = ' ' paramfile = ' ' finidat_rvic = ' ' nunit = shr_file_getunit() open (nunit,file='rvic_in',status="old",action="read") read (nunit,nml=rvic_inparm,iostat=ier) close(nunit) call shr_file_freeunit(nunit) write(iulog,format) 'FINISHED READING NAMELIST' write(iulog,format) '--------------------- NAMELIST VARIABLES ---------------------' write(iulog,*) 'rvic_tstep: ', rvic_tstep write(iulog,*) 'gridfile: ', gridfile write(iulog,*) 'paramfile: ', paramfile write(iulog,*) 'finidat_rvic: ', finidat_rvic write(iulog,*) 'nrevsn_rvic: ', nrevsn_rvic write(iulog,*) 'rvichist_ndens: ', rvichist_ndens write(iulog,*) 'rvichist_mfilt: ', rvichist_mfilt write(iulog,*) 'rvichist_nhtfrq: ', rvichist_nhtfrq write(iulog,*) 'rvichist_fincl1: ', rvichist_fincl1 write(iulog,*) 'rvichist_fincl2: ', rvichist_fincl2 write(iulog,*) 'rvichist_fincl3: ', rvichist_fincl3 write(iulog,*) 'rvichist_fexcl1: ', rvichist_fexcl1 write(iulog,*) 'rvichist_fexcl2: ', rvichist_fexcl2 write(iulog,*) 'rvichist_fexcl3: ', rvichist_fexcl3 write(iulog,*) 'rvichist_avgflag_pertape: ', rvichist_avgflag_pertape write(iulog,format) '------------------ END OF NAMELIST VARIABLES -----------------' do_rvic = .true. ice_runoff = .true. call mpi_bcast (rvic_tstep, 1, MPI_INTEGER, 0, mpicom_rof, ier) call mpi_bcast (do_rvic, 1, MPI_LOGICAL, 0, mpicom_rof, ier) call mpi_bcast (ice_runoff, 1, MPI_LOGICAL, 0, mpicom_rof, ier) call mpi_bcast (finidat_rvic, len(finidat_rvic), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (nrevsn_rvic, len(nrevsn_rvic), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_nhtfrq, size(rvichist_nhtfrq), MPI_INTEGER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_mfilt , size(rvichist_mfilt) , MPI_INTEGER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_ndens , size(rvichist_ndens) , MPI_INTEGER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_fexcl1, (max_namlen+2)*size(rvichist_fexcl1), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_fexcl2, (max_namlen+2)*size(rvichist_fexcl2), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_fexcl3, (max_namlen+2)*size(rvichist_fexcl3), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_fincl1, (max_namlen+2)*size(rvichist_fincl1), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_fincl2, (max_namlen+2)*size(rvichist_fincl2), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_fincl3, (max_namlen+2)*size(rvichist_fincl3), MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (rvichist_avgflag_pertape, size(rvichist_avgflag_pertape), MPI_CHARACTER, 0, mpicom_rof, ier) call RvicVarInit() !------------------------------------------------------- ! Initialize rvic time manager !------------------------------------------------------- ! Intiialize RVIC pio call ncd_pio_init() ! Obtain restart file if appropriate if ((nsrest == nsrContinue) .or. (nsrest == nsrBranch )) then call RvicRestGetfile( file=fnamer, path=pnamer ) finidat_rvic = fnamer else fnamer = finidat_rvic endif ! Initialize time manager if (nsrest == nsrStartup) then call timemgr_init(dtime_in=rvic_tstep) else call RvicRestTimeManager(file=fnamer) end if !------------------------------------------------------- ! Read Grid File data !------------------------------------------------------- call RvicGridFileRead(trim(gridfile) ,rviclon, rviclat) write(iulog, *) subname, 'rviclat=',rviclat, 'rviclon=', rviclon ! ------------------------------------------------------- !------------------------------------------------------- ! Read Param File data !------------------------------------------------------- call RvicParamFileRead(trim(paramfile)) if ((Rvar%nTracers /= nt_rvic).and.(Rvar%nTracers==1)) then write(iulog, format) 'using the liquid uh for all tracers' do n=1,nt_rvic nt_uh(n) = 1 end do elseif (Rvar%nTracers==nt_rvic) then write(iulog, format) 'Number of tracers in parameter file matches rvic_cpl_indices nt_rvic' else write(iulog, format) 'FATAL ERROR: mismatch between tracers and uh shape' write(iulog, *) 'nt_rvic=',nt_rvic, 'nTracers=',Rvar%nTracers call shr_sys_abort(subname) end if ! ------------------------------------------------------- ! Allocate runoff datatype and flux containers ! ------------------------------------------------------- call RunoffInit(1, rviclon*rviclat, rviclon*rviclat, rviclon*rviclat) call FluxContainInit(1, rviclon*rviclat, rviclon*rviclat) !------------------------------------------------------- ! Initialize runoff data type !------------------------------------------------------- n=0 do i=1,rviclon do j=1,rviclat n=n+1 runoff%lonc(n) = lon(i,j) runoff%latc(n) = lat(i,j) runoff%area(n) = area(i,j) runoff%gindex(n) = n end do end do runoff%rlon = runoff%lonc runoff%rlat = runoff%latc runoff%lnumr = rviclon*rviclat !------------------------------------------------------- ! Read restart/initial info !------------------------------------------------------- write(iulog,format) 'Getting initial state data' write(iulog, *) 'initial state file: ',finidat_rvic ! The call below opens and closes the file if ((nsrest == nsrContinue) .or. (nsrest == nsrBranch)) then if (finidat_rvic /= ' ') then call RvicRestFileRead( file=fnamer ) else write(iulog, format) 'Run is continue or branch and no restart file was found' call shr_sys_abort(subname) end if elseif ((nsrest == nsrStartup) .and. (finidat_rvic /= ' ')) then call RvicRestFileRead( file=fnamer ) elseif ((nsrest == nsrStartup) .and. (finidat_rvic == ' ')) then write(iulog, format) 'Run is a drystart (no init data provided)' else write(iulog, format) 'Problem determining runtype or initial state file' call shr_sys_abort(subname) end if ! ------------------------------------------------------- call RvicRvarSet(rvar_set) if (.not.rvar_set) then call shr_sys_abort(subname//' FATAL ERROR: checking runoff var returned an error') end if ! ------------------------------------------------------- ! Update the runoff type with initial conditions ! ------------------------------------------------------- call RvicUpdate() ! ------------------------------------------------------- ! get rvic agg_tsteps ! ------------------------------------------------------- call RvicTimeMode(rvic_tstep, rvar%unit_hydrograph_dt, rvar%agg_tsteps) !------------------------------------------------------- ! Initialize RVIC history handler and fields !------------------------------------------------------- call RvicHistFldsInit() if (nsrest==nsrStartup .or. nsrest==nsrBranch) then call RvicHistHtapesBuild() end if call RvicHistFldsSet() rvic_active = do_rvic write(iulog,*) 'done with rvic_init' end subroutine rvic_init subroutine rvic_run(rdate, rstwr, nlend) !--------------------------------------------------------------------------- ! DESCRIPTION: ! Run the RVIC model ! ! ARGUMENTS: character(len=*), intent(in) :: rdate ! LOCAL VARIABLES logical :: rstwr logical :: nlend character(len=72) :: state_fname character(len=*), parameter :: subname = 'rvic_run' character(len=*), parameter :: format = "('("//trim(subname)//") :',A)" !--------------------------------------------------------------------------- !--------------------------------------------------------------------------- ! Accumulate the incoming flux (totrunin) into agg_in !--------------------------------------------------------------------------- call t_startf ('rvic_run_convolve') Rvar%agg_in = Rvar%agg_in + Rvar%point_in Rvar%agg_counter = Rvar%agg_counter + 1 if (minval(Rvar%point_in) < 0.0) then write(iulog, *) subname,'minval point_in = ',minval(Rvar%point_in) write(iulog, *) subname,'maxval point_in = ',maxval(Rvar%point_in) write(iulog, *) subname,'ERROR: incoming runoff is negative!' ! call shr_sys_abort(subname//' ERROR: incoming runoff is negative!') endif ! if it is time, convolve, the fluxes, otherwise wait if (Rvar%agg_counter == Rvar%agg_tsteps) then Rvar%agg_in = Rvar%agg_in/Rvar%agg_tsteps call RvicConvolve() Rvar%agg_counter = 0 Rvar%agg_in = 0.0_r8 else write(iulog, *) "Agg_counter is at ", Rvar%agg_counter, " of ", Rvar%agg_tsteps write(iulog, *) "Waiting to convolve..." end if call t_stopf ('rvic_run_convolve') !--------------------------------------------------------------------------- ! BUDGET, compute global sums and check ! if (budget_check) then ! end if !--------------------------------------------------------------------------- call t_startf ('rvic_run_write') !--------------------------------------------------------------------------- ! Write RVIC History File !--------------------------------------------------------------------------- call t_startf('RVIChbuf') call RvicHistFldsSet() call RvicHistUpdateHbuf() call t_stopf('RVIChbuf') call t_startf('RVIChtapes') call RvicHistHtapesWrapup( rstwr, nlend ) call t_stopf('RVIChtapes') !--------------------------------------------------------------------------- ! Write the current state (restart) file !--------------------------------------------------------------------------- if (rstwr) then call t_startf('RVICrest') call RvicRestFileWrite( rdate ) call t_stopf('RVICrest') end if call t_stopf ('rvic_run_write') !--------------------------------------------------------------------------- ! Wrap up !--------------------------------------------------------------------------- write(iulog, format) 'Done with rvic_run for '//rdate end subroutine rvic_run subroutine RvicParamFileRead(paramfile) ! arguments character(len=*), intent(in) :: paramfile ! RVIC parameter file ! local variables integer :: ncId ! netCDF ID integer :: timeDimId ! netCDF time dimension ID integer :: outletDimID ! netCDF outlet dimension ID integer :: sourceDimID ! netCDF point dimension ID integer :: tracersDimID ! netCDF tracer dimension ID integer :: nTimes ! number of time steps in subset integer :: VarID ! netCDF variable ID integer :: status ! netCDF function status integer :: i,n,s ! loop counters real(r8),pointer :: uh_read(:,:,:) ! tmp for reading uh character(len=*), parameter :: subname = "RvicParamFileRead" character(len=*), parameter :: format = "('("//trim(subname)//"):',A)" write(iulog,format) 'Reading input file now and populating Rvar ...' ! open parameter netcdf status = nf90_open(paramfile, nf90_nowrite, ncId) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Opened parameter file: ",trim(paramfile) ! get dimension sizes and inds status = nf90_inq_dimid(ncId, "timesteps", timeDimID) if (status /= nf90_noerr) call handle_err(status) status = nf90_inquire_dimension(ncId, timeDimID, len = nTimes) if (status /= nf90_noerr) call handle_err(status) write(iulog,*) "Size of time dimension: ",nTimes status = nf90_inq_dimid(ncId, "outlets", outletDimID) if (status /= nf90_noerr) call handle_err(status) status = nf90_inquire_dimension(ncId, outletDimID, len = Rvar%nOutlets) if (status /= nf90_noerr) call handle_err(status) write(iulog,*) "Number of Outlet Points: ",Rvar%nOutlets status = nf90_inq_dimid(ncId, "sources", sourceDimID) if (status /= nf90_noerr) call handle_err(status) status = nf90_inquire_dimension(ncId, sourceDimID, len = Rvar%nSources) if (status /= nf90_noerr) call handle_err(status) write(iulog,*) "Number of Source Points: ",Rvar%nSources status = nf90_inq_dimid(ncId, "tracers", tracersDimID) if (status /= nf90_noerr) call handle_err(status) status = nf90_inquire_dimension(ncId, tracersDimID, len = Rvar%nTracers) if (status /= nf90_noerr) call handle_err(status) write(iulog,*) "Number of tracers in param file: ",Rvar%nTracers ! Allocate Memory for incoming variables write(iulog,format) "Allocating memory..." allocate(Rvar%unit_hydrograph(nTimes, Rvar%nSources, Rvar%ntracers), & Rvar%source_decomp_ind(Rvar%nSources), & Rvar%source_x_ind(Rvar%nSources), & Rvar%source_y_ind(Rvar%nSources), & Rvar%source_time_offset(Rvar%nSources), & Rvar%source2outlet_ind(Rvar%nSources), & Rvar%outlet_decomp_ind(Rvar%nOutlets), & Rvar%outlet_x_ind(Rvar%nOutlets), & Rvar%outlet_y_ind(Rvar%nOutlets), & Rvar%outlet_number(Rvar%nOutlets), & Rvar%outlet_mask(Rvar%nOutlets)) write(iulog,format) "Reading variables from parameter file..." ! get variables, need to swap order of dimensions allocate(uh_read(Rvar%nTracers, Rvar%nSources, nTimes)) status = nf90_inq_varid(ncId, "unit_hydrograph", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, uh_read) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read unit_hydrograph" do i = 1,nTimes do s = 1,Rvar%nSources do n = 1,Rvar%ntracers Rvar%unit_hydrograph(i,s,n) = uh_read(n,s,i) enddo enddo enddo deallocate(uh_read) status = nf90_inq_varid(ncId, "source_decomp_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%source_decomp_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read source_decomp_ind" status = nf90_inq_varid(ncId, "source_x_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%source_x_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read source_x_ind" status = nf90_inq_varid(ncId, "source_y_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%source_y_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read source_y_ind" status = nf90_inq_varid(ncId, "source_time_offset", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%source_time_offset) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read source_time_offset" status = nf90_inq_varid(ncId, "source2outlet_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%source2outlet_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read source2outlet_ind" status = nf90_inq_varid(ncId, "outlet_decomp_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%outlet_decomp_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read outlet_decomp_ind" status = nf90_inq_varid(ncId, "outlet_x_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%outlet_x_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read outlet_x_ind" status = nf90_inq_varid(ncId, "outlet_y_ind", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%outlet_y_ind) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read outlet_y_ind" status = nf90_inq_varid(ncId, "outlet_mask", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%outlet_mask) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read outlet_mask" status = nf90_inq_varid(ncId, "outlet_number", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%outlet_number) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read outlet_number" status = nf90_inq_varid(ncId, "subset_length", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%subset_length) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read subset_length: ",Rvar%subset_length status = nf90_inq_varid(ncId, "full_time_length", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%full_time_length) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read full_length: ",Rvar%full_time_length status = nf90_inq_varid(ncId, "unit_hydrograph_dt", VarID) if(status /= nf90_NoErr) call handle_err(status) status = nf90_get_var(ncId, VarID, Rvar%unit_hydrograph_dt) if(status /= nf90_NoErr) call handle_err(status) write(iulog,format) "Read timestep: ",Rvar%unit_hydrograph_dt ! close the parameter file status = nf90_close(ncId) if (status /= nf90_noerr) call handle_err(status) write(iulog,format) 'Done reading parameter file, Rvar parameter variables are populated...' ! allocate memory for state variables write(iulog,format) "Allocating memory for state variables" allocate(Rvar%ring( Rvar%full_time_length, Rvar%nOutlets, nt_rvic ), & Rvar%ring_timesteps( Rvar%full_time_length ), & Rvar%upstrm_storage( Rvar%nOutlets, nt_rvic ), & Rvar%upstrm_dstorage( Rvar%nOutlets, nt_rvic )) Rvar%ring = 0.0_r8 Rvar%ring_timesteps = 0 Rvar%upstrm_storage = 0.0_r8 Rvar%upstrm_dstorage = 0.0_r8 end subroutine RvicParamFileRead !=============================================================================== subroutine RvicGridFileRead(gridfile, nxg, nyg) !----------------------------------------------------- ! ! !DESCRIPTION: ! Read the grid file, set module data ! assume nxg,nyg,lon,lat,area,mask are module data ! ! !ARGUMENTS: use netcdf implicit none character(len=*),intent(in) :: gridfile ! input grid file integer, intent(out) :: nxg, nyg ! size of x and y grid dims !----------------------------------------------------- ! !LOCAL: integer :: ncid, vid, dids(2), ndims, rcode character(len=*),parameter :: subname = "RvicGridFileRead" !----------------------------------------------------- rcode = nf90_open(trim(gridfile),nf90_nowrite,ncid) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_open file '//trim(gridfile)) rCode = nf90_inq_varid(ncid,'LONGXY',vid) if (rcode /= nf90_noerr) rCode = nf90_inq_varid(ncid,'xc',vid) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inq_varid longxy'//trim(gridfile)) rCode = nf90_inquire_variable(ncid,vid,ndims=ndims) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_variable ndims') rCode = nf90_inquire_variable(ncid,vid,dimids=dids) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_variable dimids') rCode = nf90_inquire_dimension(ncid,dids(1),len=nxg) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension nxg') rCode = nf90_inquire_dimension(ncid,dids(2),len=nyg) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inquire_dimension nyg') allocate(lon(nxg,nyg),lat(nxg,nyg),area(nxg,nyg),mask(nxg,nyg)) rCode = nf90_inq_varid(ncid,'LONGXY',vid) if (rcode /= nf90_noerr) rCode = nf90_inq_varid(ncid,'xc',vid) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inq_varid lon') rCode = nf90_get_var(ncid,vid,lon) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var lon') rCode = nf90_inq_varid(ncid,'LATIXY',vid) if (rcode /= nf90_noerr) rCode = nf90_inq_varid(ncid,'yc',vid) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inq_varid lat') rCode = nf90_get_var(ncid,vid,lat) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var lat') rCode = nf90_inq_varid(ncid,'AREA',vid) if (rcode /= nf90_noerr) rCode = nf90_inq_varid(ncid,'area',vid) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_inq_varid area') rCode = nf90_get_var(ncid,vid,area) if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var area') ! rCode = nf90_inq_varid(ncid,'MASK',vid) ! if (rcode /= nf90_noerr) rCode = nf90_inq_varid(ncid,'mask',vid) ! if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: ! nf90_inq_varid mask') ! rCode = nf90_get_var(ncid,vid,mask) ! if (rcode /= nf90_noerr) call shr_sys_abort(subname//' ERROR: nf90_get_var ! mask') mask = 1.0_r8 write(iulog,*) subname,' min/max lon ',minval(lon) ,maxval(lon) write(iulog,*) subname,' min/max lat ',minval(lat) ,maxval(lat) write(iulog,*) subname,' min/max area',minval(area),maxval(area) write(iulog,*) subname,' min/max mask',minval(mask),maxval(mask) end subroutine RvicGridFileRead subroutine RvicConvolve() !--------------------------------------------------------------------------- ! DESCRIPTION: ! Convolution of the unit hydrograph (Rvar%unit_hydrograph) and the liquid flux ! (agg_in). The basics steps of the convolution are as follows: ! ! 0. Rotate the ring after zeroing out the first column ! 1. Loop over all source points, convolving the flux at that source point ! with the unit_hydrograph of that source point. ! 2. Add the convolution convolved hydrograph to the ring (Rvar%ring) ! for the source points corresponding outlet point. ! ! ARGUMENTS: ! None ! Local Variables integer :: d, s, i, j, outlet, nt, uhnt integer, pointer :: mask call t_startf ('rvic_convolve_shift') ! print out some debugging info write(iulog, *) "Doing convolution now" write(iulog, *) "Aggrunin: ", minval(Rvar%agg_in), maxval(Rvar%agg_in) write(iulog, *) "Ring before: ", minval(Rvar%ring), maxval(Rvar%ring) ! Do the convolution ! First freshen the convolution ring ! Rvar%ring(1, :, :) = 0.0_r8 ! Zero out current ring ! Rvar%ring = cshift(Rvar%ring, SHIFT = 1, DIM = 1) do nt = 1, nt_rvic do s = 1, Rvar%nOutlets do j = 1, Rvar%full_time_length-1 Rvar%ring(j,s,nt) = Rvar%ring(j+1,s,nt) enddo Rvar%ring(Rvar%full_time_length,s,nt) = 0.0_r8 enddo enddo call t_stopf ('rvic_convolve_shift') ! Now do the convolution call t_startf ('rvic_convolve_loop') do nt=1,nt_rvic uhnt = nt_uh(nt) do s=1,Rvar%nSources ! loop over all source points outlet = Rvar%source2outlet_ind(s) + 1 ! outlet position in ring d = Rvar%source_decomp_ind(s) + 1 ! decomp id of source point do i=1,Rvar%subset_length ! Loop over unit hydrograph j = i + Rvar%source_time_offset(s) ! Position in ring Rvar%ring(j,outlet,nt) = Rvar%ring(j,outlet,nt) + Rvar%unit_hydrograph(i, s, uhnt)*Rvar%agg_in(d,nt) end do end do end do call t_stopf ('rvic_convolve_loop') call t_startf ('rvic_convolve_update') call RvicUpdate() call t_stopf ('rvic_convolve_update') end subroutine RvicConvolve subroutine RvicUpdate() ! Update the runoff container with the fluxes from the ring ! Local Variables integer :: i, outlet, nt integer :: mask ! set all fields to zero Rvar%upstrm_dstorage = 0.0_r8 Runoff%Runoff = 0.0_r8 Runoff%dvolrdt = 0.0_r8 Runoff%volr = 0.0_r8 Runoff%runoffocn = 0.0_r8 Runoff%runofflnd = 0.0_r8 Runoff%dvolrdtlnd = 0.0_r8 Runoff%volrlnd = 0.0_r8 !Runoff%runoffguage = 0.0_r8 ! Get current storage Rvar%upstrm_dstorage = Rvar%upstrm_storage - sum(Rvar%ring, dim=1) Rvar%upstrm_storage = 0.0_r8 Rvar%upstrm_storage = sum(Rvar%ring, dim=1) ! Update the runoff fields do nt = 1,nt_rvic do i=1,Rvar%nOutlets outlet = Rvar%outlet_decomp_ind(i)+1 mask = Rvar%outlet_mask(i) Runoff%Runoff(outlet, nt) = Rvar%Ring(1, i, nt) Runoff%dvolrdt(outlet, nt) = Rvar%upstrm_dstorage(i, nt) Runoff%volr(outlet, nt) = Rvar%upstrm_storage(i, nt) if (mask == 0) then Runoff%runoffocn(outlet,nt) = Rvar%Ring(1, i, nt) else if (mask == 1) then Runoff%runofflnd(outlet,nt) = Rvar%Ring(1, i, nt) Runoff%dvolrdtlnd(outlet,nt) = Rvar%upstrm_dstorage(i, nt) Runoff%volrlnd(outlet,nt) = Rvar%upstrm_storage(i, nt) !else if (mask == 2) then ! Runoff%runoffguage(outlet,nt) = Rvar%Ring(1, i, nt) else write(iulog, *) 'unknown mask at outlet ', outlet end if end do end do end subroutine RvicUpdate subroutine handle_err(status) integer, intent ( in) :: status if(status /= nf90_noerr) then write(iulog,*) nf90_strerror(status) stop "Stopped due to error status" end if end subroutine handle_err !--------------------------------------------------------------------------- end module RvicMod