diff --git a/cime_config/namelist_definition_mosart.xml b/cime_config/namelist_definition_mosart.xml index a6ad704..86f5224 100644 --- a/cime_config/namelist_definition_mosart.xml +++ b/cime_config/namelist_definition_mosart.xml @@ -288,7 +288,7 @@ -24 - Frequency to perform budget check. Similar to nhtfrq, + Frequency to perform budget check. Similar to nhtfrq, positive means in time steps, 0=monthly, negative means hours (i.e. 24 means every 24 time-steps and -24 means every day diff --git a/cime_config/testdefs/testlist_mosart.xml b/cime_config/testdefs/testlist_mosart.xml index 9e79cbd..168d47b 100644 --- a/cime_config/testdefs/testlist_mosart.xml +++ b/cime_config/testdefs/testlist_mosart.xml @@ -21,6 +21,15 @@ + + + + + + + + + diff --git a/docs/ChangeLog b/docs/ChangeLog index de7e2da..444d6d9 100644 --- a/docs/ChangeLog +++ b/docs/ChangeLog @@ -1,3 +1,26 @@ +=============================================================== +Tag name: mosart1_1_02 +Originator(s): mvertens +Date: Jun 21, 2024 +One-line Summary: cism runoff will be now routed to ocn via mosart + +Enables CISM runoff to be routed to the ocean via mosart. + +All runoff from CISM will be routed directly to the outlet points +New fields will be advertised in the mosart cap +call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofl_glc') +call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofi_glc') + +Issues addressed: +Fixes #92 +Fixes #102 + +Testing: standard testing + izumi -- OK + cheyenne -- OK + +See https://github.com/ESCOMP/MOSART/pull/94 for more details + =============================================================== Tag name: mosart1_1_01 Originator(s): mvertens diff --git a/src/cpl/nuopc/rof_comp_nuopc.F90 b/src/cpl/nuopc/rof_comp_nuopc.F90 index bd222d6..29b75d4 100644 --- a/src/cpl/nuopc/rof_comp_nuopc.F90 +++ b/src/cpl/nuopc/rof_comp_nuopc.F90 @@ -25,14 +25,14 @@ module rof_comp_nuopc model_label_DataInitialize => label_DataInitialize, & model_label_SetRunClock => label_SetRunClock, & model_label_Finalize => label_Finalize, & - SetVM, NUOPC_ModelGet + NUOPC_ModelGet use shr_kind_mod , only : R8=>SHR_KIND_R8, CL=>SHR_KIND_CL, CS=>SHR_KIND_CS use shr_sys_mod , only : shr_sys_abort use shr_log_mod , only : shr_log_getlogunit, shr_log_setlogunit use shr_cal_mod , only : shr_cal_noleap, shr_cal_gregorian, shr_cal_ymd2date use mosart_vars , only : nsrStartup, nsrContinue, nsrBranch, & inst_index, inst_suffix, inst_name, & - mainproc, mpicom_rof, iam, npes, iulog, & + mainproc, mpicom_rof, iam, npes, iulog, vm, & nsrest, caseid, ctitle, version, hostname, username use mosart_data , only : ctl use mosart_driver , only : mosart_read_namelist, mosart_init1, mosart_init2, mosart_run @@ -49,7 +49,6 @@ module rof_comp_nuopc ! Module routines public :: SetServices - public :: SetVM private :: InitializeP0 private :: InitializeAdvertise private :: InitializeRealize @@ -159,7 +158,6 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) type(ESMF_Time) :: refTime ! Ref time type(ESMF_TimeInterval) :: timeStep ! Model timestep type(ESMF_CalKind_Flag) :: esmf_caltype ! esmf calendar type - type(ESMF_VM) :: vm ! esmf virtual machine integer :: ref_ymd ! reference date (YYYYMMDD) integer :: ref_tod ! reference time of day (sec) integer :: yy,mm,dd ! Temporaries for time query @@ -187,6 +185,8 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ! generate local mpi comm !---------------------------------------------------------------------------- + ! Note vm is in mosart_vars.F90 and can be shared among components + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 9cb67db..80cc7c7 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -27,6 +27,7 @@ module rof_import_export private :: state_getimport private :: state_setexport private :: check_for_nans + private :: fldchk type fld_list_type character(len=128) :: stdname @@ -82,9 +83,12 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) isPresent=isPresent, isSet=isSet, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (isPresent .and. isSet) read(cvalue,*) flds_r2l_stream_channel_depths + call fldlist_add(fldsFrRof_num, fldsFrRof, trim(flds_scalar_name)) call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofl') call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofi') + call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofl_glc') + call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofi_glc') call fldlist_add(fldsFrRof_num, fldsFrRof, 'Flrr_flood') call fldlist_add(fldsFrRof_num, fldsFrRof, 'Flrr_volr') call fldlist_add(fldsFrRof_num, fldsFrRof, 'Flrr_volrmch') @@ -109,6 +113,8 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_rofsub') call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_rofi') call fldlist_add(fldsToRof_num, fldsToRof, 'Flrl_irrig') + call fldlist_add(fldsToRof_num, fldsToRof, 'Fgrg_rofl') ! liq runoff from glc + call fldlist_add(fldsToRof_num, fldsToRof, 'Fgrg_rofi') ! ice runoff from glc do n = 1,fldsToRof_num call NUOPC_Advertise(importState, standardName=fldsToRof(n)%stdname, & @@ -222,6 +228,15 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'MOSART' end if + if (fldchk(importState, 'Fgrg_rofl') .and. fldchk(importState, 'Fgrg_rofl')) then + ctl%rof_from_glc = .true. + else + ctl%rof_from_glc = .false. + end if + if (mainproc) then + write(iulog,'(A,l1)') trim(subname) //' rof from glc is ',ctl%rof_from_glc + end if + end subroutine realize_fields !=============================================================================== @@ -239,7 +254,7 @@ subroutine import_fields( gcomp, begr, endr, rc ) ! Local variables type(ESMF_State) :: importState integer :: n,nt - integer :: nliq, nfrz + integer :: nliq, nice character(len=*), parameter :: subname='(rof_import_export:import_fields)' !--------------------------------------------------------------------------- @@ -250,17 +265,8 @@ subroutine import_fields( gcomp, begr, endr, rc ) call NUOPC_ModelGet(gcomp, importState=importState, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ! Set tracers - nliq = 0 - nfrz = 0 - do nt = 1,ctl%ntracers - if (trim(ctl%tracer_names(nt)) == 'LIQ') nliq = nt - if (trim(ctl%tracer_names(nt)) == 'ICE') nfrz = nt - enddo - if (nliq == 0 .or. nfrz == 0) then - write(iulog,*) trim(subname),': ERROR in tracers LIQ ICE ',nliq,nfrz,ctl%tracer_names(:) - call shr_sys_abort() - endif + nliq = ctl%nt_liq + nice = ctl%nt_ice ! determine output array and scale by unit convertsion ! NOTE: the call to state_getimport will convert from input kg/m2s to m3/s @@ -277,7 +283,7 @@ subroutine import_fields( gcomp, begr, endr, rc ) do_area_correction=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call state_getimport(importState, 'Flrl_rofi', begr, endr, ctl%area, output=ctl%qsur(:,nfrz), & + call state_getimport(importState, 'Flrl_rofi', begr, endr, ctl%area, output=ctl%qsur(:,nice), & do_area_correction=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -285,8 +291,20 @@ subroutine import_fields( gcomp, begr, endr, rc ) do_area_correction=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - ctl%qsub(begr:endr, nfrz) = 0.0_r8 - ctl%qgwl(begr:endr, nfrz) = 0.0_r8 + ctl%qsub(begr:endr, nice) = 0.0_r8 + ctl%qgwl(begr:endr, nice) = 0.0_r8 + + if (ctl%rof_from_glc) then + call state_getimport(importState, 'Fgrg_rofl', begr, endr, ctl%area, output=ctl%qglc_liq(:), & + do_area_correction=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getimport(importState, 'Fgrg_rofi', begr, endr, ctl%area, output=ctl%qglc_ice(:), & + do_area_correction=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + else + ctl%qglc_liq(:) = 0._r8 + ctl%qglc_ice(:) = 0._r8 + end if end subroutine import_fields @@ -305,9 +323,11 @@ subroutine export_fields (gcomp, begr, endr, rc) ! Local variables type(ESMF_State) :: exportState integer :: n,nt - integer :: nliq, nfrz + integer :: nliq, nice real(r8) :: rofl(begr:endr) real(r8) :: rofi(begr:endr) + real(r8) :: rofl_glc(begr:endr) + real(r8) :: rofi_glc(begr:endr) real(r8) :: flood(begr:endr) real(r8) :: volr(begr:endr) real(r8) :: volrmch(begr:endr) @@ -325,16 +345,8 @@ subroutine export_fields (gcomp, begr, endr, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return ! Set tracers - nliq = 0 - nfrz = 0 - do nt = 1,ctl%ntracers - if (trim(ctl%tracer_names(nt)) == 'LIQ') nliq = nt - if (trim(ctl%tracer_names(nt)) == 'ICE') nfrz = nt - enddo - if (nliq == 0 .or. nfrz == 0) then - write(iulog,*) trim(subname),': ERROR in tracers LIQ ICE ',nliq,nfrz,ctl%tracer_names(:) - call shr_sys_abort() - endif + nliq = ctl%nt_liq + nice = ctl%nt_ice if (first_time) then if (mainproc) then @@ -351,24 +363,29 @@ subroutine export_fields (gcomp, begr, endr, rc) ! separate liquid and ice runoff do n = begr,endr rofl(n) = ctl%direct(n,nliq) / (ctl%area(n)*0.001_r8) - rofi(n) = ctl%direct(n,nfrz) / (ctl%area(n)*0.001_r8) + rofi(n) = ctl%direct(n,nice) / (ctl%area(n)*0.001_r8) if (ctl%mask(n) >= 2) then ! liquid and ice runoff are treated separately - this is what goes to the ocean rofl(n) = rofl(n) + ctl%runoff(n,nliq) / (ctl%area(n)*0.001_r8) - rofi(n) = rofi(n) + ctl%runoff(n,nfrz) / (ctl%area(n)*0.001_r8) + rofi(n) = rofi(n) + ctl%runoff(n,nice) / (ctl%area(n)*0.001_r8) end if end do else ! liquid and ice runoff added to liquid runoff, ice runoff is zero do n = begr,endr - rofl(n) = (ctl%direct(n,nfrz) + ctl%direct(n,nliq)) / (ctl%area(n)*0.001_r8) + rofl(n) = (ctl%direct(n,nice) + ctl%direct(n,nliq)) / (ctl%area(n)*0.001_r8) if (ctl%mask(n) >= 2) then - rofl(n) = rofl(n) + (ctl%runoff(n,nfrz) + ctl%runoff(n,nliq)) / (ctl%area(n)*0.001_r8) + rofl(n) = rofl(n) + (ctl%runoff(n,nice) + ctl%runoff(n,nliq)) / (ctl%area(n)*0.001_r8) endif rofi(n) = 0._r8 end do end if + do n = begr,endr + rofl_glc(n) = ctl%direct_glc(n,nliq) / (ctl%area(n)*0.001_r8) + rofi_glc(n) = ctl%direct_glc(n,nice) / (ctl%area(n)*0.001_r8) + end do + ! Flooding back to land, sign convention is positive in land->rof direction ! so if water is sent from rof to land, the flux must be negative. ! scs: is there a reason for the wr+wt rather than volr (wr+wt+wh)? @@ -391,6 +408,12 @@ subroutine export_fields (gcomp, begr, endr, rc) call state_setexport(exportState, 'Forr_rofi', begr, endr, input=rofi, do_area_correction=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_setexport(exportState, 'Forr_rofl_glc', begr, endr, input=rofl_glc, do_area_correction=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + call state_setexport(exportState, 'Forr_rofi_glc', begr, endr, input=rofi_glc, do_area_correction=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_setexport(exportState, 'Flrr_flood', begr, endr, input=flood, do_area_correction=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return @@ -654,4 +677,25 @@ subroutine check_for_nans(array, fname, begg) end if end subroutine check_for_nans + !=============================================================================== + logical function fldchk(state, fldname) + ! ---------------------------------------------- + ! Determine if field with fldname is in the input state + ! ---------------------------------------------- + + ! input/output variables + type(ESMF_State), intent(in) :: state + character(len=*), intent(in) :: fldname + + ! local variables + type(ESMF_StateItem_Flag) :: itemFlag + ! ---------------------------------------------- + call ESMF_StateGet(state, trim(fldname), itemFlag) + if (itemflag /= ESMF_STATEITEM_NOTFOUND) then + fldchk = .true. + else + fldchk = .false. + endif + end function fldchk + end module rof_import_export diff --git a/src/riverroute/mosart_budget_type.F90 b/src/riverroute/mosart_budget_type.F90 index 258155f..3df5c8a 100644 --- a/src/riverroute/mosart_budget_type.F90 +++ b/src/riverroute/mosart_budget_type.F90 @@ -3,40 +3,41 @@ module mosart_budget_type ! Variables and routines used for ! calculating and checking tracer global and local budgets - use shr_kind_mod, only: r8 => shr_kind_r8, CL => SHR_KIND_CL - use shr_sys_mod, only: shr_sys_abort - use shr_mpi_mod, only: shr_mpi_sum, shr_mpi_max - use mosart_vars, only: re, spval, barrier_timers, iulog, & - mainproc, npes, iam, mpicom_rof + use shr_kind_mod, only: r8 => shr_kind_r8, CL => SHR_KIND_CL + use shr_sys_mod, only: shr_sys_abort + use shr_mpi_mod, only: shr_mpi_sum, shr_mpi_max + use mosart_vars, only: re, spval, barrier_timers, iulog, mainproc, npes, iam, mpicom_rof + use mosart_data, only: ctl, Tctl, Tunit, TRunoff, Tpara use mosart_timemanager, only: get_nstep, get_curr_date - use mpi implicit none private type budget_type ! accumulated budget over run (not used for now) - real(r8), pointer :: accum_grc(:, :) ! Gridcell level budget accumulator per tracer over the run (m3) - real(r8), pointer :: accum_glob(:) ! Global budget accumulator (1e6 m3) + real(r8), pointer :: accum_grc(:, :) ! Gridcell level budget accumulator per tracer over the run (m3) + real(r8), pointer :: accum_glob(:) ! Global budget accumulator (1e6 m3) ! budget terms per gridcell - real(r8), pointer :: beg_vol_grc(:, :) ! volume begining of the timestep (m3) - real(r8), pointer :: end_vol_grc(:, :) ! volume end of the timestep (m3) - real(r8), pointer :: in_grc(:, :) ! budget in terms (m3) - real(r8), pointer :: out_grc(:, :) ! budget out terms (m3) - real(r8), pointer :: net_grc(:, :) ! net budget (dvolume + inputs - outputs) (m3) - real(r8), pointer :: lag_grc(:, :) ! euler erout lagged (m3) + real(r8), pointer :: beg_vol_grc(:, :) ! volume begining of the timestep (m3) + real(r8), pointer :: end_vol_grc(:, :) ! volume end of the timestep (m3) + real(r8), pointer :: in_grc(:, :) ! budget in terms (m3) + real(r8), pointer :: out_grc(:, :) ! budget out terms (m3) + real(r8), pointer :: net_grc(:, :) ! net budget (dvolume + inputs - outputs) (m3) + real(r8), pointer :: lag_grc(:, :) ! euler erout lagged (m3) + ! budget global terms - real(r8), pointer :: beg_vol_glob(:) ! volume begining of the timestep (1e6 m3) - real(r8), pointer :: end_vol_glob(:) ! volume end of the timestep (1e6 m3) - real(r8), pointer :: in_glob(:) ! budget in terms (1e6 m3) - real(r8), pointer :: out_glob(:) ! budget out terms (1e6 m3) - real(r8), pointer :: net_glob(:) ! net budget (dvolume + inputs - outputs) (1e6 m3) - real(r8), pointer :: lag_glob(:) ! euler erout lagged (1e6 m3) + real(r8), pointer :: beg_vol_glob(:) ! volume begining of the timestep (1e6 m3) + real(r8), pointer :: end_vol_glob(:) ! volume end of the timestep (1e6 m3) + real(r8), pointer :: in_glob(:) ! budget in terms (1e6 m3) + real(r8), pointer :: out_glob(:) ! budget out terms (1e6 m3) + real(r8), pointer :: net_glob(:) ! net budget (dvolume + inputs - outputs) (1e6 m3) + real(r8), pointer :: lag_glob(:) ! euler erout lagged (1e6 m3) + ! budget parameters - real(r8) :: tolerance = 1e-6_r8 ! budget absolute tolerance - real(r8) :: rel_tolerance = 1e-6_r8 ! budget relative tolerance - logical(r8), pointer :: do_budget(:) ! if budget should be checked (per tracer) + real(r8) :: tolerance = 1e-6_r8 ! budget absolute tolerance + real(r8) :: rel_tolerance = 1e-6_r8 ! budget relative tolerance + logical(r8), pointer :: do_budget(:) ! if budget should be checked (per tracer) contains procedure, public :: Init procedure, public :: set_budget @@ -44,20 +45,28 @@ module mosart_budget_type end type budget_type public :: budget_type + integer, parameter :: index_beg_vol_grc = 1 + integer, parameter :: index_end_vol_grc = 2 + integer, parameter :: index_in_grc = 3 + integer, parameter :: index_out_grc = 4 + integer, parameter :: index_net_grc = 5 + integer, parameter :: index_lag_grc = 6 + character(*), parameter :: u_FILE_u = & __FILE__ - !----------------------------------------------------------------------- + +!----------------------------------------------------------------------- contains +!----------------------------------------------------------------------- subroutine Init(this, begr, endr, ntracers) - ! USES: + ! Initialize budget type - ! ARGUMENTS: + ! Arguments class(budget_type) :: this integer, intent(in) :: begr, endr, ntracers - - ! LOCAL VARIABLES: + !------------------------------------------------- ! gridcell level: allocate (this%accum_grc(begr:endr, ntracers)) @@ -108,22 +117,30 @@ subroutine Init(this, begr, endr, ntracers) end subroutine Init + !----------------------------------------------------------------------- + subroutine set_budget(this, begr, endr, ntracers, dt) - !USES: - use mosart_data, only: ctl, Tctl, Tunit, TRunoff, Tpara - !ARGUMENTS: + ! Arguments class(budget_type) :: this integer, intent(in) :: begr, endr, ntracers real(r8), intent(in) :: dt - !LOCAL VARIABLES: - integer nr, nt !indecies + ! local variables + integer :: nr, nt !indices + integer :: nt_liq, nt_ice + !------------------------------------------------- + nt_liq = ctl%nt_liq + nt_ice = ctl%nt_ice do nr = begr, endr do nt = 1, ntracers this%beg_vol_grc(nr, nt) = ctl%volr(nr, nt) - this%in_grc(nr, nt) = (ctl%qsur(nr, nt) + ctl%qsub(nr, nt) + ctl%qgwl(nr, nt)) * dt + if (nt == nt_ice) then + this%in_grc(nr, nt) = (ctl%qsur(nr, nt) + ctl%qsub(nr, nt) + ctl%qgwl(nr, nt) + ctl%qglc_ice(nr)) * dt + else if (nt == nt_liq) then + this%in_grc(nr, nt) = (ctl%qsur(nr, nt) + ctl%qsub(nr, nt) + ctl%qgwl(nr, nt) + ctl%qglc_liq(nr)) * dt + end if ! this was for budget_terms(17) !if (nt==1) then ! this%in_grc(nr,nt)=this%in_grc(nr,nt) +ctl%qirrig(nr) @@ -131,69 +148,71 @@ subroutine set_budget(this, begr, endr, ntracers, dt) end do end do - this%end_vol_grc = 0.0_r8 - this%out_grc = 0.0_r8 - this%net_grc = 0.0_r8 - this%lag_grc = 0.0_r8 + this%end_vol_grc(:,:) = 0.0_r8 + this%out_grc(:,:) = 0.0_r8 + this%net_grc(:,:) = 0.0_r8 + this%lag_grc(:,:) = 0.0_r8 - this%beg_vol_glob = 0.0_r8 - this%end_vol_glob = 0.0_r8 - this%in_glob = 0.0_r8 - this%out_glob = 0.0_r8 - this%net_glob = 0.0_r8 - this%lag_glob = 0.0_r8 + this%beg_vol_glob(:) = 0.0_r8 + this%end_vol_glob(:) = 0.0_r8 + this%in_glob(:) = 0.0_r8 + this%out_glob(:) = 0.0_r8 + this%net_glob(:) = 0.0_r8 + this%lag_glob(:) = 0.0_r8 end subroutine set_budget + !----------------------------------------------------------------------- + subroutine check_budget(this, begr, endr, ntracers, dt) - !USES: - use mosart_data, only: ctl, Tctl, Tunit, TRunoff, Tpara - !ARGUMENTS: + + ! Arguments class(budget_type) :: this integer, intent(in) :: begr, endr, ntracers real(r8), intent(in) :: dt - !LOCAL VARIABLES: - integer nr, nt !indecies - integer yr,mon,day,ymd,tod !time vars - real(r8) :: tmp_in(6, ntracers) ! array to pass to mpi_sum - real(r8) :: tmp_glob(6, ntracers) ! array from mpi_sum - logical :: error_budget ! flag for an error + ! Local variables + integer :: nr, nt !indecies + integer :: nt_liq ! liquid index + integer :: yr,mon,day,ymd,tod !time vars + real(r8) :: tmp_in(6, ntracers) ! array to pass to mpi_sum + real(r8) :: tmp_glob(6, ntracers) ! array from mpi_sum + logical :: error_budget ! flag for an error real(r8) :: abserr, relerr + !------------------------------------------------- call get_curr_date(yr, mon, day, tod) ymd = yr*10000 + mon*100 + day tmp_in = 0.0_r8 tmp_glob = 0.0_r8 + nt_liq = ctl%nt_liq do nr = begr, endr do nt = 1, ntracers this%end_vol_grc(nr, nt) = ctl%volr(nr, nt) - this%out_grc(nr, nt) = this%out_grc(nr, nt) + ctl%direct(nr, nt) - if (nt == 1) then + this%out_grc(nr, nt) = this%out_grc(nr, nt) + ctl%direct(nr, nt) + ctl%direct_glc(nr, nt) + if (nt == nt_liq) then this%out_grc(nr, nt) = this%out_grc(nr, nt) + ctl%flood(nr) end if if (ctl%mask(nr) >= 2) then this%out_grc(nr, nt) = this%out_grc(nr, nt) + ctl%runoff(nr, nt) else - this%lag_grc(nr, nt) = this%lag_grc(nr, nt) - ctl%erout_prev(nr, nt) & - - ctl%flow(nr, nt) + this%lag_grc(nr, nt) = this%lag_grc(nr, nt) - ctl%erout_prev(nr, nt) - ctl%flow(nr, nt) end if this%out_grc(nr,nt) = this%out_grc(nr,nt) * dt this%lag_grc(nr,nt) = this%lag_grc(nr,nt) * dt - this%net_grc(nr, nt) = this%end_vol_grc(nr, nt) - this%beg_vol_grc(nr, nt) & - - (this%in_grc(nr, nt) - this%out_grc(nr, nt)) - this%accum_grc(nr, nt) = this%accum_grc(nr, nt) + this%net_grc(nr, nt) + this%net_grc(nr,nt) = this%end_vol_grc(nr,nt) - this%beg_vol_grc(nr,nt) - (this%in_grc(nr,nt)-this%out_grc(nr,nt)) + this%accum_grc(nr,nt) = this%accum_grc(nr,nt) + this%net_grc(nr,nt) end do end do do nt = 1, ntracers - tmp_in(1, nt) = sum(this%beg_vol_grc(:, nt)) - tmp_in(2, nt) = sum(this%end_vol_grc(:, nt)) - tmp_in(3, nt) = sum(this%in_grc(:, nt)) - tmp_in(4, nt) = sum(this%out_grc(:, nt)) - tmp_in(5, nt) = sum(this%net_grc(:, nt)) - tmp_in(6, nt) = sum(this%lag_grc(:, nt)) + tmp_in(index_beg_vol_grc, nt) = sum(this%beg_vol_grc(:, nt)) + tmp_in(index_end_vol_grc, nt) = sum(this%end_vol_grc(:, nt)) + tmp_in(index_in_grc, nt) = sum(this%in_grc(:, nt)) + tmp_in(index_out_grc, nt) = sum(this%out_grc(:, nt)) + tmp_in(index_net_grc, nt) = sum(this%net_grc(:, nt)) + tmp_in(index_lag_grc, nt) = sum(this%lag_grc(:, nt)) end do tmp_in = tmp_in*1e-6_r8 !convert to million m3 @@ -203,23 +222,22 @@ subroutine check_budget(this, begr, endr, ntracers, dt) error_budget = .false. abserr = 0.0_r8 relerr = 0.0_r8 - this%beg_vol_glob(nt) = tmp_glob(1, nt) - this%end_vol_glob(nt) = tmp_glob(2, nt) - this%in_glob(nt) = tmp_glob(3, nt) - this%out_glob(nt) = tmp_glob(4, nt) - this%net_glob(nt) = tmp_glob(5, nt) - this%lag_glob(nt) = tmp_glob(6, nt) + this%beg_vol_glob(nt) = tmp_glob(index_beg_vol_grc, nt) + this%end_vol_glob(nt) = tmp_glob(index_end_vol_grc, nt) + this%in_glob(nt) = tmp_glob(index_in_grc, nt) + this%out_glob(nt) = tmp_glob(index_out_grc, nt) + this%net_glob(nt) = tmp_glob(index_net_grc, nt) + this%lag_glob(nt) = tmp_glob(index_lag_grc, nt) if (this%do_budget(nt)) then if (abs(this%net_glob(nt) - this%lag_glob(nt)*dt) > this%tolerance) then error_budget = .true. abserr = abs(this%net_glob(nt) - this%lag_glob(nt)) end if if (abs(this%net_glob(nt) + this%lag_glob(nt)) > 1e-6) then - if (abs(this%net_glob(nt) - this%lag_glob(nt)) & - /abs(this%net_glob(nt) + this%lag_glob(nt)) > this%rel_tolerance) then + if ( abs(this%net_glob(nt) - this%lag_glob(nt)) & + /abs(this%net_glob(nt) + this%lag_glob(nt)) > this%rel_tolerance) then error_budget = .true. - relerr = abs(this%net_glob(nt) - this%lag_glob(nt)) & - /abs(this%net_glob(nt) + this%lag_glob(nt)) + relerr = abs(this%net_glob(nt) - this%lag_glob(nt)) /abs(this%net_glob(nt) + this%lag_glob(nt)) end if end if if (mainproc) then diff --git a/src/riverroute/mosart_control_type.F90 b/src/riverroute/mosart_control_type.F90 index 1245545..672613f 100644 --- a/src/riverroute/mosart_control_type.F90 +++ b/src/riverroute/mosart_control_type.F90 @@ -1,17 +1,18 @@ module mosart_control_type - use shr_kind_mod, only : r8 => shr_kind_r8 - use shr_sys_mod, only : shr_sys_abort - use shr_const_mod, only : shr_const_pi, shr_const_rearth - use shr_mpi_mod, only : shr_mpi_sum - use mosart_io, only : ncd_io, ncd_pio_openfile, ncd_pio_closefile - use mosart_vars, only : mainproc, iam, npes, mpicom_rof, iulog, spval, re - use pio, only : file_desc_t, PIO_BCAST_ERROR, pio_seterrorhandling - use ESMF, only : ESMF_DistGrid, ESMF_Array, ESMF_RouteHandle, ESMF_SUCCESS, & - ESMF_DistGridCreate, ESMF_ArrayCreate, ESMF_ArrayHaloStore, & - ESMF_ArrayHalo, ESMF_ArrayGet - use perf_mod, only : t_startf, t_stopf - use nuopc_shr_methods , only : chkerr + use shr_kind_mod, only : r8 => shr_kind_r8, CS => shr_kind_cs + use shr_sys_mod, only : shr_sys_abort + use shr_const_mod, only : shr_const_pi, shr_const_rearth + use shr_string_mod, only : shr_string_listGetNum, shr_string_listGetName + use shr_mpi_mod, only : shr_mpi_sum + use mosart_io, only : ncd_io, ncd_pio_openfile, ncd_pio_closefile + use mosart_vars, only : mainproc, iam, npes, mpicom_rof, iulog, spval, re, vm + use pio, only : file_desc_t, PIO_BCAST_ERROR, pio_seterrorhandling + use ESMF, only : ESMF_DistGrid, ESMF_Array, ESMF_RouteHandle, ESMF_SUCCESS, & + ESMF_DistGridCreate, ESMF_ArrayCreate, ESMF_ArrayHaloStore, & + ESMF_ArrayHalo, ESMF_ArrayGet, ESMF_VMAllReduce, ESMF_REDUCE_SUM + use perf_mod, only : t_startf, t_stopf + use nuopc_shr_methods, only : chkerr implicit none private @@ -27,6 +28,9 @@ module mosart_control_type ! tracers integer :: ntracers = -999 ! number of tracers character(len=3), allocatable :: tracer_names(:)! tracer names + integer :: nt_liq ! index of liquid tracer in tracer_names + integer :: nt_ice ! index of ice tracer in tracer_names + logical :: rof_from_glc ! if true, will receive liq and ice runoff from glc ! decomp info integer :: begr ! local start index @@ -49,13 +53,16 @@ module mosart_control_type real(r8), pointer :: qsur(:,:) => null() ! surface runoff from coupler [m3/s] (lnd) real(r8), pointer :: qsub(:,:) => null() ! subsurfacer runoff from coupler [m3/s] (lnd) real(r8), pointer :: qgwl(:,:) => null() ! glacier/wetland/lake runoff from coupler [m3/s] (lnd) + real(r8), pointer :: qirrig(:) => null() ! irrigation flow from coupler [m3/s] + real(r8), pointer :: qglc_liq(:) => null() ! glacier liquid runoff from coupler [m3/s] (glc) + real(r8), pointer :: qglc_ice(:) => null() ! glacier ice runoff from coupler [m3/s] (glc) ! outputs from MOSART real(r8), pointer :: flood(:) => null() ! flood water to coupler [m3/s] (lnd) real(r8), pointer :: runoff(:,:) => null() ! runoff (from outlet to reach) to coupler [m3/s] - real(r8), pointer :: direct(:,:) => null() ! direct flow to coupler [m3/s] - real(r8), pointer :: qirrig(:) => null() ! irrigation flow to coupler [m3/s] + real(r8), pointer :: direct(:,:) => null() ! direct flow to outlet from land input [m3/s] real(r8), pointer :: qirrig_actual(:) => null() ! minimum of irrigation and available main channel storage [m3/s] + real(r8), pointer :: direct_glc(:,:) =>null() ! direct flow to outlet from glc input [m3/s] ! storage, runoff real(r8), pointer :: runofflnd(:,:) => null() ! runoff masked for land [m3/s] @@ -88,6 +95,7 @@ module mosart_control_type contains procedure, public :: Init + procedure, public :: init_tracer_names procedure, private :: init_decomp procedure, private :: test_halo procedure, public :: calc_gradient @@ -115,13 +123,52 @@ module mosart_control_type integer, public :: halo_w = 7 integer, public :: halo_nw = 8 + ! The following are set from + character(*), parameter :: u_FILE_u = & __FILE__ - !======================================================================== +!======================================================================== contains - !======================================================================== +!======================================================================== + + subroutine init_tracer_names(this, mosart_tracers) + + ! This sets the indices for liquid and ice runoff. In the future additional tracers + ! will be enabled so this is a starting point. + + ! Arguments + class(control_type) :: this + character(len=CS) :: mosart_tracers ! colon delimited string of tracer names + + ! Local variables + integer :: nt + character(len=*),parameter :: subname = '(mosart_control_type: init_tracer_names)' + !----------------------------------------------------------------------- + + ! Determine number of tracers and array of tracer names + this%ntracers = shr_string_listGetNum(mosart_tracers) + allocate(this%tracer_names(this%ntracers)) + do nt = 1,this%ntracers + call shr_string_listGetName(mosart_tracers, nt, this%tracer_names(nt)) + end do + + ! Set tracers + this%nt_liq = 0 + this%nt_ice = 0 + do nt = 1,this%ntracers + if (trim(this%tracer_names(nt)) == 'LIQ') this%nt_liq = nt + if (trim(this%tracer_names(nt)) == 'ICE') this%nt_ice = nt + enddo + if (this%nt_liq == 0 .or. this%nt_ice == 0) then + write(iulog,*) trim(subname),': ERROR in tracers LIQ ICE ',this%nt_liq,this%nt_ice,this%tracer_names(:) + call shr_sys_abort() + endif + + end subroutine init_tracer_names + + !======================================================================== subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) ! Arguments @@ -139,7 +186,8 @@ subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) real(r8) :: rlatn(this%nlat) ! latitude of 1d north grid cell edge (deg) real(r8) :: rlonw(this%nlon) ! longitude of 1d west grid cell edge (deg) real(r8) :: rlone(this%nlon) ! longitude of 1d east grid cell edge (deg) - real(r8) :: larea ! tmp local sum of area + real(r8) :: larea(1) ! tmp local sum of area + real(r8) :: totarea(1) ! tmp total area real(r8) :: deg2rad ! pi/180 integer :: g, n, i, j, nr, nt ! iterators real(r8) :: edgen ! North edge of the direction file @@ -303,6 +351,8 @@ subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) this%qgwl(begr:endr,ntracers), & this%qirrig(begr:endr), & this%qirrig_actual(begr:endr), & + this%qglc_liq(begr:endr), & + this%qglc_ice(begr:endr), & ! this%evel(begr:endr,ntracers), & this%flow(begr:endr,ntracers), & @@ -311,6 +361,7 @@ subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) this%erlat_avg(begr:endr,ntracers), & ! this%effvel(ntracers), & + this%direct_glc(begr:endr,2), & stat=ier) if (ier /= 0) then write(iulog,*)'mosarart_control_type allocation error' @@ -332,11 +383,14 @@ subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) this%qsur(:,:) = 0._r8 this%qsub(:,:) = 0._r8 this%qgwl(:,:) = 0._r8 + this%qglc_liq(:) = 0._r8 + this%qglc_ice(:) = 0._r8 this%fthresh(:) = abs(spval) this%flow(:,:) = 0._r8 this%erout_prev(:,:) = 0._r8 this%eroutup_avg(:,:) = 0._r8 this%erlat_avg(:,:) = 0._r8 + this%direct_glc(:,:) = 0._r8 this%effvel(:) = effvel0 ! downstream velocity (m/s) do nt = 1,ntracers @@ -354,15 +408,19 @@ subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) this%area(nr) = area_global(n) enddo - larea = 0.0_r8 + larea(1) = 0.0_r8 do nr = begr,endr - larea = larea + this%area(nr) + larea(1) = larea(1) + this%area(nr) end do if (minval(this%mask) < 1) then write(iulog,*) subname,'ERROR this mask lt 1 ',minval(this%mask),maxval(this%mask) call shr_sys_abort(subname//' ERROR this mask') endif - call shr_mpi_sum(larea, this%totarea, mpicom_rof, 'mosart totarea', all=.true.) + + call ESMF_VMAllReduce(vm, larea, totarea, 1, ESMF_REDUCE_SUM, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + this%totarea = totarea(1) + if (mainproc) then write(iulog,*) subname,' earth area ',4.0_r8*shr_const_pi*1.0e6_r8*re*re write(iulog,*) subname,' mosart area ',this%totarea @@ -1175,7 +1233,7 @@ subroutine calc_gradient(this, fld, fld_halo_array, dfld_dx, dfld_dy, rc) dfld_dx(n) = dfld_dx(n) + (fld_surrounding(ax_indices(i)) - fld_surrounding(sx_indices(i))) dfld_dy(n) = dfld_dy(n) + (fld_surrounding(ay_indices(i)) - fld_surrounding(sy_indices(i))) enddo - + dfld_dx(n) = dfld_dx(n) / (8._r8*mean_dx) dfld_dy(n) = dfld_dy(n) / (8._r8*mean_dy) diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 83fc1ee..652bf98 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -7,7 +7,6 @@ module mosart_driver use shr_kind_mod , only : r8 => shr_kind_r8, CS => shr_kind_cs, CL => shr_kind_CL use shr_sys_mod , only : shr_sys_abort use shr_const_mod , only : SHR_CONST_PI, SHR_CONST_CDAY - use shr_string_mod , only : shr_string_listGetNum, shr_string_listGetName use mosart_vars , only : re, spval, iulog, ice_runoff, & frivinp, nsrContinue, nsrBranch, nsrStartup, nsrest, & inst_index, inst_suffix, inst_name, decomp_option, & @@ -43,11 +42,11 @@ module mosart_driver public :: mosart_run ! River routing model ! mosart namelists - integer :: coupling_period ! mosart coupling period - integer :: delt_mosart ! mosart internal timestep (->nsub) - logical :: use_halo_option ! enable halo capability using ESMF - character(len=CS) :: mosart_tracers ! colon delimited string of tracer names - character(len=CS) :: mosart_euler_calc ! colon delimited string of logicals for using Euler algorithm + integer :: coupling_period ! mosart coupling period + integer :: delt_mosart ! mosart internal timestep (->nsub) + logical :: use_halo_option ! enable halo capability using ESMF + character(len=CS) :: mosart_tracers ! colon delimited string of tracer names + character(len=CS) :: mosart_euler_calc ! colon delimited string of logicals for using Euler algorithm ! subcycling integer :: nsub_save ! previous nsub @@ -62,6 +61,8 @@ module mosart_driver character(len=CL) :: nlfilename_rof = 'mosart_in' character(len=CL) :: fnamer ! name of netcdf restart file + integer :: nt_liq, nt_ice + character(*), parameter :: u_FILE_u = & __FILE__ !----------------------------------------------------------------------- @@ -148,12 +149,10 @@ subroutine mosart_read_namelist() call mpi_bcast (mosart_euler_calc, CS, MPI_CHARACTER, 0, mpicom_rof, ier) call mpi_bcast (budget_frq,1,MPI_INTEGER,0,mpicom_rof,ier) - ! Determine number of tracers and array of tracer names - ctl%ntracers = shr_string_listGetNum(mosart_tracers) - allocate(ctl%tracer_names(ctl%ntracers)) - do i = 1,ctl%ntracers - call shr_string_listGetName(mosart_tracers, i, ctl%tracer_names(i)) - end do + ! Determine number of tracers and array of tracer names and initialize module variables + call ctl%init_tracer_names(mosart_tracers) + nt_liq = ctl%nt_liq + nt_ice = ctl%nt_ice runtyp(:) = 'missing' runtyp(nsrStartup + 1) = 'initial' @@ -162,17 +161,17 @@ subroutine mosart_read_namelist() if (mainproc) then write(iulog,*) 'define run:' - write(iulog,'(a)' ) ' run type = '//trim(runtyp(nsrest+1)) - write(iulog,'(a,i8)') ' coupling_period = ',coupling_period - write(iulog,'(a,i8)') ' delt_mosart = ',delt_mosart - write(iulog,'(a)' ) ' decomp option = '//trim(decomp_option) - write(iulog,'(a,l1)') ' use_halo_option = ',use_halo_option - write(iulog,'(a)' ) ' bypass_routing option = '//trim(bypass_routing_option) - write(iulog,'(a)' ) ' qgwl runoff option = '//trim(qgwl_runoff_option) - write(iulog,'(a)' ) ' mosart tracers = '//trim(mosart_tracers) - write(iulog,'(a)' ) ' mosart euler calc = '//trim(mosart_euler_calc) + write(iulog,'(a)' ) ' run type = '//trim(runtyp(nsrest+1)) + write(iulog,'(a,i8)') ' coupling_period = ',coupling_period + write(iulog,'(a,i8)') ' delt_mosart = ',delt_mosart + write(iulog,'(a)' ) ' decomp option = '//trim(decomp_option) + write(iulog,'(a,l1)') ' use_halo_option = ',use_halo_option + write(iulog,'(a)' ) ' bypass_routing option = '//trim(bypass_routing_option) + write(iulog,'(a)' ) ' qgwl runoff option = '//trim(qgwl_runoff_option) + write(iulog,'(a)' ) ' mosart tracers = '//trim(mosart_tracers) + write(iulog,'(a)' ) ' mosart euler calc = '//trim(mosart_euler_calc) if (nsrest == nsrStartup .and. finidat /= ' ') then - write(iulog,'(a)') ' mosart initial data = '//trim(finidat) + write(iulog,'(a)') ' mosart initial data = '//trim(finidat) end if endif @@ -466,8 +465,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) call t_stopf('mosartr_budgetset') endif - - ! data for euler solver, in m3/s here + ! initialize data for euler solver, in m3/s here do nr = begr,endr do nt = 1,ntracers TRunoff%qsur(nr,nt) = ctl%qsur(nr,nt) @@ -483,7 +481,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !----------------------------------- call t_startf('mosartr_irrig') - nt = 1 ctl%qirrig_actual = 0._r8 do nr = begr,endr @@ -492,10 +489,10 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ! compare irrig_volume to main channel storage; ! add overage to subsurface runoff - if(irrig_volume > TRunoff%wr(nr,nt)) then - ctl%qsub(nr,nt) = ctl%qsub(nr,nt) + (TRunoff%wr(nr,nt) - irrig_volume) / coupling_period - TRunoff%qsub(nr,nt) = ctl%qsub(nr,nt) - irrig_volume = TRunoff%wr(nr,nt) + if(irrig_volume > TRunoff%wr(nr,nt_liq)) then + ctl%qsub(nr,nt_liq) = ctl%qsub(nr,nt_liq) + (TRunoff%wr(nr,nt_liq) - irrig_volume) / coupling_period + TRunoff%qsub(nr,nt_liq) = ctl%qsub(nr,nt_liq) + irrig_volume = TRunoff%wr(nr,nt_liq) endif ! actual irrigation rate [m3/s] @@ -504,7 +501,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ctl%qirrig_actual(nr) = - irrig_volume / coupling_period ! remove irrigation from wr (main channel) - TRunoff%wr(nr,nt) = TRunoff%wr(nr,nt) - irrig_volume + TRunoff%wr(nr,nt_liq) = TRunoff%wr(nr,nt_liq) - irrig_volume enddo call t_stopf('mosartr_irrig') @@ -517,14 +514,13 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !----------------------------------- call t_startf('mosartr_flood') - nt = 1 ctl%flood = 0._r8 do nr = begr,endr ! initialize ctl%flood to zero if (ctl%mask(nr) == 1) then - if (ctl%volr(nr,nt) > ctl%fthresh(nr)) then + if (ctl%volr(nr,nt_liq) > ctl%fthresh(nr)) then ! determine flux that is sent back to the land this is in m3/s - ctl%flood(nr) = (ctl%volr(nr,nt)-ctl%fthresh(nr)) / (delt_coupling) + ctl%flood(nr) = (ctl%volr(nr,nt_liq)-ctl%fthresh(nr)) / (delt_coupling) ! ctl%flood will be sent back to land - so must subtract this ! from the input runoff from land @@ -535,7 +531,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ! it at the end or even during the run loop as the ! new volume is computed. fluxout depends on volr, so ! how this is implemented does impact the solution. - TRunoff%qsur(nr,nt) = TRunoff%qsur(nr,nt) - ctl%flood(nr) + TRunoff%qsur(nr,nt_liq) = TRunoff%qsur(nr,nt_liq) - ctl%flood(nr) endif endif enddo @@ -564,45 +560,80 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) if (chkerr(rc,__LINE__,u_FILE_u)) return !----------------------------------------------------- - !--- all frozen runoff passed direct to outlet + !--- initialize ctl%direct + !----------------------------------------------------- + + ctl%direct(:,:) = 0._r8 + + !----------------------------------------------------- + !--- direct to outlet: all liquid and frozen runoff from glc + !----------------------------------------------------- + + if (ctl%rof_from_glc) then + src_direct(:,:) = 0._r8 + dst_direct(:,:) = 0._r8 + + cnt = 0 + do nr = begr,endr + cnt = cnt + 1 + src_direct(nt_liq,cnt) = ctl%qglc_liq(nr) + src_direct(nt_ice,cnt) = ctl%qglc_ice(nr) + enddo + + call ESMF_FieldSMM(Tunit%srcfield, Tunit%dstfield, Tunit%rh_direct, termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + + ! copy direct transfer water to output field + cnt = 0 + do nr = begr,endr + cnt = cnt + 1 + ctl%direct_glc(nr,nt_liq) = dst_direct(nt_liq,cnt) + ctl%direct_glc(nr,nt_ice) = dst_direct(nt_ice,cnt) + enddo + else + ctl%direct_glc(:,:) = 0._r8 + ctl%direct_glc(:,:) = 0._r8 + end if + + !----------------------------------------------------- + !--- direct to outlet: all frozen runoff from lnd !----------------------------------------------------- - nt = 2 src_direct(:,:) = 0._r8 dst_direct(:,:) = 0._r8 - ! set euler_calc = false for frozen runoff - ! TODO: will be reworked after addition of multiple tracers - Tunit%euler_calc(nt) = .false. - cnt = 0 do nr = begr,endr cnt = cnt + 1 - src_direct(nt,cnt) = TRunoff%qsur(nr,nt) + TRunoff%qsub(nr,nt) + TRunoff%qgwl(nr,nt) - TRunoff%qsur(nr,nt) = 0._r8 - TRunoff%qsub(nr,nt) = 0._r8 - TRunoff%qgwl(nr,nt) = 0._r8 + src_direct(nt_ice,cnt) = TRunoff%qsur(nr,nt_ice) + TRunoff%qsub(nr,nt_ice) + TRunoff%qgwl(nr,nt_ice) enddo call ESMF_FieldSMM(Tunit%srcfield, Tunit%dstfield, Tunit%rh_direct, termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return ! copy direct transfer water to output field - ctl%direct = 0._r8 cnt = 0 do nr = begr,endr cnt = cnt + 1 - ctl%direct(nr,nt) = ctl%direct(nr,nt) + dst_direct(nt,cnt) + ctl%direct(nr,nt_ice) = ctl%direct(nr,nt_ice) + dst_direct(nt_ice,cnt) enddo + ! set euler_calc = false for frozen runoff + ! TODO: will be reworked after addition of multiple tracers + Tunit%euler_calc(nt_ice) = .false. + + ! Set Trunoff%qsur, TRunoff%qsub and Trunoff%qgwl to zero for nt_ice + TRunoff%qsur(:,nt_ice) = 0._r8 + TRunoff%qsub(:,nt_ice) = 0._r8 + TRunoff%qgwl(:,nt_ice) = 0._r8 + !----------------------------------------------------- - !--- direct to outlet qgwl + !--- direct to outlet: qgwl !----------------------------------------------------- !-- liquid runoff components if (trim(bypass_routing_option) == 'direct_to_outlet') then - nt = 1 src_direct(:,:) = 0._r8 dst_direct(:,:) = 0._r8 @@ -611,12 +642,12 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) do nr = begr,endr cnt = cnt + 1 if (trim(qgwl_runoff_option) == 'all') then - src_direct(nt,cnt) = TRunoff%qgwl(nr,nt) - TRunoff%qgwl(nr,nt) = 0._r8 + src_direct(nt_liq,cnt) = TRunoff%qgwl(nr,nt_liq) + TRunoff%qgwl(nr,nt_liq) = 0._r8 else if (trim(qgwl_runoff_option) == 'negative') then - if(TRunoff%qgwl(nr,nt) < 0._r8) then - src_direct(nt,cnt) = TRunoff%qgwl(nr,nt) - TRunoff%qgwl(nr,nt) = 0._r8 + if(TRunoff%qgwl(nr,nt_liq) < 0._r8) then + src_direct(nt_liq,cnt) = TRunoff%qgwl(nr,nt_liq) + TRunoff%qgwl(nr,nt_liq) = 0._r8 endif endif enddo @@ -628,63 +659,57 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) cnt = 0 do nr = begr,endr cnt = cnt + 1 - ctl%direct(nr,nt) = ctl%direct(nr,nt) + dst_direct(nt,cnt) + ctl%direct(nr,nt_liq) = ctl%direct(nr,nt_liq) + dst_direct(nt_liq,cnt) enddo endif !----------------------------------------------------- - !--- direct in place qgwl + !--- direct in place qgwl, qgwl !----------------------------------------------------- if (trim(bypass_routing_option) == 'direct_in_place') then - - nt = 1 do nr = begr,endr - if (trim(qgwl_runoff_option) == 'all') then - ctl%direct(nr,nt) = TRunoff%qgwl(nr,nt) - TRunoff%qgwl(nr,nt) = 0._r8 + ctl%direct(nr,nt_liq) = TRunoff%qgwl(nr,nt_liq) + TRunoff%qgwl(nr,nt_liq) = 0._r8 else if (trim(qgwl_runoff_option) == 'negative') then - if(TRunoff%qgwl(nr,nt) < 0._r8) then - ctl%direct(nr,nt) = TRunoff%qgwl(nr,nt) - TRunoff%qgwl(nr,nt) = 0._r8 + if(TRunoff%qgwl(nr,nt_liq) < 0._r8) then + ctl%direct(nr,nt_liq) = TRunoff%qgwl(nr,nt_liq) + TRunoff%qgwl(nr,nt_liq) = 0._r8 endif else if (trim(qgwl_runoff_option) == 'threshold') then ! --- calculate volume of qgwl flux during timestep - qgwl_volume = TRunoff%qgwl(nr,nt) * ctl%area(nr) * coupling_period + qgwl_volume = TRunoff%qgwl(nr,nt_liq) * ctl%area(nr) * coupling_period river_volume_minimum = river_depth_minimum * ctl%area(nr) ! if qgwl is negative, and adding it to the main channel ! would bring main channel storage below a threshold, ! send qgwl directly to ocean - if (((qgwl_volume + TRunoff%wr(nr,nt)) < river_volume_minimum) .and. (TRunoff%qgwl(nr,nt) < 0._r8)) then - ctl%direct(nr,nt) = TRunoff%qgwl(nr,nt) - TRunoff%qgwl(nr,nt) = 0._r8 + if (((qgwl_volume + TRunoff%wr(nr,nt_liq)) < river_volume_minimum) .and. (TRunoff%qgwl(nr,nt_liq) < 0._r8)) then + ctl%direct(nr,nt_liq) = TRunoff%qgwl(nr,nt_liq) + TRunoff%qgwl(nr,nt_liq) = 0._r8 endif endif enddo - endif !------------------------------------------------------- - !--- add other direct terms, e.g. inputs outside of - !--- mosart mask, negative qsur + !--- direct in place: add other direct terms, e.g. inputs outside of mosart mask, negative qsur !------------------------------------------------------- if (trim(bypass_routing_option) == 'direct_in_place') then do nt = 1,ntracers do nr = begr,endr - if (TRunoff%qsub(nr,nt) < 0._r8) then ctl%direct(nr,nt) = ctl%direct(nr,nt) + TRunoff%qsub(nr,nt) TRunoff%qsub(nr,nt) = 0._r8 endif - if (TRunoff%qsur(nr,nt) < 0._r8) then ctl%direct(nr,nt) = ctl%direct(nr,nt) + TRunoff%qsur(nr,nt) TRunoff%qsur(nr,nt) = 0._r8 endif - + ! Note Tunit%mask is set in Tunit%init and is obtained from reading in fdir + ! if fdir<0 then mask=0 (ocean), if fdir=0 then mask=2 (outlet) and if fdir>0 then mask=1 (land) if (Tunit%mask(nr) > 0) then ! mosart euler else @@ -697,11 +722,13 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) enddo endif - if (trim(bypass_routing_option) == 'direct_to_outlet') then + !------------------------------------------------------- + !--- direct to outlet: add other direct terms, e.g. inputs outside of mosart mask, negative qsur + !------------------------------------------------------- + if (trim(bypass_routing_option) == 'direct_to_outlet') then src_direct(:,:) = 0._r8 dst_direct(:,:) = 0._r8 - cnt = 0 do nr = begr,endr cnt = cnt + 1 @@ -720,15 +747,20 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !---- water outside the basin --- !---- *** DO NOT TURN THIS ONE OFF, conservation will fail *** --- + + ! Note Tunit%mask is set in Tunit%init and is obtained from reading in fdir + ! if fdir<0 then mask=0 (ocean), if fdir=0 then mask=2 (outlet) and if fdir>0 then mask=1 (land) if (Tunit%mask(nr) > 0) then ! mosart euler else - src_direct(nt,cnt) = src_direct(nt,cnt) + TRunoff%qsub(nr,nt) + TRunoff%qsur(nr,nt) & - + TRunoff%qgwl(nr,nt) + ! NOTE: that when nt = nt_ice, the TRunoff terms + ! below have already been set to zero in the frozen + ! runoff calculation above - where frozen runoff is always set to the outlet + src_direct(nt,cnt) = src_direct(nt,cnt) + TRunoff%qsub(nr,nt) + TRunoff%qsur(nr,nt) + TRunoff%qgwl(nr,nt) TRunoff%qsub(nr,nt) = 0._r8 TRunoff%qsur(nr,nt) = 0._r8 TRunoff%qgwl(nr,nt) = 0._r8 - endif + end if enddo enddo @@ -817,7 +849,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ctl%erlat_avg = ctl%erlat_avg / float(nsub) ! update states when subsycling completed - ! TODO: move of this to hist_set_flds ctl%runoff = 0._r8 ctl%runofflnd = spval ctl%runoffocn = spval @@ -841,12 +872,19 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) endif enddo enddo + + ! final update from glc input + do nr = begr,endr + ctl%runofftot(nr,nt_liq) = ctl%runofftot(nr,nt_liq) + ctl%direct_glc(nr,nt_liq) + ctl%runofftot(nr,nt_ice) = ctl%runofftot(nr,nt_ice) + ctl%direct_glc(nr,nt_ice) + end do + call t_stopf('mosartr_subcycling') !----------------------------------- ! BUDGET !----------------------------------- - if (budget_check) then + if (budget_check) then call t_startf('mosartr_budgetcheck') call budget%check_budget(begr,endr,ntracers,delt_coupling) call t_stopf('mosartr_budgetcheck') diff --git a/src/riverroute/mosart_histfile.F90 b/src/riverroute/mosart_histfile.F90 index 4872b34..c0eaf23 100644 --- a/src/riverroute/mosart_histfile.F90 +++ b/src/riverroute/mosart_histfile.F90 @@ -644,16 +644,12 @@ subroutine htape_create (t, histrest) call ncd_putatt(lnfid, ncd_global, 'username' , trim(username)) call ncd_putatt(lnfid, ncd_global, 'version' , trim(version)) call ncd_putatt(lnfid, ncd_global, 'model_doi_url', trim(model_doi_url)) - write(6,*)'DEBUG: I am here7' call ncd_putatt(lnfid, ncd_global, 'case_title', trim(ctitle)) - write(6,*)'DEBUG: I am here8' call ncd_putatt(lnfid, ncd_global, 'case_id', trim(caseid)) - write(6,*)'DEBUG: I am here9' str = get_filename(frivinp) call ncd_putatt(lnfid, ncd_global, 'input_dataset', trim(str)) - write(6,*)'DEBUG: I am here10' ! ! add global attribute time_period_freq @@ -679,7 +675,6 @@ subroutine htape_create (t, histrest) 999 format(a,i0) call ncd_putatt(lnfid, ncd_global, 'time_period_freq', trim(time_period_freq)) - write(6,*)'DEBUG: I am here6' ! Define dimensions. ! Time is an unlimited dimension. Character string is treated as an array of characters. @@ -689,12 +684,10 @@ subroutine htape_create (t, histrest) call ncd_defdim(lnfid, 'lat' , ctl%nlat , dimid) call ncd_defdim(lnfid, 'allrof', ctl%numr , dimid) call ncd_defdim(lnfid, 'string_length', 8, strlen_dimid) - write(6,*)'DEBUG: I am here7' if ( .not. lhistrest )then call ncd_defdim(lnfid, 'hist_interval', 2, hist_interval_dimid) call ncd_defdim(lnfid, 'time', ncd_unlimited, time_dimid) - write(6,*)'DEBUG: I am here8' if (mainproc)then write(iulog,*) trim(subname),' : Successfully defined netcdf history file ',t end if diff --git a/src/riverroute/mosart_histflds.F90 b/src/riverroute/mosart_histflds.F90 index 31287df..18534f1 100644 --- a/src/riverroute/mosart_histflds.F90 +++ b/src/riverroute/mosart_histflds.F90 @@ -23,6 +23,7 @@ module mosart_histflds type(hist_pointer_type), allocatable :: h_runoffocn(:) type(hist_pointer_type), allocatable :: h_runofftot(:) type(hist_pointer_type), allocatable :: h_direct(:) + type(hist_pointer_type), allocatable :: h_direct_glc(:) type(hist_pointer_type), allocatable :: h_dvolrdtlnd(:) type(hist_pointer_type), allocatable :: h_dvolrdtocn(:) type(hist_pointer_type), allocatable :: h_volr(:) @@ -31,6 +32,8 @@ module mosart_histflds type(hist_pointer_type), allocatable :: h_qgwl(:) real(r8), pointer :: h_volr_mch(:) + real(r8), pointer :: h_qglc_liq_input(:) + real(r8), pointer :: h_qglc_ice_input(:) !------------------------------------------------------------------------ contains @@ -60,6 +63,7 @@ subroutine mosart_histflds_init(begr, endr, ntracers) allocate(h_qsur(ntracers)) allocate(h_qsub(ntracers)) allocate(h_qgwl(ntracers)) + allocate(h_direct_glc(2)) do nt = 1,ntracers allocate(h_runofflnd(nt)%data(begr:endr)) @@ -73,8 +77,12 @@ subroutine mosart_histflds_init(begr, endr, ntracers) allocate(h_qsub(nt)%data(begr:endr)) allocate(h_qgwl(nt)%data(begr:endr)) end do + allocate(h_direct_glc(ctl%nt_liq)%data(begr:endr)) + allocate(h_direct_glc(ctl%nt_ice)%data(begr:endr)) allocate(h_volr_mch(begr:endr)) + allocate(h_qglc_liq_input(begr:endr)) + allocate(h_qglc_ice_input(begr:endr)) !------------------------------------------------------- ! Build master field list of all possible fields in a history file. @@ -91,7 +99,7 @@ subroutine mosart_histflds_init(begr, endr, ntracers) call mosart_hist_addfld (fname='RIVER_DISCHARGE_TO_OCEAN'//'_'//trim(ctl%tracer_names(nt)), units='m3/s', & avgflag='A', long_name='MOSART river discharge into ocean: '//trim(ctl%tracer_names(nt)), & - ptr_rof=h_runoffocn(nt)%data, default='inactive') + ptr_rof=h_runoffocn(nt)%data, default='active') call mosart_hist_addfld (fname='TOTAL_DISCHARGE_TO_OCEAN'//'_'//trim(ctl%tracer_names(nt)), units='m3/s', & avgflag='A', long_name='MOSART total discharge into ocean: '//trim(ctl%tracer_names(nt)), & @@ -101,6 +109,10 @@ subroutine mosart_histflds_init(begr, endr, ntracers) avgflag='A', long_name='MOSART direct discharge into ocean: '//trim(ctl%tracer_names(nt)), & ptr_rof=h_direct(nt)%data, default='active') + call mosart_hist_addfld (fname='DIRECT_DISCHARGE_TO_OCEAN_GLC'//'_'//trim(ctl%tracer_names(nt)), units='m3/s', & + avgflag='A', long_name='MOSART direct discharge into ocean from glc: '//trim(ctl%tracer_names(nt)), & + ptr_rof=h_direct_glc(nt)%data, default='active') + call mosart_hist_addfld (fname='STORAGE'//'_'//trim(ctl%tracer_names(nt)), units='m3', & avgflag='A', long_name='MOSART storage: '//trim(ctl%tracer_names(nt)), & ptr_rof=h_volr(nt)%data, default='inactive') @@ -138,6 +150,14 @@ subroutine mosart_histflds_init(begr, endr, ntracers) avgflag='A', long_name='Actual irrigation (if limited by river storage)', & ptr_rof=ctl%qirrig_actual, default='inactive') + call mosart_hist_addfld (fname='QGLC_LIQ_INPUT', units='m3', & + avgflag='A', long_name='liquid runoff from glc input', & + ptr_rof=h_qglc_liq_input, default='active') + + call mosart_hist_addfld (fname='QGLC_ICE_INPUT', units='m3', & + avgflag='A', long_name='ice runoff from glc input', & + ptr_rof=h_qglc_ice_input, default='active') + ! print masterlist of history fields call mosart_hist_printflds() @@ -156,6 +176,10 @@ subroutine mosart_histflds_set(ntracers) ! Local variables integer :: nt + integer :: nt_liq, nt_ice + + nt_liq = ctl%nt_liq + nt_ice = ctl%nt_ice do nt = 1,ntracers h_runofflnd(nt)%data(:) = ctl%runofflnd(:,nt) @@ -169,6 +193,10 @@ subroutine mosart_histflds_set(ntracers) h_qgwl(nt)%data(:) = ctl%qgwl(:,nt) end do h_volr_mch(:) = Trunoff%wr(:,1) + h_qglc_liq_input(:) = ctl%qglc_liq(:) + h_qglc_ice_input(:) = ctl%qglc_ice(:) + h_direct_glc(nt_liq)%data(:) = ctl%direct_glc(:,nt_liq) + h_direct_glc(nt_ice)%data(:) = ctl%direct_glc(:,nt_ice) end subroutine mosart_histflds_set diff --git a/src/riverroute/mosart_physics.F90 b/src/riverroute/mosart_physics.F90 index 4350f64..3700d42 100644 --- a/src/riverroute/mosart_physics.F90 +++ b/src/riverroute/mosart_physics.F90 @@ -28,6 +28,15 @@ module mosart_physics public :: subnetworkrouting public :: mainchannelrouting + private :: Routing_KW + private :: CRVRMAN_nosqrt + private :: CREHT_nosqrt + private :: GRMR + private :: GRHT + private :: GRPT + private :: GRRR + private :: GRPR + real(r8), parameter :: TINYVALUE = 1.0e-14_r8 ! double precision variable has a significance of about 16 decimal digits real(r8), parameter :: SLOPE1def = 0.1_r8 ! here give it a small value in order to avoid the abrupt change of hydraulic radidus etc. real(r8) :: sinatanSLOPE1defr ! 1.0/sin(atan(slope1)) @@ -270,14 +279,8 @@ subroutine mainchannelRouting(nr, nt, theDeltaT) if(Tctl%RoutingMethod == 1) then call Routing_KW(nr, nt, theDeltaT) - else if(Tctl%RoutingMethod == 2) then - call Routing_MC(nr, nt, theDeltaT) - else if(Tctl%RoutingMethod == 3) then - call Routing_THREW(nr, nt, theDeltaT) - else if(Tctl%RoutingMethod == 4) then - call Routing_DW(nr, nt, theDeltaT) else - call shr_sys_abort( "mosart: Please check the routing method! There are only 4 methods available." ) + call shr_sys_abort( "mosart: Please check the routing method! There is only 1 method currently available." ) end if end subroutine mainchannelRouting @@ -346,39 +349,6 @@ end subroutine Routing_KW !----------------------------------------------------------------------- - subroutine Routing_MC(nr, nt, theDeltaT) - ! Muskingum-Cunge routing method - - ! Arguments - integer, intent(in) :: nr, nt - real(r8), intent(in) :: theDeltaT - - end subroutine Routing_MC - - !----------------------------------------------------------------------- - - subroutine Routing_THREW(nr, nt, theDeltaT) - ! kinematic wave routing method from THREW model - - ! Arguments - integer, intent(in) :: nr, nt - real(r8), intent(in) :: theDeltaT - - end subroutine Routing_THREW - - !----------------------------------------------------------------------- - - subroutine Routing_DW(nr, nt, theDeltaT) - ! classic diffusion wave routing method - - ! Arguments - integer, intent(in) :: nr, nt - real(r8), intent(in) :: theDeltaT - - end subroutine Routing_DW - - !----------------------------------------------------------------------- - subroutine updateState_hillslope(nr,nt) ! update the state variables at hillslope diff --git a/src/riverroute/mosart_tctl_type.F90 b/src/riverroute/mosart_tctl_type.F90 index ce35168..3571086 100644 --- a/src/riverroute/mosart_tctl_type.F90 +++ b/src/riverroute/mosart_tctl_type.F90 @@ -10,7 +10,7 @@ module mosart_tctl_type integer :: DLevelH2R ! The base number of channel routing sub-time-steps within one hillslope routing step. ! Usually channel routing requires small time steps than hillslope routing. integer :: DLevelR ! The number of channel routing sub-time-steps at a higher level within one channel routing step at a lower level. - integer :: RoutingMethod ! Flag for routing methods. 1 --> variable storage method from SWAT model; 2 --> Muskingum method? + integer :: RoutingMethod ! Flag for routing methods. 1 --> variable storage method from SWAT model contains procedure :: Init end type Tctl_type diff --git a/src/riverroute/mosart_timemanager.F90 b/src/riverroute/mosart_timemanager.F90 index 1d35ae7..df53ba5 100644 --- a/src/riverroute/mosart_timemanager.F90 +++ b/src/riverroute/mosart_timemanager.F90 @@ -3,7 +3,7 @@ module mosart_timemanager use shr_kind_mod , only: r8 => shr_kind_r8, CS => shr_kind_CS use shr_sys_mod , only: shr_sys_abort use shr_string_mod , only: shr_string_toUpper - use mosart_vars , only: isecspday, iulog, nsrest, nsrContinue, mpicom_rof, mainproc + use mosart_vars , only: isecspday, iulog, nsrest, nsrContinue, mainproc use ESMF , only: ESMF_MAXSTR, ESMF_Calendar, ESMF_Clock, ESMF_Time, ESMF_TimeInterval, & ESMF_TimeIntervalSet, ESMF_TimeIntervalGet, ESMF_TimeSet, ESMF_TimeGet, & ESMF_ClockCreate, ESMF_ClockGet, ESMF_ClockAdvance, & diff --git a/src/riverroute/mosart_vars.F90 b/src/riverroute/mosart_vars.F90 index cc0cc48..6712c4d 100644 --- a/src/riverroute/mosart_vars.F90 +++ b/src/riverroute/mosart_vars.F90 @@ -3,6 +3,7 @@ module mosart_vars use shr_kind_mod , only : r8 => shr_kind_r8, CL => SHR_KIND_CL, CS => shr_kind_CS use shr_const_mod , only : SHR_CONST_CDAY,SHR_CONST_REARTH use shr_sys_mod , only : shr_sys_abort + use ESMF , only : ESMF_VM implicit none public @@ -13,6 +14,7 @@ module mosart_vars integer :: npes ! number of processors for mosart integer :: mpicom_rof ! communicator group for mosart logical :: barrier_timers = .false. ! barrier timers + type(ESMF_VM) :: vm ! ESMF VM ! Constants integer , parameter :: iundef = -9999999 @@ -31,12 +33,12 @@ module mosart_vars integer :: nsrest = iundef ! Type of run ! Namelist variables - character(len=CL) :: frivinp ! MOSART input data file name - logical :: ice_runoff ! true => runoff is split into liquid and ice, otherwise just liquid - character(len=CS) :: decomp_option ! decomp option - character(len=CS) :: bypass_routing_option ! bypass routing model method - character(len=CS) :: qgwl_runoff_option ! method for handling qgwl runoff - integer :: budget_frq = -24 ! budget check frequency + character(len=CL) :: frivinp ! MOSART input data file name + logical :: ice_runoff ! true => runoff is split into liquid and ice, otherwise just liquid + character(len=CS) :: decomp_option ! decomp option + character(len=CS) :: bypass_routing_option ! bypass routing model method + character(len=CS) :: qgwl_runoff_option ! method for handling qgwl runoff + integer :: budget_frq = -24 ! budget check frequency ! Metadata variables used in history and restart generation character(len=CL) :: caseid = ' ' ! case id