module river_mod !----------------------------------------------------------------------- ! GNU General Public License ! ! This program is free software; you can redistribute it and/or modify it and ! are expected to follow the terms of the GNU General Public License ! as published by the Free Software Foundation; either version 2 of ! the License, or (at your option) any later version. ! ! MOM is distributed in the hope that it will be useful, but WITHOUT ! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY ! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public ! License for more details. ! ! For the full text of the GNU General Public License, ! write to: Free Software Foundation, Inc., ! 675 Mass Ave, Cambridge, MA 02139, USA. ! or see: http://www.gnu.org/licenses/gpl.html !----------------------------------------------------------------------- ! Kirsten Findell ! Zhi Liang ! ! ! Processor domain layout for river model. If layout(1)*layout(2) is not equal ! to mpp_npes, the river model layout will be assigned the layout of land model ! passed through river_init. ! ! ! set true to run river model ( default is true). If FALSE, rivers are ! essentially turned off to save computing time ! ! ! slow time step for river model. dt_slow must be integer multiplier of dt_fast passed ! from land model ( land model time step). ! ! ! Number of slow time steps between sending out diagnositics data(default is 1). Please ! note that diagnostic output frequency ( specified in diag_table ) must be divided by ! diag_freq*dt_slow. ! ! #ifdef INTERNAL_FILE_NML use mpp_mod, only: input_nml_file #else use fms_mod, only: open_namelist_file #endif use mpp_mod, only : CLOCK_SUBCOMPONENT, CLOCK_ROUTINE use mpp_mod, only : mpp_error, mpp_chksum, FATAL, WARNING, NOTE, stdlog, mpp_npes use mpp_mod, only : mpp_pe, stdout, mpp_chksum, mpp_max use mpp_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end, MPP_CLOCK_DETAILED use mpp_domains_mod, only : domain2d, mpp_get_compute_domain, mpp_get_global_domain use mpp_domains_mod, only : mpp_get_data_domain, mpp_update_domains, mpp_get_ntile_count use fms_mod, only : write_version_number, check_nml_error, string use fms_mod, only : close_file, file_exist, field_size, read_data, write_data, lowercase use fms_mod, only : field_exist, CLOCK_FLAG_DEFAULT use fms_io_mod, only : get_mosaic_tile_file, get_instance_filename use diag_manager_mod, only : diag_axis_init, register_diag_field, register_static_field, send_data use time_manager_mod, only : time_type, increment_time, get_time use river_type_mod, only : river_type, Leo_Mad_trios, NO_RIVER_FLAG use river_physics_mod, only : river_physics_step, river_physics_init, river_impedes_lake, & river_impedes_large_lake use constants_mod, only : PI, RADIAN, tfreeze, DENS_H2O, hlf use stock_constants_mod, only : ISTOCK_WATER, ISTOCK_HEAT use land_tile_mod, only : land_tile_type, land_tile_enum_type, & first_elmt, tail_elmt, next_elmt, current_tile, get_elmt_indices, & operator(/=) use land_data_mod, only : land_data_type, land_state_type, lnd use lake_tile_mod, only : num_l implicit none private !--- version information --------------------------------------------- character(len=128) :: version = '$Id: river.F90,v 20.0 2013/12/13 23:29:41 fms Exp $' character(len=128) :: tagname = '$Name: tikal $' !--- public interface ------------------------------------------------ public :: river_init, river_end, river_type, update_river, river_stock_pe public :: save_river_restart !--- namelist interface ---------------------------------------------- logical :: do_rivers = .TRUE. ! if FALSE, rivers are essentially turned off to save computing time real :: dt_slow integer :: diag_freq = 1 ! Number of slow time steps between sending out diagnositics data. logical :: debug_river = .FALSE. logical :: do_age = .false. real :: Somin = 0.00005 ! There are 7 points with So = -9.999 but basinid > 0.... real :: outflowmean_min = 1. ! temporary fix, should not allow zero in input file integer :: num_c=0, num_species logical :: land_area_called_cellarea = .false. logical :: all_big_outlet_ctn0 = .false. !Balaji public :: num_species !public for test_river_solo character(len=6), dimension(10) :: rt_c_name character(len=128), dimension(10) :: rt_source_conc_file, rt_source_flux_file character(len=128), dimension(10) :: rt_source_conc_name, rt_source_flux_name real, dimension(10) :: rt_t_ref, rt_vf_ref, rt_q10, rt_kinv character(len=128),allocatable, dimension(:) :: source_conc_file, source_flux_file character(len=128),allocatable, dimension(:) :: source_conc_name, source_flux_name real, dimension(3) :: ave_DHG_exp = (/0.49,0.33,0.18/) ! (/B, F, M for avg of many rivers, 15Nov05/) real, dimension(3) :: ave_AAS_exp = (/0.19,0.39,0.42/) ! (/b, f, m for avg of many rivers, 15Nov05/) real, dimension(3) :: ave_DHG_coef = (/4.62,0.26,0.82/) ! (/A, C, K for avg of many rivers, 15Nov05/) real :: sinuosity = 1.3 real :: channel_tau = 86400*365.25*10 ! channel geometry reflects average flow over O(10 y) logical :: lake_area_bug = .FALSE. ! if set to true, reverts to buggy (quebec) ! behavior, where by mistake cell area was used instead of land area to ! compute the area of lakes. logical :: stop_on_mask_mismatch = .TRUE. ! if set to false, then the data mismatches (mmismatch ! of land and river masks, and discharges in pouints where there is no ocean) are reported, ! but don't cause the abort of the program. namelist /river_nml/ dt_slow, diag_freq, debug_river, do_age, & Somin, outflowmean_min, num_c, rt_c_name, rt_t_ref, & rt_vf_ref, rt_q10, rt_kinv, rt_source_conc_file, & rt_source_flux_file, rt_source_conc_name, & rt_source_flux_name, ave_DHG_exp, ave_AAS_exp, & ave_DHG_coef, do_rivers, sinuosity, channel_tau, & land_area_called_cellarea, all_big_outlet_ctn0, & lake_area_bug, stop_on_mask_mismatch character(len=128) :: river_src_file = 'INPUT/river_data.nc' character(len=128) :: river_Omean_file = 'INPUT/river_Omean.nc' !--------------------------------------------------------------------- logical :: module_is_initialized = .FALSE. integer :: isc, iec, jsc, jec ! compute domain decomposition integer :: isd, ied, jsd, jed ! data domain decomposition integer :: nlon, nlat ! size of computational river grid integer :: num_lake_lev integer :: id_outflowmean, id_lake_depth_sill integer :: id_dx, id_basin, id_So, id_depth, id_width, id_vel integer :: id_LWSr, id_FWSr, id_HSr, id_meltr integer :: i_species integer :: id_travel, id_elev, id_tocell ! *** ! the following id_* retained temporarily for compatibility with older diag tables integer :: id_storage_old, id_stordis_old, id_infloc_old, id_inflow_old, id_outflow_old integer :: id_lake_outflow_old, id_disw2o_old, id_outflowmean_old integer :: id_depth_old, id_width_old, id_vel_old integer :: id_r_t_rivr, id_i_t_rivr, id_o_t_rivr, id_lot_rivr, id_s_t_rivr, id_dot_rivr integer :: id_r_frazil, id_i_frazil, id_o_frazil, id_lofrazil, id_s_frazil, id_dofrazil ! *** integer :: maxtravel real :: missing = -1.e8 real, parameter :: CONST_OMEAN = 80000 integer, parameter :: num_phys = 2 real, parameter :: epsln = 1.e-6 real, parameter :: sec_in_day = 86400. real, allocatable, dimension(:,:) :: discharge2ocean_next ! store discharge value real, allocatable, dimension(:,:,:) :: discharge2ocean_next_c ! store discharge value integer, allocatable, dimension(:) :: id_infloc, id_storage, id_stordis, id_inflow integer, allocatable, dimension(:) :: id_run_stor integer, allocatable, dimension(:) :: id_outflow, id_removal, id_dis integer, allocatable, dimension(:) :: id_lake_outflow character(len=4), allocatable, dimension(:) :: c_name character(len=8), allocatable, dimension(:) :: if_name, of_name, lo_name, do_name character(len=8), allocatable, dimension(:) :: st_name, sd_name, rf_name, rm_name, sr_name character(len=8), allocatable, dimension(:) :: c_desc character(len=24), allocatable, dimension(:) :: if_desc, of_desc, lo_desc, do_desc character(len=64), allocatable, dimension(:) :: st_desc, sd_desc, rf_desc, rm_desc, sr_desc character(len=7), allocatable, dimension(:) :: flux_units, store_units character(len=5), allocatable, dimension(:) :: conc_units integer :: num_fast_calls integer :: slow_step = 0 ! record number of slow time step run. type(domain2d), save :: domain type(river_type) , save :: River !--- these variables are for communication purpose integer :: pe integer, allocatable :: ncells(:) ! number of points with each travel value !--- clock id variable integer :: slowclock, bndslowclock, physicsclock, diagclock, riverclock contains !##################################################################### subroutine river_init( land_lon, land_lat, time, dt_fast, land_domain, & land_frac, id_lon, id_lat, river_land_mask ) real, intent(in) :: land_lon(:,:) ! geographical lontitude of cell center real, intent(in) :: land_lat(:,:) ! geographical lattitude of cell center type(time_type), intent(in) :: time ! current time type(time_type), intent(in) :: dt_fast ! fast time step type(domain2d), intent(in) :: land_domain ! land domain real, intent(in) :: land_frac(:,:) ! land area fraction from land model integer, intent(in) :: id_lon, id_lat ! IDs of diagnostic axes logical, intent(out):: river_land_mask(:,:) ! land mask seen by rivers integer :: unit, io_status, ierr integer :: sec, day, i, j integer :: nxc, nyc character(len=128) :: filename type(Leo_Mad_trios) :: DHG_exp ! downstream equation exponents type(Leo_Mad_trios) :: DHG_coef ! downstream equation coefficients type(Leo_Mad_trios) :: AAS_exp ! at-a-station equation exponents riverclock = mpp_clock_id('update_river' , CLOCK_FLAG_DEFAULT, CLOCK_SUBCOMPONENT) slowclock = mpp_clock_id('update_river_slow' , CLOCK_FLAG_DEFAULT, CLOCK_ROUTINE) bndslowclock = mpp_clock_id('update_river_bnd_slow', CLOCK_FLAG_DEFAULT, CLOCK_ROUTINE) physicsclock = mpp_clock_id('river phys' , CLOCK_FLAG_DEFAULT, CLOCK_ROUTINE) diagclock = mpp_clock_id('river diag' , CLOCK_FLAG_DEFAULT, CLOCK_ROUTINE) !--- read namelist ------------------------------------------------- #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=river_nml, iostat=io_status) ierr = check_nml_error(io_status, 'river_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file() ierr = 1; do while (ierr /= 0) read (unit, nml=river_nml, iostat=io_status, end=10) ierr = check_nml_error(io_status,'river_nml') enddo 10 continue call close_file (unit) endif #endif !--- write version and namelist info to logfile -------------------- call write_version_number(version, tagname) unit=stdlog() write(unit, river_nml) if(.not.do_rivers) return ! do nothing further if the rivers are turned off !--- check name list variables if(diag_freq .le. 0) call mpp_error(FATAL,'river_mod: diag_freq should be a positive integer') pe = mpp_pe() ! set up time-related values River % time = time call get_time(dt_fast, sec, day) River%dt_fast = day*sec_in_day+sec River%dt_slow = dt_slow River%channel_tau = channel_tau num_fast_calls = River%dt_slow/River%dt_fast num_species = num_phys + num_c if (do_age) num_species = num_species + 1 River%num_species = num_species River%num_c = num_c River%do_age = do_age River%num_phys = num_phys if(River%dt_slow .lt. River%dt_fast) call mpp_error(FATAL, & 'river_mod: river slow time step dt_slow should be no less than land model fast time step dt_fast') if ( mod(River%dt_slow,River%dt_fast) .ne. 0 ) call mpp_error(FATAL, & 'river_mod: river slow time step dt_slow should be multiple of land model fast time step dt_fast') !--- get the domain decompsition, river and land will be on the same grid and have the same domain decomposition. domain = land_domain call mpp_get_global_domain (domain, xsize=River%nlon, ysize=River%nlat) call mpp_get_compute_domain(domain, isc, iec, jsc, jec) call mpp_get_data_domain (domain, isd, ied, jsd, jed) !---- make sure the halo size is 1 if( ied-iec .NE. 1 .OR. isc-isd .NE. 1 .OR. jed-jec .NE. 1 .OR. jsc-jsd .NE. 1 ) & call mpp_error(FATAL, "river_mod: halo size in four direction should all be 1") nxc = iec - isc + 1; nyc = jec - jsc + 1 !--- make sure land_lon, land_lat, land_frac is on the compute domain if(size(land_lon,1) .NE. nxc .OR. size(land_lon,2) .NE. nyc ) call mpp_error(FATAL, & "river_mod: land_lon should be on the compute domain") if(size(land_lat,1) .NE. nxc .OR. size(land_lat,2) .NE. nyc ) call mpp_error(FATAL, & "river_mod: land_lat should be on the compute domain") if(size(land_frac,1) .NE. nxc .OR. size(land_frac,2) .NE. nyc ) call mpp_error(FATAL, & "river_mod: land_frac should be on the compute domain") allocate(discharge2ocean_next (isc:iec,jsc:jec )) allocate(discharge2ocean_next_c(isc:iec,jsc:jec,num_species)) allocate(id_infloc (0:num_species), id_storage(0:num_species)) allocate(id_inflow (0:num_species), id_outflow(0:num_species)) allocate(id_dis (0:num_species), id_lake_outflow (0:num_species)) allocate(id_removal(0:num_species), id_stordis(0:num_species), id_run_stor(0:num_species)) allocate(c_name (0:num_species)) allocate(if_name (0:num_species), of_name (0:num_species)) allocate(rf_name (0:num_species), rm_name (0:num_species)) allocate(st_name (0:num_species), sd_name (0:num_species), sr_name(0:num_species)) allocate(lo_name (0:num_species), do_name (0:num_species)) allocate(c_desc (0:num_species)) allocate(if_desc (0:num_species), of_desc (0:num_species)) allocate(rf_desc (0:num_species), rm_desc (0:num_species)) allocate(st_desc (0:num_species), sd_desc (0:num_species), sr_desc(0:num_species)) allocate(lo_desc (0:num_species), do_desc (0:num_species)) allocate(store_units (0:num_species), flux_units (0:num_species)) allocate(conc_units (num_species)) allocate(source_conc_file(num_species-num_c+1:num_species)) allocate(source_flux_file(num_species-num_c+1:num_species)) allocate(source_conc_name(num_species-num_c+1:num_species)) allocate(source_flux_name(num_species-num_c+1:num_species)) !--- read the data from the file river_src_file -- has all static river network data call get_river_data(land_lon, land_lat, land_frac) river_land_mask = River%mask River%t_ref = rt_t_ref (1:num_c) River%vf_ref = rt_vf_ref(1:num_c) River%q10 = rt_q10 (1:num_c) River%kinv = rt_kinv (1:num_c) source_conc_file(num_species-num_c+1:num_species) = rt_source_conc_file(1:num_c) source_flux_file(num_species-num_c+1:num_species) = rt_source_flux_file(1:num_c) source_conc_name(num_species-num_c+1:num_species) = rt_source_conc_name(1:num_c) source_flux_name(num_species-num_c+1:num_species) = rt_source_flux_name(1:num_c) c_name(0)='_h2o' c_desc(0)='h2o mass' c_name(1)='_ice' c_desc(1)='ice mass' c_name(2)='_hEt' c_desc(2)='sen.heat' if (do_age) c_name(3)='_age' if (do_age) c_desc(3)='mass.age' do i_species = num_species-num_c+1, num_species c_name(i_species)=trim(rt_c_name(i_species-(num_species-num_c))) c_desc(i_species)=c_name(i_species) enddo conc_units(1) = 'kg/kg' conc_units(2) = ' K ' if (do_age) conc_units(3) = 'days ' conc_units(num_species-num_c+1:num_species) = 'kg/kg' ! check this one flux_units(0) = 'kg/m2/s' flux_units(1) = 'kg/m2/s' flux_units(2) = 'W/m2 ' if (do_age) flux_units(3) = 'kg/m2 ' flux_units(num_species-num_c+1:num_species) = 'kg/m2/s' store_units(0) = 'kg/m2 ' store_units(1) = 'kg/m2 ' store_units(2) = 'J/m2 ' if (do_age) store_units(3) = 'kg-s/m2' store_units(num_species-num_c+1:num_species)= 'kg/m2 ' do i_species = 0, num_species if_name(i_species)='rv_i'//trim(c_name(i_species)) of_name(i_species)='rv_o'//trim(c_name(i_species)) do_name(i_species)='rv_d'//trim(c_name(i_species)) lo_name(i_species)='rv_l'//trim(c_name(i_species)) rf_name(i_species)='rv_r'//trim(c_name(i_species)) rm_name(i_species)='rv_m'//trim(c_name(i_species)) st_name(i_species)='rv_s'//trim(c_name(i_species)) sd_name(i_species)='rv_n'//trim(c_name(i_species)) sr_name(i_species)='rv_u'//trim(c_name(i_species)) enddo do i_species = 0, num_species if_desc(i_species)='river inflow, ' //trim(c_desc(i_species)) of_desc(i_species)='river outflow, ' //trim(c_desc(i_species)) do_desc(i_species)='ocean discharge,' //trim(c_desc(i_species)) lo_desc(i_species)='lake outflow, ' //trim(c_desc(i_species)) rf_desc(i_species)='local runoff, ' //trim(c_desc(i_species)) rm_desc(i_species)='river removal, ' //trim(c_desc(i_species)) st_desc(i_species)='river storage ' //trim(c_desc(i_species)) sd_desc(i_species)='river discharge lag (numerical) storage, ' & //trim(c_desc(i_species)) sr_desc(i_species)='river runoff lag (numerical) storage, ' & //trim(c_desc(i_species)) enddo !--- register diag field call river_diag_init (id_lon, id_lat) !--- read restart file call get_instance_filename('INPUT/river.res.nc', filename) call get_mosaic_tile_file(trim(filename), filename, .false., domain) if(file_exist(trim(filename),domain) ) then call mpp_error(NOTE, 'river_init : Read restart files '//trim(filename)) call read_data(filename,'storage', River%storage, domain) call read_data(filename,'storage_c', River%storage_c, domain) call read_data(filename,'discharge2ocean', discharge2ocean_next, domain) call read_data(filename,'discharge2ocean_c',discharge2ocean_next_c, domain) call read_data(filename,'Omean', River%outflowmean, domain) if (field_exist(filename,'depth',domain)) then ! call mpp_error(WARNING, 'river_init : Reading field "depth" from '//trim(filename)) call read_data(filename,'depth', River%depth, domain) else ! call mpp_error(WARNING, 'river_init : "depth" is not present in '//trim(filename)) endif else call mpp_error(NOTE, 'river_init : cold start, set data to 0') River%storage = 0.0 River%storage_c = 0.0 discharge2ocean_next = 0.0 discharge2ocean_next_c = 0.0 if(file_exist(river_Omean_file)) then call read_data(river_Omean_file, 'Omean', River%outflowmean, domain) else River%outflowmean = CONST_OMEAN end if endif River%stordis_c = River%dt_slow * discharge2ocean_next_c/DENS_H2O River%stordis = River%dt_slow *(discharge2ocean_next + & discharge2ocean_next_c(:,:,1))/DENS_H2O where(River%outflowmean .le. outflowmean_min) River%outflowmean=outflowmean_min maxtravel = maxval(River%travel) call mpp_max(maxtravel) call river_physics_init(River, domain, id_lon, id_lat) call get_Leo_Mad_params(DHG_exp, DHG_coef, AAS_exp) River%o_exp = 1./ (AAS_exp%on_w + AAS_exp%on_d) do j = jsc, jec do i = isc, iec if ( River%reach_length(i,j) > 0.0) then River%o_coef(i,j) = River%outflowmean(i,j) / & ((sinuosity*River%reach_length(i,j))*DHG_coef%on_w*DHG_coef%on_d & *(River%outflowmean(i,j)**(DHG_exp%on_w+DHG_exp%on_d)))**River%o_exp endif enddo enddo River%d_exp = AAS_exp%on_d River%d_coef = DHG_coef%on_d & *(River%outflowmean**(DHG_exp%on_d-AAS_exp%on_d)) River%w_exp = AAS_exp%on_w River%w_coef = DHG_coef%on_w & *(River%outflowmean**(DHG_exp%on_w-AAS_exp%on_w)) num_lake_lev = num_l module_is_initialized = .TRUE. end subroutine river_init !##################################################################### subroutine update_river ( runoff, runoff_c, discharge2ocean, & discharge2ocean_c ) real, dimension(:,:), intent(in) :: runoff real, dimension(:,:,:), intent(in) :: runoff_c real, dimension(:,:), intent(out) :: discharge2ocean real, dimension(:,:,:), intent(out) :: discharge2ocean_c integer, save :: n = 0 ! fast time step with each slow time step call mpp_clock_begin(riverclock) if (.not.do_rivers) then discharge2ocean = 0; discharge2ocean_c = 0 call mpp_clock_end(riverclock) !needed for early exit when do_rivers=.false. return endif discharge2ocean = discharge2ocean_next discharge2ocean_c = discharge2ocean_next_c ! deplete the discharge storage pools River%stordis_c = River%stordis_c & - River%dt_fast * discharge2ocean_c/DENS_H2O River%stordis = River%stordis & - River%dt_fast *(discharge2ocean + & discharge2ocean_c(:,:,1))/DENS_H2O ! increment time River%Time = increment_time(River%Time, River%dt_fast, 0) n = n + 1 !--- accumulate runoff --------------------- River%run_stor = River%run_stor + runoff River%run_stor_c = River%run_stor_c + runoff_c if(n == num_fast_calls) then call mpp_clock_begin(slowclock) call update_river_slow(River%run_stor(:,:)/real(num_fast_calls), & River%run_stor_c(:,:,:)/real(num_fast_calls) ) call mpp_clock_end(slowclock) call mpp_clock_begin(bndslowclock) call update_river_bnd_slow call mpp_clock_end(bndslowclock) n = 0 River%run_stor = 0 River%run_stor_c = 0 endif call mpp_clock_end(riverclock) end subroutine update_river !##################################################################### subroutine update_river_slow(runoff, runoff_c) real, dimension(:,:), intent(in) :: runoff real, dimension(:,:,:), intent(in) :: runoff_c real, dimension(isd:ied,jsd:jed) :: & lake_sfc_A, lake_sfc_bot, lake_conn real, dimension(isd:ied,jsd:jed,num_lake_lev) :: & lake_wl, lake_ws real, dimension(isc:iec,jsc:jec) :: & lake_depth_sill, lake_width_sill, lake_backwater, & lake_backwater_1, & lake_whole_area, & rivr_LMASS, & ! mass of liquid water in rivers in cell rivr_FMASS, & ! mass of ice in rivers in cell rivr_MELT, & ! net mass melt in rivers in cell rivr_HEAT ! sensible heat content of rivers in cell real, dimension(isc:iec,jsc:jec,num_lake_lev) :: & lake_T integer :: travelnow, lev type(Leo_Mad_trios) :: DHG_exp type(Leo_Mad_trios) :: DHG_coef type(Leo_Mad_trios) :: AAS_exp integer i,j,k, i_next, j_next type(land_tile_enum_type) :: te,ce ! last and current tile list elements type(land_tile_type), pointer :: tile ! pointer to current tile logical :: used slow_step = slow_step + 1 River%infloc = River%land_area*runoff /DENS_H2O River%infloc_c = 0 do i_species = 1, River%num_phys River%infloc_c(:,:,i_species) = & River%land_area*runoff_c(:,:,i_species)/DENS_H2O enddo if (River%do_age) then i_species = 3 River%infloc_c(:,:,i_species) = & River%land_area*runoff_c(:,:,i_species)/DENS_H2O endif do i_species = num_species-River%num_c+1, num_species ! create mass flux inputs from c data where (River%land_area.gt.0.) & River%infloc_c(:,:,i_species) = & River%infloc*River%source_conc(:,:,i_species) & + River%source_flux(:,:,i_species) enddo River%inflow = 0 River%inflow_c = 0 River%lake_outflow = 0 River%lake_outflow_c = 0 River%disw2o = 0. River%disc2o = 0. River%melt = 0. lake_sfc_A = 0 lake_sfc_bot= 0 lake_T = 0 lake_wl = 0 lake_ws = 0 lake_depth_sill = 0 lake_width_sill = 0 lake_whole_area = 0 lake_conn = 0 lake_backwater = 0 lake_backwater_1 = 0 ce = first_elmt(lnd%tile_map, is=isc, js=jsc) te = tail_elmt (lnd%tile_map) do while(ce /= te) call get_elmt_indices(ce,i,j,k) tile=>current_tile(ce) ! get pointer to current tile ce=next_elmt(ce) ! advance position to the next tile if (.not.associated(tile%lake)) cycle if (lake_area_bug) then lake_sfc_A (i,j) = tile%frac * lnd%cellarea(i,j) else lake_sfc_A (i,j) = tile%frac * lnd%area(i,j) endif do lev = 1, num_lake_lev lake_T (i,j,lev) = tile%lake%prog(lev)%T lake_wl(i,j,lev) = tile%lake%prog(lev)%wl lake_ws(i,j,lev) = tile%lake%prog(lev)%ws enddo lake_sfc_bot(i,j) = (sum(tile%lake%prog(:)%wl+tile%lake%prog(:)%ws) & -tile%lake%prog(1)%wl-tile%lake%prog(1)%ws ) & / DENS_H2O lake_depth_sill(i,j) = tile%lake%pars%depth_sill lake_width_sill(i,j) = tile%lake%pars%width_sill lake_whole_area(i,j) = tile%lake%pars%whole_area lake_conn (i,j) = tile%lake%pars%connected_to_next lake_backwater(i,j) = tile%lake%pars%backwater lake_backwater_1(i,j) = tile%lake%pars%backwater_1 enddo call mpp_update_domains (lake_sfc_A, domain) call mpp_update_domains (lake_sfc_bot,domain) call mpp_update_domains (lake_wl, domain) call mpp_update_domains (lake_ws, domain) call mpp_update_domains (lake_conn, domain) do i=isc,iec do j=jsc,jec if (River%i_tocell(i,j) /= NO_RIVER_FLAG) then i_next = River%i_tocell(i,j) j_next = River%j_tocell(i,j) else ! to avoid indices out of bounds in the lake_sfc_A check i_next = i; j_next=j endif if (lake_backwater(i,j).gt.0.5 .and. lake_sfc_A(i,j).gt.0. .and. & lake_sfc_A(i_next,j_next).gt.0. ) then ! because of river backwater, lake in this cell relaxes toward level of ! lake in next cell downstream. (river depth is still simple function ! of discharge though.) lake_depth_sill(i,j) = lake_sfc_bot(i_next,j_next) & +(lake_wl(i_next,j_next,1)+lake_ws(i_next,j_next,1))/DENS_H2O elseif (lake_backwater_1(i,j).gt.0.5) then ! to determine depth of backwater, lake at coastal cell has base level ! set to river depth in same cell lake_depth_sill(i,j) = lake_depth_sill(i,j) + River%depth(i,j) elseif (lake_conn(i,j).gt.0.5 ) then ! for all but furthest dowstream cell of a multi-cell lake, ! relax toward level in next cell (same lake) downstream if (lake_conn(i_next,j_next).gt.0.5 .or. all_big_outlet_ctn0) then lake_depth_sill(i,j) = lake_sfc_bot(i_next,j_next) & +(lake_wl(i_next,j_next,1)+lake_ws(i_next,j_next,1))/DENS_H2O endif elseif (river_impedes_lake) then if (lake_width_sill(i,j).lt.0..or.river_impedes_large_lake) then ! lake level in cell relaxes toward river level in cell lake_depth_sill(i,j) = lake_depth_sill(i,j) + River%depth(i,j) endif endif enddo enddo ! leftovers from horizontal mixing option, now gone !call mpp_update_domains (lake_T, domain) !call mpp_update_domains (lake_depth_sill, domain) !call mpp_update_domains (lake_tau, domain) travelnow = maxtravel do travelnow = maxtravel, 0, -1 call mpp_clock_begin(physicsclock) !*************************************************************** call river_physics_step (River, travelnow, & lake_sfc_A, lake_sfc_bot, lake_depth_sill, & lake_width_sill, lake_whole_area, & lake_T, lake_wl, lake_ws ) !*************************************************************** call mpp_clock_end(physicsclock) enddo ce = first_elmt(lnd%tile_map, is=isc, js=jsc) te = tail_elmt (lnd%tile_map) do while(ce /= te) call get_elmt_indices(ce,i,j,k) tile=>current_tile(ce) ! get pointer to current tile ce=next_elmt(ce) ! advance position to the next tile if (.not.associated(tile%lake)) cycle do lev = 1, num_lake_lev tile%lake%prog(lev)%T = lake_T (i,j,lev) tile%lake%prog(lev)%wl = lake_wl(i,j,lev) tile%lake%prog(lev)%ws = lake_ws(i,j,lev) enddo enddo River%outflowmean = River%outflowmean + & (River%outflow-River%outflowmean)*River%dt_slow/River%channel_tau where(River%outflowmean .le. outflowmean_min) River%outflowmean=outflowmean_min call get_Leo_Mad_params(DHG_exp, DHG_coef, AAS_exp) do j = jsc, jec do i = isc, iec if ( River%reach_length(i,j) > 0.0) then River%o_coef(i,j) = River%outflowmean(i,j) / & ((sinuosity*River%reach_length(i,j))*DHG_coef%on_w*DHG_coef%on_d & *(River%outflowmean(i,j)**(DHG_exp%on_w+DHG_exp%on_d)))**River%o_exp endif enddo enddo River%d_coef = DHG_coef%on_d & *(River%outflowmean**(DHG_exp%on_d-AAS_exp%on_d)) River%w_coef = DHG_coef%on_w & *(River%outflowmean**(DHG_exp%on_w-AAS_exp%on_w)) River%stordis = River%dt_slow*River%disw2o do i_species = 1, num_species River%stordis_c(:,:,i_species) = River%dt_slow*River%disc2o(:,:,i_species) enddo rivr_FMASS = DENS_H2O * (River%storage_c(:,:,1) + River%stordis_c(:,:,1)) rivr_LMASS = DENS_H2O * (River%storage + River%stordis) - rivr_FMASS rivr_MELT = DENS_H2O * River%melt / River%dt_slow rivr_HEAT = DENS_H2O * (River%storage_c(:,:,2) + River%stordis_c(:,:,2)) & - hlf*rivr_FMASS call mpp_clock_begin(diagclock) ! convert area-integrated river outputs to unit-area quantities, ! using land area for stores, cell area for fluxes to ocean if (id_LWSr > 0) then where (lnd%area > 0) & rivr_LMASS = rivr_LMASS / lnd%area used = send_data (id_LWSr, rivr_LMASS, River%time, mask=lnd%area>0) endif if (id_FWSr > 0) then where (lnd%area > 0) & rivr_FMASS = rivr_FMASS / lnd%area used = send_data (id_FWSr, rivr_FMASS, River%time, mask=lnd%area>0) endif if (id_HSr > 0) then where (lnd%area > 0) & rivr_HEAT = rivr_HEAT / lnd%area used = send_data (id_HSr, rivr_HEAT, River%time, mask=lnd%area>0) endif if (id_meltr > 0) then where (lnd%area > 0) & rivr_MELT = rivr_MELT / lnd%area used = send_data (id_meltr, rivr_MELT, River%time, mask=lnd%area>0) end if if(mod(slow_step, diag_freq) == 0) call river_diag(lake_depth_sill) call mpp_clock_end(diagclock) end subroutine update_river_slow !##################################################################### subroutine update_river_bnd_slow ! note that land_area is not the total area of the cell, but just the land area ! within the cell, so it cannot be used to normalize fluxes to all-ocean cells. ! we need a true cell area for normalization, so river will ! just return mass flux per unit time and let land_model divide by area discharge2ocean_next = DENS_H2O*(River%disw2o - River%disc2o(:,:,1)) do i_species = 1, num_species discharge2ocean_next_c(:,:,i_species) = DENS_H2O*River%disc2o(:,:,i_species) enddo end subroutine update_river_bnd_slow !##################################################################### subroutine river_end integer :: outunit ! unit number for stdout if(.not.do_rivers) return ! do nothing further if rivers are turned off !--- write out checksum outunit=stdout() write(outunit,*)"Chksum for storage ==> ", mpp_chksum(River%storage(isc:iec,jsc:jec)) write(outunit,*)"Chksum for storage_c ==> ", mpp_chksum(River%storage_c(isc:iec,jsc:jec,:)) write(outunit,*)"Chksum for discharge2ocean_next ==> ", mpp_chksum(discharge2ocean_next(isc:iec,jsc:jec)) write(outunit,*)"Chksum for discharge2ocean_next_c ==> ", mpp_chksum(discharge2ocean_next_c(isc:iec,jsc:jec,:)) !--- release memory deallocate(discharge2ocean_next, discharge2ocean_next_c,& River%run_stor, River%run_stor_c) deallocate( River%lon, River%lat) deallocate(River%land_area , River%basinid ) deallocate(River%landfrac ) deallocate(River%tocell ) deallocate(River%travel ) deallocate(River%outflow ) deallocate(River%inflow ) deallocate(River%lake_outflow) deallocate(River%lake_outflow_c) deallocate(River%storage ) deallocate(River%stordis ) deallocate(River%melt ) deallocate(River%disw2o ) deallocate(River%disc2o ) deallocate(River%infloc ) deallocate(River%reach_length ) deallocate(River%mask ) deallocate(River%So ) deallocate(River%depth ) deallocate(River%width ) deallocate(River%vel ) deallocate(River%infloc_c , River%storage_c , River%stordis_c ) deallocate(River%inflow_c, River%outflow_c ) deallocate(River%removal_c ) deallocate(River%vf_ref,River%t_ref,River%q10,River%kinv) deallocate(River%d_coef,River%o_coef,River%w_coef) module_is_initialized = .FALSE. end subroutine river_end !##################################################################### !--- write to restart file subroutine save_river_restart(timestamp) character(*), intent(in) :: timestamp character(len=128) :: filename if(.not.do_rivers) return ! do nothing further if rivers are turned off filename = 'RESTART/'//trim(timestamp)//'river.res.nc' call write_data(filename,'storage', River%storage(isc:iec,jsc:jec), domain) call write_data(filename,'storage_c', River%storage_c(isc:iec,jsc:jec,:), domain) !--- write out discharge data call write_data(filename,'discharge2ocean' ,discharge2ocean_next (isc:iec,jsc:jec), domain) call write_data(filename,'discharge2ocean_c',discharge2ocean_next_c(isc:iec,jsc:jec,:), domain) call write_data(filename,'Omean', River%outflowmean, domain) call write_data(filename,'depth', River%depth, domain) end subroutine save_river_restart !##################################################################### subroutine get_river_data(land_lon, land_lat, land_frac) real, intent(in) :: land_lon(isc:,jsc:) ! geographical lontitude of cell center real, intent(in) :: land_lat(isc:,jsc:) ! geographical lattitude of cell center real, intent(in) :: land_frac(isc:,jsc:) ! land area fraction of land grid. integer :: ni, nj, i, j, siz(4), ntiles real, dimension(:,:), allocatable :: xt, yt, frac, glon, glat, lake_frac integer :: nerrors ! number of errors detected during initialization ntiles = mpp_get_ntile_count(domain) call field_size(river_src_file, 'basin', siz, domain=domain) ni = siz(1) nj = siz(2) if(ni .NE. River%nlon .OR. nj .NE. River%nlat) call mpp_error(FATAL, & "river_mod: size mismatch between river grid and land grid") allocate(glon(ni,nj), glat(ni, nj)) allocate(xt(isc:iec, jsc:jec), yt(isc:iec, jsc:jec), frac(isc:iec, jsc:jec) ) allocate(lake_frac(isc:iec, jsc:jec)) if (ntiles == 1) then call read_data(river_src_file, 'x', glon, no_domain=.true.) call read_data(river_src_file, 'y', glat, no_domain=.true.) endif call read_data(river_src_file, 'x', xt, domain) call read_data(river_src_file, 'y', yt, domain) call read_data(river_src_file, 'land_frac', frac, domain) call read_data(river_src_file, 'lake_frac', lake_frac, domain) !--- the following will be changed when the river data sets is finalized. xt = land_lon yt = land_lat !--- transform to radians, since land model grid use radians and compare with land grid. allocate(River%lon_1d (1:ni ) ) allocate(River%lat_1d (1:nj ) ) allocate(River%lon (isc:iec, jsc:jec) ) allocate(River%lat (isc:iec, jsc:jec) ) allocate(River%land_area (isc:iec, jsc:jec) ) allocate(River%basinid (isc:iec, jsc:jec) ) allocate(River%landfrac (isc:iec, jsc:jec) ) allocate(River%mask (isc:iec, jsc:jec) ) allocate(River%tocell (isc:iec, jsc:jec) ) allocate(River%i_tocell (isc:iec, jsc:jec) ) allocate(River%j_tocell (isc:iec, jsc:jec) ) allocate(River%travel (isd:ied, jsd:jed) ) allocate(River%inflow (isc:iec, jsc:jec) ) allocate(River%outflow (isc:iec, jsc:jec) ) allocate(River%lake_outflow(isc:iec, jsc:jec) ) allocate(River%storage (isc:iec, jsc:jec) ) allocate(River%stordis (isc:iec, jsc:jec) ) allocate(River%run_stor (isc:iec, jsc:jec) ) allocate(River%melt (isc:iec, jsc:jec) ) allocate(River%disw2o (isc:iec, jsc:jec) ) allocate(River%infloc (isc:iec, jsc:jec)) allocate(River%reach_length(isc:iec, jsc:jec) ) allocate(River%So (isc:iec, jsc:jec) ) allocate(River%depth (isc:iec, jsc:jec) ) allocate(River%width (isc:iec, jsc:jec) ) allocate(River%vel (isc:iec, jsc:jec) ) allocate(River%infloc_c (isc:iec, jsc:jec, num_species) ) allocate(River%storage_c (isc:iec, jsc:jec, num_species) ) allocate(River%stordis_c (isc:iec, jsc:jec, num_species) ) allocate(River%run_stor_c (isc:iec, jsc:jec, num_species) ) allocate(River%outflow_c (isc:iec, jsc:jec, num_species) ) allocate(River%lake_outflow_c (isc:iec, jsc:jec, num_species) ) allocate(River%removal_c (isc:iec, jsc:jec, num_species) ) allocate(River%inflow_c (isc:iec, jsc:jec, num_species) ) allocate(River%disc2o (isc:iec, jsc:jec, num_species)) allocate(River%d_coef (isc:iec, jsc:jec) ) allocate(River%o_coef (isc:iec, jsc:jec) ) allocate(River%w_coef (isc:iec, jsc:jec) ) allocate(River%outflowmean(isc:iec, jsc:jec) ) allocate(River%t_ref(4:num_species),River%vf_ref(4:num_species)) allocate(River%q10 (4:num_species),River%kinv (4:num_species)) allocate(River%source_conc(isc:iec, jsc:jec,num_species-num_c+1:num_species)) allocate(River%source_flux(isc:iec, jsc:jec,num_species-num_c+1:num_species)) if(ntiles == 1) then ! lat-lon grid, use actual grid location River%lon_1d(:) = glon(:,1) River%lat_1d(:) = glat(1,:) else ! cubic grid, use index. River%lon_1d(:) = (/ (i, i=1,River%nlon) /) River%lat_1d(:) = (/ (i, i=1,River%nlat) /) end if deallocate(glon, glat) River%lon(:,:) = land_lon(:,:) River%lat(:,:) = land_lat(:,:) !!$ River%landfrac(:,:) = land_frac(:,:) River%landfrac(:,:) = frac(:,:) River%infloc = 0.0 River%infloc_c = 0.0 River%storage = 0.0 River%storage_c = 0.0 River%stordis = 0.0 River%run_stor = 0.0 River%stordis_c = 0.0 River%run_stor_c= 0.0 River%removal_c = 0.0 River%depth = 0. River%width = 0. River%vel = 0. River%outflow = 0. River%outflow_c = 0. River%inflow = 0. River%inflow_c = 0. !--- read the data from the source file call read_data(river_src_file, 'tocell', River%tocell, domain) where (River%tocell(:,:).eq. 4) River%tocell(:,:)=3 where (River%tocell(:,:).eq. 8) River%tocell(:,:)=4 where (River%tocell(:,:).eq. 16) River%tocell(:,:)=5 where (River%tocell(:,:).eq. 32) River%tocell(:,:)=6 where (River%tocell(:,:).eq. 64) River%tocell(:,:)=7 where (River%tocell(:,:).eq.128) River%tocell(:,:)=8 nerrors = 0 do j = jsc, jec do i = isc, iec !!$ if(abs(xt(i,j) - land_lon(i,j)) > epsln) call mpp_error(FATAL, & !!$ "get_river_data: longitude mismatch between river grid and land grid") !!$ if(abs(yt(i,j) - land_lat(i,j)) > epsln) call mpp_error(FATAL, & !!$ "get_river_data: latitude mismatch between river grid and land grid") !!$ if(abs(frac(i,j) - land_frac(i,j)) > epsln) call mpp_error(FATAL, & !!$ "get_river_data: area fraction mismatch between river grid and land grid") ! check that river and land masks match if ((frac(i,j)>0).neqv.(land_frac(i,j)>0)) then call mpp_error(WARNING,'get_river_data: land and river masks do not match at '//& trim(coordinates(i,j))) nerrors = nerrors+1 endif ! check that the rivers do not discarge in the middle of the continents if ((River%tocell(i,j)==0).and.(land_frac(i,j)>1.0-epsln)) then call mpp_error(WARNING, & 'get_river_data: river discharges into a land point '& //trim(coordinates(i,j))//' where there is no ocean') nerrors = nerrors+1 endif end do end do if (nerrors>0.and.stop_on_mask_mismatch) call mpp_error(FATAL,& 'get_river_data: river/land mask-related mismatch detected during river data initialization') call read_data(river_src_file, 'basin', River%basinid, domain) where (River%basinid >0) River%mask = .true. elsewhere River%mask = .false. endwhere River%travel = 0 call read_data(river_src_file, 'travel', River%travel(isc:iec,jsc:jec), domain) call mpp_update_domains(River%travel, domain) call read_data(river_src_file, 'celllength', River%reach_length, domain) River%reach_length = River%reach_length * River%landfrac * (1.-lake_frac) if (land_area_called_cellarea) then call read_data(river_src_file, 'cellarea', River%land_area, domain) else call read_data(river_src_file, 'land_area', River%land_area, domain) endif ! call read_data(river_src_file, 'So', River%So, domain) River%So = 0.0 where (River%So .LT. 0.0) River%So = Somin do i_species = num_species-num_c+1, num_species if (trim(source_conc_file(i_species)).eq.'') then River%source_conc(:,:,i_species)=0 if (trim(source_conc_name(i_species)).eq.'one') River%source_conc(:,:,i_species)=1 else if (trim(source_conc_name(i_species)).ne.'') then call read_data(trim(source_conc_file(i_species)), trim(source_conc_name(i_species)), & River%source_conc(:,:,i_species), no_domain=.true.) else River%source_conc(:,:,i_species) = 0 endif if (trim(source_flux_file(i_species)).eq.'') then River%source_flux(:,:,i_species)=0 if (trim(source_flux_name(i_species)).eq.'one') River%source_flux(:,:,i_species)=1 else if (trim(source_flux_name(i_species)).ne.'') then call read_data(trim(source_flux_file(i_species)), & trim(source_flux_name(i_species)), & River%source_flux(:,:,i_species), no_domain=.true.) else River%source_flux(:,:,i_species) = 0 endif enddo River%source_conc = max(River%source_conc, 0.) River%source_flux = max(River%source_flux, 0.) deallocate(lake_frac) end subroutine get_river_data !##################################################################### subroutine river_diag_init(id_lon, id_lat) integer, intent(in) :: id_lon ! ID of land longitude (X) diag axis integer, intent(in) :: id_lat ! ID of land latitude (Y) diag axis character(len=11) :: mod_name = 'river' real, dimension(isc:iec,jsc:jec) :: tmp logical :: sent integer :: i ! regular diagnostic fields do i_species = 0, num_species id_inflow(i_species) = register_diag_field ( mod_name, if_name(i_species), & (/id_lon, id_lat/), River%Time, if_desc(i_species), flux_units(i_species), & missing_value=missing ) id_outflow(i_species) = register_diag_field ( mod_name, of_name(i_species), & (/id_lon, id_lat/), River%Time, of_desc(i_species), flux_units(i_species), & missing_value=missing ) id_dis(i_species) = register_diag_field ( mod_name, do_name(i_species), & (/id_lon, id_lat/), River%Time, do_desc(i_species), flux_units(i_species), & missing_value=missing ) id_lake_outflow(i_species) = register_diag_field ( mod_name, lo_name(i_species), & (/id_lon, id_lat/), River%Time, lo_desc(i_species), flux_units(i_species), & missing_value=missing ) id_infloc(i_species) = register_diag_field ( mod_name, rf_name(i_species), & (/id_lon, id_lat/), River%Time, rf_desc(i_species), flux_units(i_species), & missing_value=missing ) id_removal(i_species) = register_diag_field ( mod_name, rm_name(i_species), & (/id_lon, id_lat/), River%Time, rm_desc(i_species), flux_units(i_species), & missing_value=missing ) id_storage(i_species) = register_diag_field ( mod_name, st_name(i_species), & (/id_lon, id_lat/), River%Time, st_desc(i_species), store_units(i_species), & missing_value=missing ) id_stordis(i_species) = register_diag_field ( mod_name, sd_name(i_species), & (/id_lon, id_lat/), River%Time, sd_desc(i_species), store_units(i_species), & missing_value=missing ) id_run_stor(i_species) = register_diag_field ( mod_name, sr_name(i_species), & (/id_lon, id_lat/), River%Time, sr_desc(i_species), store_units(i_species), & missing_value=missing ) enddo id_lake_depth_sill= register_diag_field ( mod_name, 'rv_dsill', (/id_lon, id_lat/), & River%Time, 'effective lake sill depth', 'm', missing_value=missing ) id_outflowmean = register_diag_field ( mod_name, 'rv_Qavg', (/id_lon, id_lat/), & River%Time, 'long-time average vol. flow', 'm3/s', missing_value=missing ) id_depth = register_diag_field ( mod_name, 'rv_depth', (/id_lon, id_lat/), & River%Time, 'river flow depth', 'm', missing_value=missing ) id_width = register_diag_field ( mod_name, 'rv_width', (/id_lon, id_lat/), & River%Time, 'river flow width', 'm', missing_value=missing ) id_vel = register_diag_field ( mod_name, 'rv_veloc', (/id_lon, id_lat/), & River%Time, 'river flow velocity', 'm/s', missing_value=missing ) id_storage_old=register_diag_field(mod_name,'storage',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_s_h2o instead', store_units(0),missing_value=missing ) id_stordis_old=register_diag_field(mod_name,'stordis',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_s_h2o instead', store_units(0),missing_value=missing ) id_s_frazil=register_diag_field(mod_name,'s_frazil',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_s_ice instead', store_units(1),missing_value=missing ) id_s_t_rivr=register_diag_field(mod_name,'s_t_rivr',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_s_hEt instead', store_units(2),missing_value=missing ) id_infloc_old=register_diag_field(mod_name,'infloc',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_r_h2o instead', flux_units(0),missing_value=missing ) id_r_frazil=register_diag_field(mod_name,'r_frazil',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_r_ice instead', flux_units(1),missing_value=missing ) id_r_t_rivr=register_diag_field(mod_name,'r_t_rivr',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_r_hEt instead', flux_units(2),missing_value=missing ) id_inflow_old=register_diag_field(mod_name,'inflow',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_r_h2o instead', flux_units(0),missing_value=missing ) id_i_frazil=register_diag_field(mod_name,'i_frazil',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_r_ice instead', flux_units(1),missing_value=missing ) id_i_t_rivr=register_diag_field(mod_name,'i_t_rivr',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_r_hEt instead', flux_units(2),missing_value=missing ) id_outflow_old=register_diag_field(mod_name,'outflow',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_o_h2o instead', flux_units(0),missing_value=missing ) id_o_frazil=register_diag_field(mod_name,'o_frazil',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_o_ice instead', flux_units(1),missing_value=missing ) id_o_t_rivr=register_diag_field(mod_name,'o_t_rivr',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_o_hEt instead', flux_units(2),missing_value=missing ) id_lake_outflow_old=register_diag_field(mod_name,'lake_outflow',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_l_h2o instead', flux_units(0),missing_value=missing ) id_lofrazil=register_diag_field(mod_name,'lofrazil',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_l_ice instead', flux_units(1),missing_value=missing ) id_lot_rivr=register_diag_field(mod_name,'lot_rivr',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_l_hEt instead', flux_units(2),missing_value=missing ) id_disw2o_old=register_diag_field(mod_name,'disw2o',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_d_h2o instead', flux_units(0),missing_value=missing ) id_dofrazil=register_diag_field(mod_name,'dofrazil',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_d_ice instead', flux_units(1),missing_value=missing ) id_dot_rivr=register_diag_field(mod_name,'dot_rivr',(/id_lon, id_lat/), River%Time, & 'obsolete, pls use rv_d_hEt instead', flux_units(2),missing_value=missing ) id_outflowmean_old = register_diag_field ( mod_name, 'flowmean', (/id_lon, id_lat/), & River%Time, 'obsolete, pls use rv_Qavg', 'm3/s', missing_value=missing ) id_depth_old = register_diag_field ( mod_name, 'depth', (/id_lon, id_lat/), & River%Time, 'obsolete, pls use rv_depth', 'm', missing_value=missing ) id_width_old = register_diag_field ( mod_name, 'width', (/id_lon, id_lat/), & River%Time, 'obsolete, pls use rv_width', 'm', missing_value=missing ) id_vel_old = register_diag_field ( mod_name, 'vel', (/id_lon, id_lat/), & River%Time, 'obsolete, pls use rv_veloc', 'm/s', missing_value=missing ) ! fields that historically were in the the land_model.F90. They are registered ! for module 'land' to preserve compatibility with older diag tables id_LWSr = register_diag_field ( 'land', 'LWSr', (/id_lon, id_lat/), & River%Time, 'river liquid mass storage', 'kg/m2', missing_value=-1.0e+20 ) id_FWSr = register_diag_field ( 'land', 'FWSr', (/id_lon, id_lat/), & River%Time, 'river ice mass storage', 'kg/m2', missing_value=-1.0e+20 ) id_HSr = register_diag_field ( 'land', 'HSr', (/id_lon, id_lat/), & River%Time, 'river heat storage', 'J/m2', missing_value=-1.0e+20 ) id_meltr = register_diag_field ( 'land', 'meltr', (/id_lon, id_lat/), & River%Time, 'melt in river system', 'kg/m2/s', missing_value=-1.0e+20 ) ! static fields id_dx = register_static_field ( mod_name, 'rv_length', (/id_lon, id_lat/), & 'river reach length', 'm', missing_value=missing ) id_basin = register_static_field ( mod_name, 'rv_basin', (/id_lon, id_lat/), & 'river basin id', 'none', missing_value=missing ) id_So = register_static_field ( mod_name, 'So', (/id_lon, id_lat/), & 'Slope', 'none', missing_value=missing ) id_travel = register_static_field ( mod_name, 'rv_trav', (/id_lon, id_lat/), & 'cells left to travel before reaching ocean', 'none', missing_value=missing ) id_tocell = register_static_field ( mod_name, 'rv_dir', (/id_lon, id_lat/), & 'outflow direction code', 'none', missing_value=missing ) if (id_dx>0) sent=send_data(id_dx, River%reach_length, River%Time, mask=River%mask ) if (id_basin>0) then tmp = River%basinid(isc:iec,jsc:jec) sent=send_data(id_basin, tmp, River%Time, mask=River%mask ) end if if (id_So>0) sent=send_data(id_So, River%So, River%Time, mask=River%mask ) if (id_travel>0) then tmp = River%travel(isc:iec,jsc:jec) sent=send_data(id_travel, tmp, River%Time, mask=River%mask ) end if if (id_tocell>0) then tmp = River%tocell(isc:iec,jsc:jec) sent=send_data(id_tocell, tmp, River%Time, mask=River%mask ) end if end subroutine river_diag_init !##################################################################### subroutine river_diag(lake_depth_sill) real, dimension(isc:iec,jsc:jec), intent(in) :: lake_depth_sill logical :: used ! logical for send_data real diag_factor (isc:iec,jsc:jec) real diag_factor_2(isc:iec,jsc:jec) diag_factor = 0. diag_factor_2 = 0. where (River%land_area(isc:iec,jsc:jec).gt.0.) & diag_factor=DENS_H2O/River%land_area(isc:iec,jsc:jec) where (River%land_area(isc:iec,jsc:jec).gt.0.) & diag_factor_2=1./(River%land_area(isc:iec,jsc:jec)*River%dt_slow) if (id_inflow(0) > 0) used = send_data (id_inflow(0), & diag_factor*River%inflow(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_outflow(0) > 0) used = send_data (id_outflow(0), & diag_factor*River%outflow(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_storage(0) > 0) used = send_data (id_storage(0), & diag_factor*River%storage(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_stordis(0) > 0) used = send_data (id_stordis(0), & diag_factor*River%stordis(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_run_stor(0) > 0) used = send_data (id_run_stor(0), & River%dt_fast*River%run_stor(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_infloc(0) > 0) used = send_data (id_infloc(0), & diag_factor*River%infloc(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_dis(0) > 0) used = send_data (id_dis(0), & diag_factor*River%disw2o(isc:iec,jsc:jec), River%Time) if (id_lake_outflow(0) > 0) used = send_data (id_lake_outflow(0), & diag_factor_2*River%lake_outflow(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_storage_old > 0) used = send_data (id_storage_old, & diag_factor*River%storage(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_stordis_old > 0) used = send_data (id_stordis_old, & diag_factor*River%stordis(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_infloc_old > 0) used = send_data (id_infloc_old, & diag_factor*River%infloc(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_inflow_old > 0) used = send_data (id_inflow_old, & diag_factor*River%inflow(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_outflow_old > 0) used = send_data (id_outflow_old, & diag_factor*River%outflow(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_lake_outflow_old > 0) used = send_data (id_lake_outflow_old, & diag_factor_2*River%lake_outflow(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_disw2o_old > 0) used = send_data (id_disw2o_old, & diag_factor*River%disw2o(isc:iec,jsc:jec), River%Time) if (id_s_frazil > 0) used = send_data (id_s_frazil, & diag_factor*River%storage_c(isc:iec,jsc:jec,1), River%Time, mask=River%mask ) if (id_r_frazil > 0) used = send_data (id_r_frazil, & diag_factor*River%infloc_c(isc:iec,jsc:jec,1), River%Time, mask=River%mask ) if (id_i_frazil > 0) used = send_data (id_i_frazil, & diag_factor*River%inflow_c(isc:iec,jsc:jec,1), River%Time, mask=River%mask ) if (id_o_frazil > 0) used = send_data (id_o_frazil, & diag_factor*River%outflow_c(isc:iec,jsc:jec,1), River%Time, mask=River%mask ) if (id_lofrazil > 0) used = send_data (id_lofrazil, & diag_factor_2*River%lake_outflow_c(isc:iec,jsc:jec,1), River%Time, mask=River%mask ) if (id_dofrazil > 0) used = send_data (id_dofrazil, & diag_factor*River%disc2o(isc:iec,jsc:jec,1), River%Time) if (id_s_t_rivr > 0) used = send_data (id_s_t_rivr, & diag_factor*River%storage_c(isc:iec,jsc:jec,2), River%Time, mask=River%mask ) if (id_r_t_rivr > 0) used = send_data (id_r_t_rivr, & diag_factor*River%infloc_c(isc:iec,jsc:jec,2), River%Time, mask=River%mask ) if (id_i_t_rivr > 0) used = send_data (id_i_t_rivr, & diag_factor*River%inflow_c(isc:iec,jsc:jec,2), River%Time, mask=River%mask ) if (id_o_t_rivr > 0) used = send_data (id_o_t_rivr, & diag_factor*River%outflow_c(isc:iec,jsc:jec,2), River%Time, mask=River%mask ) if (id_lot_rivr > 0) used = send_data (id_lot_rivr, & diag_factor_2*River%lake_outflow_c(isc:iec,jsc:jec,2), River%Time, mask=River%mask ) if (id_dot_rivr > 0) used = send_data (id_dot_rivr, & diag_factor*River%disc2o(isc:iec,jsc:jec,2), River%Time) do i_species = 1, num_species if (id_outflow(i_species) > 0) used = send_data (id_outflow(i_species), & diag_factor*River%outflow_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_lake_outflow(i_species) > 0) used = send_data (id_lake_outflow(i_species), & diag_factor_2*River%lake_outflow_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_inflow(i_species) > 0) used = send_data (id_inflow(i_species), & diag_factor*River%inflow_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_storage(i_species) > 0) used = send_data (id_storage(i_species), & diag_factor*River%storage_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_stordis(i_species) > 0) used = send_data (id_stordis(i_species), & diag_factor*River%stordis_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_run_stor(i_species) > 0) used = send_data (id_run_stor(i_species), & River%dt_fast*River%run_stor_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_infloc(i_species) > 0) used = send_data (id_infloc(i_species), & diag_factor*River%infloc_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_removal(i_species) > 0) used = send_data (id_removal(i_species), & diag_factor*River%removal_c(isc:iec,jsc:jec,i_species), River%Time, mask=River%mask ) if (id_dis(i_species) > 0) used = send_data (id_dis(i_species), & diag_factor*River%disc2o(isc:iec,jsc:jec,i_species), River%Time) enddo if (id_lake_depth_sill > 0) used = send_data (id_lake_depth_sill, & lake_depth_sill, River%Time, mask=River%mask ) if (id_outflowmean > 0) used = send_data (id_outflowmean, & River%outflowmean(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_width > 0) used = send_data (id_width, & River%width(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_depth > 0) used = send_data (id_depth, & River%depth(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_vel > 0) used = send_data (id_vel, & River%vel(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_outflowmean_old > 0) used = send_data (id_outflowmean_old, & River%outflowmean(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_width_old > 0) used = send_data (id_width_old, & River%width(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_depth_old > 0) used = send_data (id_depth_old, & River%depth(isc:iec,jsc:jec), River%Time, mask=River%mask ) if (id_vel_old > 0) used = send_data (id_vel_old, & River%vel(isc:iec,jsc:jec), River%Time, mask=River%mask ) end subroutine river_diag !##################################################################### subroutine get_Leo_Mad_params(DHG_exp, DHG_coef, AAS_exp) type(Leo_Mad_trios), intent(inout) :: DHG_exp ! Exponents for downstream equations type(Leo_Mad_trios), intent(inout) :: DHG_coef ! Coefficients for downstream equations type(Leo_Mad_trios), intent(inout) :: AAS_exp ! Exponents for at-a-station equations !!! Exponents for the downstream hydraulic geometry equations DHG_exp%on_w = ave_DHG_exp(1) DHG_exp%on_d = ave_DHG_exp(2) DHG_exp%on_V = ave_DHG_exp(3) !!! Coefficients for the downstream hydraulic geometry equations DHG_coef%on_w = ave_DHG_coef(1) DHG_coef%on_d = ave_DHG_coef(2) DHG_coef%on_V = ave_DHG_coef(3) !!! Exponents for the at-a-station hydraulic geometry equations AAS_exp%on_w = ave_AAS_exp(1) AAS_exp%on_d = ave_AAS_exp(2) AAS_exp%on_V = ave_AAS_exp(3) end subroutine get_Leo_Mad_params !##################################################################### subroutine river_stock_pe(index, value) integer, intent(in) :: index real , intent(out) :: value ! Domain water (Kg) or heat (Joules) value = 0.0 if (.not.do_rivers) return select case(index) case(ISTOCK_WATER) value = DENS_H2O*(sum(River%storage)+sum(River%stordis)) & + sum(River%run_stor*River%land_area)*River%dt_fast case(ISTOCK_HEAT) ! heat stock not yet implemented value = 0 case default ! Lnd_stock_pe issues a FATAL error message if index is invalid end select end subroutine river_stock_pe !##################################################################### ! returns string indicating the coordiantes of the point i,j function coordinates(i,j) result(s); character(128) :: s integer, intent(in) :: i,j s ='('//trim(string(i))//','//trim(string(j))//')' if (lnd%nfaces>1) s=trim(s)//' on cubic sphere face '//string(lnd%face) end function coordinates end module river_mod #ifdef test_river_solo program river_solo use mpp_mod, only : mpp_error, mpp_pe, mpp_root_pe, mpp_npes, FATAL use mpp_mod, only : mpp_clock_id, mpp_clock_begin, mpp_clock_end use mpp_domains_mod, only : mpp_define_layout, mpp_define_domains use mpp_domains_mod, only : mpp_get_compute_domain, domain2d, CYCLIC_GLOBAL_DOMAIN use mpp_domains_mod, only : mpp_get_current_ntile, mpp_get_tile_id use mpp_io_mod, only : mpp_open, MPP_RDONLY, MPP_NETCDF, MPP_SINGLE use mpp_io_mod, only : MPP_ASCII, MPP_OVERWR, mpp_close use fms_mod, only : fms_init, fms_end, stdlog, open_namelist_file use fms_mod, only : check_nml_error, close_file, file_exist, stdout, read_data use fms_io_mod, only : fms_io_exit use time_manager_mod, only : time_type, increment_time, set_date, increment_date, set_time use time_manager_mod, only : set_calendar_type, JULIAN, NOLEAP, THIRTY_DAY_MONTHS, NO_CALENDAR use time_manager_mod, only : operator(/), operator(-), operator( + ), month_name, get_date use diag_manager_mod, only : diag_manager_init, diag_manager_end use river_mod, only : river_init, river_end, update_river use constants_mod, only : constants_init, PI, radius use grid_mod, only : get_grid_size, get_grid_cell_vertices use grid_mod, only : get_grid_cell_centers, get_grid_cell_area, get_grid_comp_area use grid_mod, only : define_cube_mosaic, get_grid_ntiles use river_mod, only : num_species implicit none real, parameter :: CONST_RUNOFF = 200.0 !--- namelist ----------------------------------------- integer, dimension(6) :: current_date = (/ 0, 0, 0, 0, 0, 0 /) character(len=16) :: calendar = 'julian' integer :: years=0, months=0, days=0, hours=0, minutes=0, seconds=0 integer :: dt_fast = 0 integer ::layout(2) = (/1,0/) namelist /river_solo_nml/ current_date, dt_fast, years, months, days, & hours, minutes, seconds, calendar, layout !-------------------------------------------------------------------- type(time_type) :: Time, Time_start, Time_end, Run_len, Time_step_fast integer :: nf, num_fast_step, unit integer :: yr,mon,day,hr,min,sec, calendar_type=-1 integer :: outunit real, allocatable :: runoff(:,:), discharge(:,:) integer :: initClock, mainClock, termClock, updateClock !Balaji real, allocatable :: runoff_c(:,:,:) real, allocatable :: discharge2ocean(:,:) real, allocatable :: discharge2ocean_c(:,:,:) real, allocatable :: liq(:,:), sol(:,:), mel(:,:), hea(:,:) call fms_init initClock = mpp_clock_id( 'Initialization' ) mainClock = mpp_clock_id( 'Main loop' ) termClock = mpp_clock_id( 'Termination' ) updateClock = mpp_clock_id( 'update river') call mpp_clock_begin(initClock) call river_solo_init call mpp_clock_end (initClock) !end initialization call mpp_clock_begin(mainClock) !begin main loop outunit=stdout() do nf = 1, num_fast_step write(outunit,*)' at river fast time step ', nf call mpp_clock_begin(updateClock) !Balaji call update_river( runoff, runoff_c, discharge2ocean, discharge2ocean_c, & liq, sol, mel, hea ) ! call update_river ( runoff, discharge ) call mpp_clock_end (updateClock) Time = Time + Time_step_fast enddo call mpp_clock_end(mainClock) call mpp_clock_begin(termClock) call river_end call diag_manager_end(Time) call get_date(Time,yr,mon,day,hr,min,sec) if (mpp_pe() == mpp_root_pe()) then call mpp_open(unit, 'RESTART/river_solo.res',form=MPP_ASCII,& action=MPP_OVERWR,threading=MPP_SINGLE,fileset=MPP_SINGLE,nohdrs=.true.) write(unit,*) yr, mon, day, hr, min, sec write(unit,*) calendar_type call mpp_close(unit) endif call fms_io_exit call fms_end contains !####################################################################### subroutine river_solo_init integer :: unit, ierr, io integer :: ni, nj, npes, isc, iec, jsc, jec, ntiles, tile type(domain2d) :: Domain integer, allocatable :: tile_ids(:) ! mosaic tile IDs for the current PE real, allocatable :: lon(:,:), lat(:,:) real, allocatable :: area_lnd(:,:), area_lnd_cell(:,:), gfrac(:,:) integer :: date(6) character(len=9) :: month call constants_init #ifdef INTERNAL_FILE_NML read (input_nml_file, nml=river_solo_nml, iostat=io) ierr = check_nml_error(io, 'river_solo_nml') #else if (file_exist('input.nml')) then unit = open_namelist_file () ierr=1 do while (ierr /= 0) read (unit, nml=river_solo_nml, iostat=io, end=10) ierr = check_nml_error (io, 'river_solo_nml') enddo 10 continue call close_file (unit) endif #endif unit=stdlog() write(unit, nml= river_solo_nml) ! set the calendar if (calendar(1:6) == 'julian') then calendar_type = julian else if (calendar(1:6) == 'NOLEAP') then calendar_type = NOLEAP else if (calendar(1:10) == 'thirty_day') then calendar_type = THIRTY_DAY_MONTHS else if (calendar(1:11) == 'no_calendar') then calendar_type = NO_CALENDAR else if (calendar(1:1) /= ' ') then call mpp_error (FATAL,'==>Error from ocean_solo_mod: invalid namelist value for calendar') else call mpp_error (FATAL,'==>Error from ocean_solo_mod: no namelist value for calendar') endif ! get river_solo restart if (file_exist('INPUT/river_solo.res')) then call mpp_open(unit,'INPUT/river_solo.res',form=MPP_ASCII,action=MPP_RDONLY) read(unit,*) date read(unit,*) calendar_type call close_file(unit) endif call set_calendar_type (calendar_type) call diag_manager_init if (sum(current_date) <= 0) then call mpp_error(FATAL,'==>Error from river_solo_mod: no namelist value for current date') else Time_start = set_date(current_date(1),current_date(2), current_date(3), & current_date(4),current_date(5),current_date(6)) endif if (file_exist('INPUT/river_solo.res')) then Time_start = set_date(date(1),date(2),date(3),date(4),date(5),date(6)) else Time_start = Time_start date = current_date endif Time = Time_start Time_end = increment_date(Time_start, years, months, days, hours, minutes, seconds) Run_len = Time_end - Time_start Time_step_fast = set_time(dt_fast, 0) num_fast_step = Run_len/Time_step_fast call mpp_open (unit, 'time_stamp.out', form=MPP_ASCII, action=MPP_OVERWR,threading=MPP_SINGLE) month = month_name(current_date(2)) if ( mpp_pe() == mpp_root_pe() ) write (unit,'(6i4,2x,a3)') date, month(1:3) call get_date (Time_end, date(1), date(2), date(3), date(4), date(5), date(6)) month = month_name(date(2)) if ( mpp_pe() == mpp_root_pe() ) write (unit,'(6i4,2x,a3)') date, month(1:3) call close_file (unit) !--- get the land grid and set up domain decomposition call get_grid_size('LND', 1, ni, nj) npes = mpp_npes() !--- define domain ------------------------------------------------ call get_grid_ntiles('LND',ntiles) if(layout(1)*layout(2)*ntiles .NE. npes) call mpp_define_layout((/1,ni,1,nj/),npes/ntiles,layout) if (ntiles==1) then call mpp_define_domains ((/1,ni, 1, nj/), layout, domain, & xflags = CYCLIC_GLOBAL_DOMAIN, whalo=1, ehalo=1, shalo=1, nhalo=1, name = 'LAND MODEL') else call define_cube_mosaic ('LND', domain, layout, halo=1 ) endif call mpp_get_compute_domain(Domain, isc, iec, jsc, jec) allocate(tile_ids(mpp_get_current_ntile(Domain))) tile_ids = mpp_get_tile_id(domain) tile = tile_ids(1) deallocate(tile_ids) !--- get grid information allocate( lon(isc:iec,jsc:jec), lat(isc:iec,jsc:jec) ) allocate( area_lnd(ni,nj), area_lnd_cell(ni,nj), gfrac(ni,nj) ) call get_grid_cell_centers ('LND', tile, lon, lat, domain) !!$ call get_grid_cell_area ('LND',tile, area_lnd_cell) !!$ call get_grid_comp_area ('LND',tile, area_lnd) lon = lon * PI/180. lat = lat * PI/180. !!$ gfrac = area_lnd/area_lnd_cell gfrac = 1 npes = mpp_npes() call river_init( lon, lat, Time_start, Time_step_fast, Domain, gfrac(isc:iec,jsc:jec) ) allocate(runoff_c(isc:iec,jsc:jec,num_species) ) allocate(runoff(isc:iec,jsc:jec), discharge(isc:iec,jsc:jec) ) if(file_exist("INPUT/runoff.nc")) then call read_data("INPUT/runoff.nc", "runoff", runoff ) else runoff = CONST_RUNOFF end if if(file_exist("INPUT/runoff.nc")) then call read_data("INPUT/runoff.nc", "runoff_c", runoff_c ) else runoff_c = CONST_RUNOFF end if allocate( discharge2ocean(isc:iec,jsc:jec) ) allocate( discharge2ocean_c(isc:iec,jsc:jec,num_species) ) end subroutine river_solo_init !##################################################################### end program river_solo #endif test_river_solo