From 6fc1de7cfc0a6d8fc3024b6a59d4a6e958ecf04a Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Mon, 27 May 2024 02:43:09 -0600 Subject: [PATCH 01/18] first set of changes to have glc->rof coupling --- src/cpl/nuopc/rof_import_export.F90 | 10 ++++++++++ src/riverroute/mosart_control_type.F90 | 10 ++++++++-- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 9cb67db..4c84346 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -109,6 +109,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_liq') ! liq runoff from glc + call fldlist_add(fldsToRof_num, fldsToRof, 'Fgrg_ice') ! ice runoff from glc do n = 1,fldsToRof_num call NUOPC_Advertise(importState, standardName=fldsToRof(n)%stdname, & @@ -288,6 +290,14 @@ subroutine import_fields( gcomp, begr, endr, rc ) ctl%qsub(begr:endr, nfrz) = 0.0_r8 ctl%qgwl(begr:endr, nfrz) = 0.0_r8 + call state_getimport(importState, 'Fgrg_liq', 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_ice', begr, endr, ctl%area, output=ctl%qglc_ice(:), & + do_area_correction=.true., rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + end subroutine import_fields !==================================================================================== diff --git a/src/riverroute/mosart_control_type.F90 b/src/riverroute/mosart_control_type.F90 index 9d98de3..ddb17fd 100644 --- a/src/riverroute/mosart_control_type.F90 +++ b/src/riverroute/mosart_control_type.F90 @@ -49,12 +49,14 @@ 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 :: qirrig_actual(:) => null() ! minimum of irrigation and available main channel storage [m3/s] ! storage, runoff @@ -303,6 +305,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), & @@ -332,6 +336,8 @@ 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 @@ -1176,7 +1182,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) From f0fe9ac5b4da823a4945a49c227f05b81a7b133c Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 29 May 2024 02:14:28 -0600 Subject: [PATCH 02/18] more updates for cism->mosart --- src/cpl/nuopc/rof_import_export.F90 | 8 +- src/riverroute/mosart_budget_type.F90 | 131 +++++++++++++------------ src/riverroute/mosart_driver.F90 | 132 +++++++++++++++----------- src/riverroute/mosart_histfile.F90 | 7 -- src/riverroute/mosart_histflds.F90 | 14 +++ 5 files changed, 169 insertions(+), 123 deletions(-) diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 4c84346..1846f08 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -109,8 +109,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_liq') ! liq runoff from glc - call fldlist_add(fldsToRof_num, fldsToRof, 'Fgrg_ice') ! ice runoff from glc + 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, & @@ -290,11 +290,11 @@ subroutine import_fields( gcomp, begr, endr, rc ) ctl%qsub(begr:endr, nfrz) = 0.0_r8 ctl%qgwl(begr:endr, nfrz) = 0.0_r8 - call state_getimport(importState, 'Fgrg_liq', begr, endr, ctl%area, output=ctl%qglc_liq(:), & + 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_ice', begr, endr, ctl%area, output=ctl%qglc_ice(:), & + 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 diff --git a/src/riverroute/mosart_budget_type.F90 b/src/riverroute/mosart_budget_type.F90 index 258155f..f25c21f 100644 --- a/src/riverroute/mosart_budget_type.F90 +++ b/src/riverroute/mosart_budget_type.F90 @@ -3,11 +3,11 @@ 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 @@ -16,27 +16,29 @@ module mosart_budget_type 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 @@ -46,18 +48,19 @@ module mosart_budget_type 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 +111,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 = 1 + nt_ice = 2 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,53 +142,56 @@ 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 = 1 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 + 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 @@ -215,11 +229,10 @@ subroutine check_budget(this, begr, endr, ntracers, dt) 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_driver.F90 b/src/riverroute/mosart_driver.F90 index eacfe25..dde8d81 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -394,6 +394,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ! ! Local variables integer :: i, j, n, nr, ns, nt, n2, nf ! indices + integer :: nt_ice, nt_liq logical :: budget_check ! if budget check needs to be performed real(r8) :: volr_init ! temporary storage to compute dvolrdt integer :: yr, mon, day, ymd, tod ! time information @@ -419,6 +420,9 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) rc = ESMF_SUCCESS + nt_liq = 1 + nt_ice = 2 + !----------------------------------------------------- ! Get date info !----------------------------------------------------- @@ -468,7 +472,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) 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) @@ -484,7 +488,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 @@ -493,10 +496,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] @@ -505,7 +508,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') @@ -518,14 +521,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 @@ -536,7 +538,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 @@ -565,24 +567,47 @@ 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 + !--- direct to outlet: all frozen runoff (from lnd and glc) !----------------------------------------------------- - 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 + ! 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,cnt) = TRunoff%qsur(nr,nt_ice) + TRunoff%qsub(nr,nt_ice) + TRunoff%qgwl(nr,nt_ice) + ctl%qglc_ice(nr) + TRunoff%qsur(nr,nt_ice) = 0._r8 + TRunoff%qsub(nr,nt_ice) = 0._r8 + TRunoff%qgwl(nr,nt_ice) = 0._r8 + 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_ice) = ctl%direct(nr,nt_ice) + dst_direct(nt,cnt) + enddo + + !----------------------------------------------------- + !--- direct to outlet: all liquid runoff from glc + !----------------------------------------------------- + + 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) enddo call ESMF_FieldSMM(Tunit%srcfield, Tunit%dstfield, Tunit%rh_direct, termorderflag=ESMF_TERMORDER_SRCSEQ, rc=rc) @@ -593,17 +618,16 @@ 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 !----------------------------------------------------- - !--- 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 @@ -612,12 +636,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 @@ -629,63 +653,59 @@ 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, qglc_liq and qflc_ice !----------------------------------------------------- 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 + ! Add glc->rof liquid going directly to outlet + src_direct(nt,cnt) = src_direct(nt,cnt) + ctl%qglc_liq(nr) 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 @@ -698,11 +718,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 @@ -721,15 +743,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 @@ -818,7 +845,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 @@ -847,7 +873,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !----------------------------------- ! 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..e042f40 100644 --- a/src/riverroute/mosart_histflds.F90 +++ b/src/riverroute/mosart_histflds.F90 @@ -31,6 +31,8 @@ module mosart_histflds type(hist_pointer_type), allocatable :: h_qgwl(:) real(r8), pointer :: h_volr_mch(:) + real(r8), pointer :: h_qglc_liq(:) + real(r8), pointer :: h_qglc_ice(:) !------------------------------------------------------------------------ contains @@ -75,6 +77,8 @@ subroutine mosart_histflds_init(begr, endr, ntracers) end do allocate(h_volr_mch(begr:endr)) + allocate(h_qglc_liq(begr:endr)) + allocate(h_qglc_ice(begr:endr)) !------------------------------------------------------- ! Build master field list of all possible fields in a history file. @@ -138,6 +142,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', units='m3', & + avgflag='A', long_name='liquid runoff from glc input', & + ptr_rof=h_qglc_liq, default='inactive') + + call mosart_hist_addfld (fname='QGLC_ICE', units='m3', & + avgflag='A', long_name='ice runoff from glc input', & + ptr_rof=h_qglc_ice, default='inactive') + ! print masterlist of history fields call mosart_hist_printflds() @@ -169,6 +181,8 @@ subroutine mosart_histflds_set(ntracers) h_qgwl(nt)%data(:) = ctl%qgwl(:,nt) end do h_volr_mch(:) = Trunoff%wr(:,1) + h_qglc_liq(:) = ctl%qglc_liq(:) + h_qglc_ice(:) = ctl%qglc_ice(:) end subroutine mosart_histflds_set From 2a0f772a12b0a4cac79389f86c800625b801dc6a Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 30 May 2024 04:31:25 -0600 Subject: [PATCH 03/18] fixed ice discharge --- src/riverroute/mosart_driver.F90 | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index dde8d81..9a5c153 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -394,7 +394,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ! ! Local variables integer :: i, j, n, nr, ns, nt, n2, nf ! indices - integer :: nt_ice, nt_liq + integer :: nt_ice, nt_liq ! incices logical :: budget_check ! if budget check needs to be performed real(r8) :: volr_init ! temporary storage to compute dvolrdt integer :: yr, mon, day, ymd, tod ! time information @@ -566,6 +566,12 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) call ESMF_FieldGet(Tunit%dstfield, farrayPtr=dst_direct, rc=rc) if (chkerr(rc,__LINE__,u_FILE_u)) return + !----------------------------------------------------- + !--- initialize ctl%direct + !----------------------------------------------------- + + ctl%direct = 0._r8 + !----------------------------------------------------- !--- direct to outlet: all frozen runoff (from lnd and glc) !----------------------------------------------------- @@ -580,7 +586,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) cnt = 0 do nr = begr,endr cnt = cnt + 1 - src_direct(nt,cnt) = TRunoff%qsur(nr,nt_ice) + TRunoff%qsub(nr,nt_ice) + TRunoff%qgwl(nr,nt_ice) + ctl%qglc_ice(nr) + src_direct(nt_ice,cnt) = TRunoff%qsur(nr,nt_ice) + TRunoff%qsub(nr,nt_ice) + TRunoff%qgwl(nr,nt_ice) + ctl%qglc_ice(nr) TRunoff%qsur(nr,nt_ice) = 0._r8 TRunoff%qsub(nr,nt_ice) = 0._r8 TRunoff%qgwl(nr,nt_ice) = 0._r8 @@ -590,11 +596,10 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, 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_ice) = ctl%direct(nr,nt_ice) + dst_direct(nt,cnt) + ctl%direct(nr,nt_ice) = ctl%direct(nr,nt_ice) + dst_direct(nt_ice,cnt) enddo !----------------------------------------------------- @@ -614,7 +619,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, 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 @@ -685,7 +689,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) endif endif ! Add glc->rof liquid going directly to outlet - src_direct(nt,cnt) = src_direct(nt,cnt) + ctl%qglc_liq(nr) + src_direct(nt_liq,cnt) = src_direct(nt_liq,cnt) + ctl%qglc_liq(nr) enddo endif From ffd51738a9eb09f0635fd07574e66c942592a2b4 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Sun, 2 Jun 2024 10:25:20 -0600 Subject: [PATCH 04/18] separation of input glc runoff and output direct rof due to glc input --- src/cpl/nuopc/rof_import_export.F90 | 47 ++++++----------- src/riverroute/mosart_budget_type.F90 | 13 +++-- src/riverroute/mosart_control_type.F90 | 73 ++++++++++++++++++++------ src/riverroute/mosart_driver.F90 | 61 +++++++++++---------- src/riverroute/mosart_histflds.F90 | 36 +++++++++---- 5 files changed, 139 insertions(+), 91 deletions(-) diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 1846f08..d31cb60 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -241,7 +241,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)' !--------------------------------------------------------------------------- @@ -252,17 +252,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 @@ -279,7 +270,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 @@ -287,8 +278,8 @@ 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 call state_getimport(importState, 'Fgrg_rofl', begr, endr, ctl%area, output=ctl%qglc_liq(:), & do_area_correction=.true., rc=rc) @@ -315,7 +306,7 @@ 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) :: flood(begr:endr) @@ -335,16 +326,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 @@ -361,23 +344,27 @@ 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(n) = rofl(n) + ctl%direct_glc(n,nliq) / (ctl%area(n)*0.001_r8) + rofi(n) = rofl(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. diff --git a/src/riverroute/mosart_budget_type.F90 b/src/riverroute/mosart_budget_type.F90 index f25c21f..b15571f 100644 --- a/src/riverroute/mosart_budget_type.F90 +++ b/src/riverroute/mosart_budget_type.F90 @@ -125,8 +125,8 @@ subroutine set_budget(this, begr, endr, ntracers, dt) integer :: nt_liq, nt_ice !------------------------------------------------- - nt_liq = 1 - nt_ice = 2 + 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) @@ -180,11 +180,11 @@ subroutine check_budget(this, begr, endr, ntracers, dt) tmp_in = 0.0_r8 tmp_glob = 0.0_r8 - nt_liq = 1 + 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) + 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 @@ -195,9 +195,8 @@ subroutine check_budget(this, begr, endr, ntracers, dt) 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 diff --git a/src/riverroute/mosart_control_type.F90 b/src/riverroute/mosart_control_type.F90 index ddb17fd..76d6452 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, CL => SHR_KIND_CL - 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, shr_mpi_max - 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, CL => SHR_KIND_CL, 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, shr_mpi_max + 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 implicit none private @@ -27,6 +28,8 @@ 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 ! decomp info integer :: begr ! local start index @@ -56,8 +59,9 @@ module mosart_control_type ! 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 :: 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] @@ -90,6 +94,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 @@ -117,13 +122,49 @@ 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) + + ! 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 @@ -315,6 +356,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' @@ -343,6 +385,7 @@ subroutine Init(this, locfn, decomp_option, use_halo_option, IDkey, rc) 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 diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 9a5c153..18fccbf 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -8,7 +8,6 @@ module mosart_driver use shr_sys_mod , only : shr_sys_abort use shr_mpi_mod , only : shr_mpi_sum, shr_mpi_max 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, & @@ -63,6 +62,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__ !----------------------------------------------------------------------- @@ -149,12 +150,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' @@ -394,7 +393,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ! ! Local variables integer :: i, j, n, nr, ns, nt, n2, nf ! indices - integer :: nt_ice, nt_liq ! incices logical :: budget_check ! if budget check needs to be performed real(r8) :: volr_init ! temporary storage to compute dvolrdt integer :: yr, mon, day, ymd, tod ! time information @@ -420,9 +418,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) rc = ESMF_SUCCESS - nt_liq = 1 - nt_ice = 2 - !----------------------------------------------------- ! Get date info !----------------------------------------------------- @@ -570,26 +565,21 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !--- initialize ctl%direct !----------------------------------------------------- - ctl%direct = 0._r8 + ctl%direct(:,:) = 0._r8 + ctl%direct_glc(:,:) = 0._r8 !----------------------------------------------------- - !--- direct to outlet: all frozen runoff (from lnd and glc) + !--- direct to outlet: all liquid and frozen runoff from glc !----------------------------------------------------- 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_ice,cnt) = TRunoff%qsur(nr,nt_ice) + TRunoff%qsub(nr,nt_ice) + TRunoff%qgwl(nr,nt_ice) + ctl%qglc_ice(nr) - TRunoff%qsur(nr,nt_ice) = 0._r8 - TRunoff%qsub(nr,nt_ice) = 0._r8 - TRunoff%qgwl(nr,nt_ice) = 0._r8 + 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) @@ -599,11 +589,12 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) cnt = 0 do nr = begr,endr cnt = cnt + 1 - ctl%direct(nr,nt_ice) = ctl%direct(nr,nt_ice) + dst_direct(nt_ice,cnt) + ctl%direct_glc(nr,nt_liq) = ctl%direct_glc(nr,nt_liq) + dst_direct(nt_liq,cnt) + ctl%direct_glc(nr,nt_ice) = ctl%direct_glc(nr,nt_ice) + dst_direct(nt_ice,cnt) enddo !----------------------------------------------------- - !--- direct to outlet: all liquid runoff from glc + !--- direct to outlet: all frozen runoff from lnd !----------------------------------------------------- src_direct(:,:) = 0._r8 @@ -612,7 +603,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) cnt = 0 do nr = begr,endr cnt = cnt + 1 - src_direct(nt_liq,cnt) = ctl%qglc_liq(nr) + 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) @@ -622,9 +613,18 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) cnt = 0 do nr = begr,endr cnt = cnt + 1 - ctl%direct(nr,nt_liq) = ctl%direct(nr,nt_liq) + dst_direct(nt_liq,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 !----------------------------------------------------- @@ -662,7 +662,7 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) endif !----------------------------------------------------- - !--- direct in place qgwl, qgwl, qglc_liq and qflc_ice + !--- direct in place qgwl, qgwl !----------------------------------------------------- if (trim(bypass_routing_option) == 'direct_in_place') then @@ -688,8 +688,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) TRunoff%qgwl(nr,nt_liq) = 0._r8 endif endif - ! Add glc->rof liquid going directly to outlet - src_direct(nt_liq,cnt) = src_direct(nt_liq,cnt) + ctl%qglc_liq(nr) enddo endif @@ -872,6 +870,13 @@ 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%runoff(nr,nt_liq) + ctl%direct_glc(nr,nt_liq) + ctl%runofftot(nr,nt_ice) = ctl%runoff(nr,nt_ice) + ctl%direct_glc(nr,nt_ice) + end do + call t_stopf('mosartr_subcycling') !----------------------------------- diff --git a/src/riverroute/mosart_histflds.F90 b/src/riverroute/mosart_histflds.F90 index e042f40..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,8 +32,8 @@ module mosart_histflds type(hist_pointer_type), allocatable :: h_qgwl(:) real(r8), pointer :: h_volr_mch(:) - real(r8), pointer :: h_qglc_liq(:) - real(r8), pointer :: h_qglc_ice(:) + real(r8), pointer :: h_qglc_liq_input(:) + real(r8), pointer :: h_qglc_ice_input(:) !------------------------------------------------------------------------ contains @@ -62,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)) @@ -75,10 +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(begr:endr)) - allocate(h_qglc_ice(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. @@ -95,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)), & @@ -105,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') @@ -142,13 +150,13 @@ 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', units='m3', & + call mosart_hist_addfld (fname='QGLC_LIQ_INPUT', units='m3', & avgflag='A', long_name='liquid runoff from glc input', & - ptr_rof=h_qglc_liq, default='inactive') + ptr_rof=h_qglc_liq_input, default='active') - call mosart_hist_addfld (fname='QGLC_ICE', units='m3', & + call mosart_hist_addfld (fname='QGLC_ICE_INPUT', units='m3', & avgflag='A', long_name='ice runoff from glc input', & - ptr_rof=h_qglc_ice, default='inactive') + ptr_rof=h_qglc_ice_input, default='active') ! print masterlist of history fields call mosart_hist_printflds() @@ -168,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) @@ -181,8 +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(:) = ctl%qglc_liq(:) - h_qglc_ice(:) = ctl%qglc_ice(:) + 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 From 38b64c84ffed689f1f99fe69bf4bd695dc6f7584 Mon Sep 17 00:00:00 2001 From: mvertens Date: Thu, 6 Jun 2024 04:57:21 -0600 Subject: [PATCH 05/18] updates to get build working with nag --- src/riverroute/mosart_control_type.F90 | 8 ++++---- src/riverroute/mosart_driver.F90 | 3 +-- src/riverroute/mosart_io.F90 | 4 ++-- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/riverroute/mosart_control_type.F90 b/src/riverroute/mosart_control_type.F90 index 76d6452..ec98507 100644 --- a/src/riverroute/mosart_control_type.F90 +++ b/src/riverroute/mosart_control_type.F90 @@ -1159,10 +1159,10 @@ subroutine calc_gradient(this, fld, fld_halo_array, dfld_dx, dfld_dy, rc) integer :: i, n, nr ! local indices real(r8) :: deg2rad real(r8) :: mean_dx, mean_dy, dlon, dlat - real(r8) :: ax_indices(4) ! x indices to add - real(r8) :: sx_indices(4) ! x indices to subtract - real(r8) :: ay_indices(4) ! y indices to add - real(r8) :: sy_indices(4) ! y indices to subtract + integer :: ax_indices(4) ! x indices to add + integer :: sx_indices(4) ! x indices to subtract + integer :: ay_indices(4) ! y indices to add + integer :: sy_indices(4) ! y indices to subtract real(r8) :: fld_surrounding(max_num_halo) real(r8) :: dx(max_num_halo) real(r8) :: dy(max_num_halo) diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 18fccbf..2ec8161 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -6,7 +6,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_mpi_mod , only : shr_mpi_sum, shr_mpi_max use shr_const_mod , only : SHR_CONST_PI, SHR_CONST_CDAY use mosart_vars , only : re, spval, iulog, ice_runoff, & frivinp, nsrContinue, nsrBranch, nsrStartup, nsrest, & @@ -166,7 +165,7 @@ subroutine mosart_read_namelist() 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,l)' ) ' use_halo_option = ',use_halo_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) diff --git a/src/riverroute/mosart_io.F90 b/src/riverroute/mosart_io.F90 index 15ae9f9..98d6e77 100644 --- a/src/riverroute/mosart_io.F90 +++ b/src/riverroute/mosart_io.F90 @@ -200,8 +200,8 @@ subroutine ncd_pio_openfile(file, fname, mode) character(len=*),parameter :: subname='ncd_pio_openfile' ! subroutine name !----------------------------------------------------------------------- + call pio_seterrorhandling(file, PIO_BCAST_ERROR) ierr = pio_openfile(pio_subsystem, file, io_type, fname, mode) - if(ierr/= PIO_NOERR) then call shr_sys_abort(subname//'ERROR: Failed to open file') else if(pio_iotask_rank(pio_subsystem)==0 .and. mainproc) then @@ -247,8 +247,8 @@ subroutine ncd_pio_createfile(file, fname) if(io_type == PIO_IOTYPE_NETCDF .or. io_type == PIO_IOTYPE_PNETCDF) then iomode = ior(iomode, io_format) endif + call pio_seterrorhandling(file, PIO_BCAST_ERROR) ierr = pio_createfile(pio_subsystem, file, io_type, fname, iomode) - if(ierr/= PIO_NOERR) then call shr_sys_abort( subname//' ERROR: Failed to open file to write: '//trim(fname)) else if(pio_iotask_rank(pio_subsystem)==0 .and. mainproc) then From b8f3aefdd2c7b57e6b1fb157b098414810412b8d Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Mon, 3 Jun 2024 03:33:55 -0600 Subject: [PATCH 06/18] cleaned up glc direct claculation --- src/riverroute/mosart_driver.F90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 2ec8161..917a5ca 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -465,7 +465,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) call t_stopf('mosartr_budgetset') endif - ! initialize data for euler solver, in m3/s here do nr = begr,endr do nt = 1,ntracers @@ -565,7 +564,6 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !----------------------------------------------------- ctl%direct(:,:) = 0._r8 - ctl%direct_glc(:,:) = 0._r8 !----------------------------------------------------- !--- direct to outlet: all liquid and frozen runoff from glc @@ -588,8 +586,8 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) cnt = 0 do nr = begr,endr cnt = cnt + 1 - ctl%direct_glc(nr,nt_liq) = ctl%direct_glc(nr,nt_liq) + dst_direct(nt_liq,cnt) - ctl%direct_glc(nr,nt_ice) = ctl%direct_glc(nr,nt_ice) + dst_direct(nt_ice,cnt) + ctl%direct_glc(nr,nt_liq) = dst_direct(nt_liq,cnt) + ctl%direct_glc(nr,nt_ice) = dst_direct(nt_ice,cnt) enddo !----------------------------------------------------- From d729df9db1f186360b78dcf1b15e8f9a143fd1aa Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Wed, 5 Jun 2024 22:47:00 -0600 Subject: [PATCH 07/18] new ability to set separate of fluxes to mediator as option --- cime_config/namelist_definition_mosart.xml | 16 ++++++++- src/cpl/nuopc/rof_import_export.F90 | 33 +++++++++++++++--- src/riverroute/mosart_driver.F90 | 39 ++++++++++++---------- src/riverroute/mosart_vars.F90 | 13 ++++---- 4 files changed, 72 insertions(+), 29 deletions(-) diff --git a/cime_config/namelist_definition_mosart.xml b/cime_config/namelist_definition_mosart.xml index a6ad704..24916b8 100644 --- a/cime_config/namelist_definition_mosart.xml +++ b/cime_config/namelist_definition_mosart.xml @@ -288,10 +288,24 @@ -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 + + logical + mosart + mosart_inparm + + .false. + + + Default: .false. + If .true., glc2ocn fluxes that are passed through mosart will be sent + as a separate fields to the mediator. + + + diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index d31cb60..00938a8 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -9,7 +9,7 @@ module rof_import_export use NUOPC_Model , only : NUOPC_ModelGet use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_abort - use mosart_vars , only : iulog, mainproc, mpicom_rof, ice_runoff + use mosart_vars , only : iulog, mainproc, mpicom_rof, ice_runoff, separate_glc2ocn_fluxes use mosart_data , only : ctl, TRunoff, TUnit use mosart_timemanager , only : get_nstep use nuopc_shr_methods , only : chkerr @@ -82,9 +82,14 @@ 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') + if (separate_glc2ocn_fluxes) then + call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofl_glc') + call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofi_glc') + end if call fldlist_add(fldsFrRof_num, fldsFrRof, 'Flrr_flood') call fldlist_add(fldsFrRof_num, fldsFrRof, 'Flrr_volr') call fldlist_add(fldsFrRof_num, fldsFrRof, 'Flrr_volrmch') @@ -309,6 +314,8 @@ subroutine export_fields (gcomp, begr, endr, rc) 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) @@ -361,10 +368,18 @@ subroutine export_fields (gcomp, begr, endr, rc) rofi(n) = 0._r8 end do end if - do n = begr,endr - rofl(n) = rofl(n) + ctl%direct_glc(n,nliq) / (ctl%area(n)*0.001_r8) - rofi(n) = rofl(n) + ctl%direct_glc(n,nice) / (ctl%area(n)*0.001_r8) - end do + + if (separate_glc2ocn_fluxes) then + 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 + else + do n = begr,endr + rofl(n) = rofl(n) + ctl%direct_glc(n,nliq) / (ctl%area(n)*0.001_r8) + rofi(n) = rofi(n) + ctl%direct_glc(n,nice) / (ctl%area(n)*0.001_r8) + end do + end if ! 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. @@ -388,6 +403,14 @@ 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 + if (separate_glc2ocn_fluxes) then + 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 + end if + call state_setexport(exportState, 'Flrr_flood', begr, endr, input=flood, do_area_correction=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 917a5ca..9bdcb9b 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -11,7 +11,8 @@ module mosart_driver frivinp, nsrContinue, nsrBranch, nsrStartup, nsrest, & inst_index, inst_suffix, inst_name, decomp_option, & bypass_routing_option, qgwl_runoff_option, barrier_timers, & - mainproc, npes, iam, mpicom_rof, budget_frq, isecspday + mainproc, npes, iam, mpicom_rof, budget_frq, isecspday, & + separate_glc2ocn_fluxes use mosart_data , only : ctl, Tctl, Tunit, TRunoff, Tpara use mosart_budget_type , only : budget_type use mosart_fileutils , only : getfil @@ -42,11 +43,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 @@ -91,7 +92,8 @@ subroutine mosart_read_namelist() namelist /mosart_inparm / frivinp, finidat, nrevsn, coupling_period, ice_runoff, & ndens, mfilt, nhtfrq, fincl1, fincl2, fincl3, fexcl1, fexcl2, fexcl3, & avgflag_pertape, decomp_option, bypass_routing_option, qgwl_runoff_option, & - use_halo_option, delt_mosart, mosart_tracers, mosart_euler_calc, budget_frq + use_halo_option, delt_mosart, mosart_tracers, mosart_euler_calc, budget_frq, & + separate_glc2ocn_fluxes ! Preset values ice_runoff = .true. @@ -105,6 +107,7 @@ subroutine mosart_read_namelist() use_halo_option = .false. mosart_tracers = 'LIQ:ICE' mosart_euler_calc = 'T:F' + separate_glc2ocn_fluxes = .false. nlfilename_rof = "mosart_in" // trim(inst_suffix) inquire (file = trim(nlfilename_rof), exist = lexist) @@ -148,6 +151,7 @@ subroutine mosart_read_namelist() call mpi_bcast (mosart_tracers, CS, MPI_CHARACTER, 0, mpicom_rof, ier) 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) + call mpi_bcast (separate_glc2ocn_fluxes, 1, MPI_LOGICAL, 0, mpicom_rof, ier) ! Determine number of tracers and array of tracer names and initialize module variables call ctl%init_tracer_names(mosart_tracers) @@ -161,17 +165,18 @@ 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) + write(iulog,'(a,l1)') ' separate_glc2ocn_fluxes = ',separate_glc2ocn_fluxes 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 diff --git a/src/riverroute/mosart_vars.F90 b/src/riverroute/mosart_vars.F90 index cc0cc48..efb40f0 100644 --- a/src/riverroute/mosart_vars.F90 +++ b/src/riverroute/mosart_vars.F90 @@ -31,12 +31,13 @@ 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 + logical :: separate_glc2cn_fluxes ! true => send fluxes from glc through mozart separately to mediator ! Metadata variables used in history and restart generation character(len=CL) :: caseid = ' ' ! case id From 410b04a6998dc9ca88b4143d58f8efa65e9535fd Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 6 Jun 2024 08:21:02 -0600 Subject: [PATCH 08/18] updates for pio error checking --- src/riverroute/mosart_io.F90 | 33 +++++++++++++++++++++------------ src/riverroute/mosart_vars.F90 | 2 +- 2 files changed, 22 insertions(+), 13 deletions(-) diff --git a/src/riverroute/mosart_io.F90 b/src/riverroute/mosart_io.F90 index 98d6e77..06945ea 100644 --- a/src/riverroute/mosart_io.F90 +++ b/src/riverroute/mosart_io.F90 @@ -196,17 +196,19 @@ subroutine ncd_pio_openfile(file, fname, mode) integer , intent(in) :: mode ! file mode ! Local variables + integer :: oldmethod integer :: ierr character(len=*),parameter :: subname='ncd_pio_openfile' ! subroutine name !----------------------------------------------------------------------- - call pio_seterrorhandling(file, PIO_BCAST_ERROR) + call pio_seterrorhandling(pio_subsystem, PIO_BCAST_ERROR, oldmethod) ierr = pio_openfile(pio_subsystem, file, io_type, fname, mode) if(ierr/= PIO_NOERR) then call shr_sys_abort(subname//'ERROR: Failed to open file') else if(pio_iotask_rank(pio_subsystem)==0 .and. mainproc) then write(iulog,*) 'Opened existing file ', trim(fname), file%fh end if + call pio_seterrorhandling(pio_subsystem, oldmethod) end subroutine ncd_pio_openfile @@ -237,6 +239,7 @@ subroutine ncd_pio_createfile(file, fname) character(len=*), intent(in) :: fname ! File name to create ! Local variables + integer :: oldmethod integer :: ierr integer :: iomode character(len=*),parameter :: subname='ncd_pio_createfile' ! subroutine name @@ -247,13 +250,14 @@ subroutine ncd_pio_createfile(file, fname) if(io_type == PIO_IOTYPE_NETCDF .or. io_type == PIO_IOTYPE_PNETCDF) then iomode = ior(iomode, io_format) endif - call pio_seterrorhandling(file, PIO_BCAST_ERROR) + call pio_seterrorhandling(pio_subsystem, PIO_BCAST_ERROR, oldmethod) ierr = pio_createfile(pio_subsystem, file, io_type, fname, iomode) if(ierr/= PIO_NOERR) then call shr_sys_abort( subname//' ERROR: Failed to open file to write: '//trim(fname)) else if(pio_iotask_rank(pio_subsystem)==0 .and. mainproc) then write(iulog,*) 'Opened file ', trim(fname), ' to write', file%fh end if + call pio_seterrorhandling(pio_subsystem, oldmethod) end subroutine ncd_pio_createfile @@ -272,6 +276,7 @@ subroutine check_var(ncid, varname, vardesc, readvar, print_err ) logical, optional, intent(in) :: print_err ! If should print about error ! Local variables + integer :: oldmethod integer :: ret ! return value logical :: log_err ! if should log error character(len=*),parameter :: subname='check_var' ! subroutine name @@ -284,14 +289,15 @@ subroutine check_var(ncid, varname, vardesc, readvar, print_err ) log_err = .true. end if readvar = .true. - call pio_seterrorhandling(ncid, PIO_BCAST_ERROR) + call pio_seterrorhandling(ncid, PIO_BCAST_ERROR, oldmethod) ret = pio_inq_varid (ncid, varname, vardesc) if (ret /= PIO_NOERR) then readvar = .false. - if (mainproc .and. log_err) & - write(iulog,*) subname//': variable ',trim(varname),' is not on dataset' + if (mainproc .and. log_err) then + write(iulog,*) subname//': variable ',trim(varname),' is not on dataset' + end if end if - call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR) + call pio_seterrorhandling(ncid, oldmethod) end subroutine check_var @@ -355,11 +361,12 @@ subroutine ncd_inqdid(ncid,name,dimid,dimexist) logical,optional, intent(out):: dimexist ! if this dimension exists or not ! Local variables + integer :: oldmethod integer :: status !----------------------------------------------------------------------- if ( present(dimexist) )then - call pio_seterrorhandling(ncid, PIO_BCAST_ERROR) + call pio_seterrorhandling(ncid, PIO_BCAST_ERROR, oldmethod) end if status = PIO_inq_dimid(ncid,name,dimid) if ( present(dimexist) )then @@ -368,7 +375,7 @@ subroutine ncd_inqdid(ncid,name,dimid,dimexist) else dimexist = .false. end if - call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR) + call pio_seterrorhandling(ncid, oldmethod) end if end subroutine ncd_inqdid @@ -430,6 +437,7 @@ subroutine ncd_inqfdims(ncid, isgrid2d, ni, nj, ns) integer , intent(out) :: nj integer , intent(out) :: ns ! Local variables + integer :: oldmethod integer :: dimid ! netCDF id integer :: ier ! error status character(len=CS) :: subname = 'surfrd_filedims' ! subroutine name @@ -438,7 +446,7 @@ subroutine ncd_inqfdims(ncid, isgrid2d, ni, nj, ns) ni = 0 nj = 0 - call pio_seterrorhandling(ncid, PIO_BCAST_ERROR) + call pio_seterrorhandling(ncid, PIO_BCAST_ERROR, oldmethod) ier = pio_inq_dimid (ncid, 'lon', dimid) if (ier == PIO_NOERR) ier = pio_inq_dimlen(ncid, dimid, ni) ier = pio_inq_dimid (ncid, 'lat', dimid) @@ -460,7 +468,7 @@ subroutine ncd_inqfdims(ncid, isgrid2d, ni, nj, ns) if (ier == PIO_NOERR) nj = 1 end if - call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR) + call pio_seterrorhandling(ncid, oldmethod) if (ni == 0 .or. nj == 0) then write(iulog,*) trim(subname),' ERROR: ni,nj = ',ni,nj,' cannot be zero ' @@ -492,13 +500,14 @@ subroutine ncd_inqvid(ncid,name,varid,vardesc,readvar) logical, optional, intent(out) :: readvar ! does variable exist ! Local variables + integer :: oldmethod integer :: ret ! return code character(len=*),parameter :: subname='ncd_inqvid' ! subroutine name !----------------------------------------------------------------------- if (present(readvar)) then readvar = .false. - call pio_seterrorhandling(ncid, PIO_BCAST_ERROR) + call pio_seterrorhandling(pio_subsystem, PIO_BCAST_ERROR, oldmethod) ret = pio_inq_varid(ncid,name,vardesc) if (ret /= PIO_NOERR) then if (mainproc) write(iulog,*) subname//': variable ',trim(name),' is not on dataset' @@ -506,7 +515,7 @@ subroutine ncd_inqvid(ncid,name,varid,vardesc,readvar) else readvar = .true. end if - call pio_seterrorhandling(ncid, PIO_INTERNAL_ERROR) + call pio_seterrorhandling(ncid, oldmethod) else ret = pio_inq_varid(ncid,name,vardesc) endif diff --git a/src/riverroute/mosart_vars.F90 b/src/riverroute/mosart_vars.F90 index efb40f0..114fa2e 100644 --- a/src/riverroute/mosart_vars.F90 +++ b/src/riverroute/mosart_vars.F90 @@ -37,7 +37,7 @@ module mosart_vars 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 - logical :: separate_glc2cn_fluxes ! true => send fluxes from glc through mozart separately to mediator + logical :: separate_glc2ocn_fluxes ! true => send fluxes from glc through mozart separately to mediator ! Metadata variables used in history and restart generation character(len=CL) :: caseid = ' ' ! case id From c95417e6b12a8f5a25e0e47c7c388a0d0bf725ab Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Fri, 7 Jun 2024 05:59:42 -0600 Subject: [PATCH 09/18] addressed issue https://github.com/ESCOMP/MOSART/issues/95 --- src/riverroute/mosart_budget_type.F90 | 32 ++++++++++++++++----------- src/riverroute/mosart_timemanager.F90 | 2 +- 2 files changed, 20 insertions(+), 14 deletions(-) diff --git a/src/riverroute/mosart_budget_type.F90 b/src/riverroute/mosart_budget_type.F90 index b15571f..3df5c8a 100644 --- a/src/riverroute/mosart_budget_type.F90 +++ b/src/riverroute/mosart_budget_type.F90 @@ -9,7 +9,6 @@ module mosart_budget_type 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 @@ -46,6 +45,13 @@ 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__ @@ -201,12 +207,12 @@ subroutine check_budget(this, begr, endr, ntracers, dt) 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 @@ -216,12 +222,12 @@ 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. 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, & From 9013f77683054005aa0a3912c8a98a8212a973e0 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Fri, 14 Jun 2024 09:52:29 -0600 Subject: [PATCH 10/18] various updates --- src/cpl/nuopc/rof_comp_nuopc.F90 | 8 ++--- src/cpl/nuopc/rof_import_export.F90 | 45 +++++++++++++----------- src/riverroute/mosart_control_type.F90 | 18 ++++++---- src/riverroute/mosart_physics.F90 | 48 +++++--------------------- src/riverroute/mosart_vars.F90 | 2 ++ 5 files changed, 52 insertions(+), 69 deletions(-) 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 00938a8..85650eb 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -9,7 +9,7 @@ module rof_import_export use NUOPC_Model , only : NUOPC_ModelGet use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_abort - use mosart_vars , only : iulog, mainproc, mpicom_rof, ice_runoff, separate_glc2ocn_fluxes + use mosart_vars , only : iulog, mainproc, ice_runoff, separate_glc2ocn_fluxes, vm use mosart_data , only : ctl, TRunoff, TUnit use mosart_timemanager , only : get_nstep use nuopc_shr_methods , only : chkerr @@ -131,8 +131,8 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) use ESMF , only : ESMF_GridComp, ESMF_StateGet use ESMF , only : ESMF_Mesh, ESMF_MeshGet use ESMF , only : ESMF_Field, ESMF_FieldGet, ESMF_FieldRegridGetArea + use ESMF , only : ESMF_GridCompGet, ESMF_VMAllReduce, ESMF_REDUCE_MAX, ESMF_REDUCE_MIN use shr_const_mod , only : shr_const_rearth - use shr_mpi_mod , only : shr_mpi_min, shr_mpi_max ! input/output variables type(ESMF_GridComp) , intent(inout) :: gcomp @@ -151,14 +151,14 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) real(r8), allocatable :: model_areas(:) real(r8), pointer :: dataptr(:) real(r8) :: re = shr_const_rearth*0.001_r8 ! radius of earth (km) - real(r8) :: max_mod2med_areacor - real(r8) :: max_med2mod_areacor - real(r8) :: min_mod2med_areacor - real(r8) :: min_med2mod_areacor - real(r8) :: max_mod2med_areacor_glob - real(r8) :: max_med2mod_areacor_glob - real(r8) :: min_mod2med_areacor_glob - real(r8) :: min_med2mod_areacor_glob + real(r8) :: max_mod2med_areacor(1) + real(r8) :: max_med2mod_areacor(1) + real(r8) :: min_mod2med_areacor(1) + real(r8) :: min_med2mod_areacor(1) + real(r8) :: max_mod2med_areacor_glob(1) + real(r8) :: max_med2mod_areacor_glob(1) + real(r8) :: min_mod2med_areacor_glob(1) + real(r8) :: min_med2mod_areacor_glob(1) character(len=*), parameter :: subname='(rof_import_export:realize_fields)' !--------------------------------------------------------------------------- @@ -213,20 +213,25 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) deallocate(model_areas) deallocate(mesh_areas) - min_mod2med_areacor = minval(mod2med_areacor) - max_mod2med_areacor = maxval(mod2med_areacor) - min_med2mod_areacor = minval(med2mod_areacor) - max_med2mod_areacor = maxval(med2mod_areacor) - call shr_mpi_max(max_mod2med_areacor, max_mod2med_areacor_glob, mpicom_rof) - call shr_mpi_min(min_mod2med_areacor, min_mod2med_areacor_glob, mpicom_rof) - call shr_mpi_max(max_med2mod_areacor, max_med2mod_areacor_glob, mpicom_rof) - call shr_mpi_min(min_med2mod_areacor, min_med2mod_areacor_glob, mpicom_rof) + min_mod2med_areacor(1) = minval(mod2med_areacor) + max_mod2med_areacor(1) = maxval(mod2med_areacor) + min_med2mod_areacor(1) = minval(med2mod_areacor) + max_med2mod_areacor(1) = maxval(med2mod_areacor) + + call ESMF_VMAllReduce(vm, max_mod2med_areacor, max_mod2med_areacor_glob, 1, ESMF_REDUCE_MAX, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllReduce(vm, min_mod2med_areacor, min_mod2med_areacor_glob, 1, ESMF_REDUCE_MIN, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllReduce(vm, max_med2mod_areacor, max_med2mod_areacor_glob, 1, ESMF_REDUCE_MAX, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllReduce(vm, min_med2mod_areacor, min_med2mod_areacor_glob, 1, ESMF_REDUCE_MIN, rc=rc) + if (chkerr(rc,__LINE__,u_FILE_u)) return if (mainproc) then write(iulog,'(2A,2g23.15,A )') trim(subname),' : min_mod2med_areacor, max_mod2med_areacor ',& - min_mod2med_areacor_glob, max_mod2med_areacor_glob, 'MOSART' + min_mod2med_areacor_glob(1), max_mod2med_areacor_glob(1), 'MOSART' write(iulog,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& - min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'MOSART' + min_med2mod_areacor_glob(1), max_med2mod_areacor_glob(1), 'MOSART' end if end subroutine realize_fields diff --git a/src/riverroute/mosart_control_type.F90 b/src/riverroute/mosart_control_type.F90 index 792751c..d5ebdf5 100644 --- a/src/riverroute/mosart_control_type.F90 +++ b/src/riverroute/mosart_control_type.F90 @@ -3,13 +3,14 @@ module mosart_control_type 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 + 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_ArrayHalo, ESMF_ArrayGet, ESMF_VMAllReduce, ESMF_REDUCE_SUM use perf_mod, only : t_startf, t_stopf use nuopc_shr_methods, only : chkerr @@ -181,7 +182,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 @@ -402,15 +404,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 diff --git a/src/riverroute/mosart_physics.F90 b/src/riverroute/mosart_physics.F90 index 4350f64..a03c61f 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,12 +279,6 @@ 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." ) end if @@ -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_vars.F90 b/src/riverroute/mosart_vars.F90 index 114fa2e..101d5d7 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 From 0336d53f28ad54c08e566d55d2723a16157184d7 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Fri, 14 Jun 2024 14:35:30 -0600 Subject: [PATCH 11/18] more simplicification for glc input runoff --- src/cpl/nuopc/rof_import_export.F90 | 47 ++++++++++++++++++++++---- src/riverroute/mosart_control_type.F90 | 1 + src/riverroute/mosart_driver.F90 | 45 +++++++++++++----------- 3 files changed, 66 insertions(+), 27 deletions(-) diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 85650eb..3731326 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 @@ -234,6 +235,12 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) min_med2mod_areacor_glob(1), max_med2mod_areacor_glob(1), '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 + end subroutine realize_fields !=============================================================================== @@ -291,13 +298,18 @@ subroutine import_fields( gcomp, begr, endr, rc ) ctl%qsub(begr:endr, nice) = 0.0_r8 ctl%qgwl(begr:endr, nice) = 0.0_r8 - 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 + 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 + write(6,*)'DEBUG: ctl%rof_from_glc = ',ctl%rof_from_glc end subroutine import_fields @@ -679,4 +691,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_control_type.F90 b/src/riverroute/mosart_control_type.F90 index d5ebdf5..04b297c 100644 --- a/src/riverroute/mosart_control_type.F90 +++ b/src/riverroute/mosart_control_type.F90 @@ -30,6 +30,7 @@ module mosart_control_type 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 diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 9bdcb9b..b161bf6 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -574,26 +574,31 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) !--- direct to outlet: all liquid and frozen runoff from glc !----------------------------------------------------- - 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 + 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 From 7c0f3fc7cfcb8cf02530e9b9fc12e8c9da7d819c Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Sat, 15 Jun 2024 02:41:29 -0600 Subject: [PATCH 12/18] removed vm updates for now from rof_import_export.F90 --- src/cpl/nuopc/rof_import_export.F90 | 49 ++++++++++++++--------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index 3731326..dc682a4 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -9,7 +9,7 @@ module rof_import_export use NUOPC_Model , only : NUOPC_ModelGet use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_abort - use mosart_vars , only : iulog, mainproc, ice_runoff, separate_glc2ocn_fluxes, vm + use mosart_vars , only : iulog, mainproc, mpicom_rof, ice_runoff, separate_glc2ocn_fluxes use mosart_data , only : ctl, TRunoff, TUnit use mosart_timemanager , only : get_nstep use nuopc_shr_methods , only : chkerr @@ -132,8 +132,8 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) use ESMF , only : ESMF_GridComp, ESMF_StateGet use ESMF , only : ESMF_Mesh, ESMF_MeshGet use ESMF , only : ESMF_Field, ESMF_FieldGet, ESMF_FieldRegridGetArea - use ESMF , only : ESMF_GridCompGet, ESMF_VMAllReduce, ESMF_REDUCE_MAX, ESMF_REDUCE_MIN use shr_const_mod , only : shr_const_rearth + use shr_mpi_mod , only : shr_mpi_min, shr_mpi_max ! input/output variables type(ESMF_GridComp) , intent(inout) :: gcomp @@ -152,14 +152,14 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) real(r8), allocatable :: model_areas(:) real(r8), pointer :: dataptr(:) real(r8) :: re = shr_const_rearth*0.001_r8 ! radius of earth (km) - real(r8) :: max_mod2med_areacor(1) - real(r8) :: max_med2mod_areacor(1) - real(r8) :: min_mod2med_areacor(1) - real(r8) :: min_med2mod_areacor(1) - real(r8) :: max_mod2med_areacor_glob(1) - real(r8) :: max_med2mod_areacor_glob(1) - real(r8) :: min_mod2med_areacor_glob(1) - real(r8) :: min_med2mod_areacor_glob(1) + real(r8) :: max_mod2med_areacor + real(r8) :: max_med2mod_areacor + real(r8) :: min_mod2med_areacor + real(r8) :: min_med2mod_areacor + real(r8) :: max_mod2med_areacor_glob + real(r8) :: max_med2mod_areacor_glob + real(r8) :: min_mod2med_areacor_glob + real(r8) :: min_med2mod_areacor_glob character(len=*), parameter :: subname='(rof_import_export:realize_fields)' !--------------------------------------------------------------------------- @@ -214,25 +214,20 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) deallocate(model_areas) deallocate(mesh_areas) - min_mod2med_areacor(1) = minval(mod2med_areacor) - max_mod2med_areacor(1) = maxval(mod2med_areacor) - min_med2mod_areacor(1) = minval(med2mod_areacor) - max_med2mod_areacor(1) = maxval(med2mod_areacor) - - call ESMF_VMAllReduce(vm, max_mod2med_areacor, max_mod2med_areacor_glob, 1, ESMF_REDUCE_MAX, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllReduce(vm, min_mod2med_areacor, min_mod2med_areacor_glob, 1, ESMF_REDUCE_MIN, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllReduce(vm, max_med2mod_areacor, max_med2mod_areacor_glob, 1, ESMF_REDUCE_MAX, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return - call ESMF_VMAllReduce(vm, min_med2mod_areacor, min_med2mod_areacor_glob, 1, ESMF_REDUCE_MIN, rc=rc) - if (chkerr(rc,__LINE__,u_FILE_u)) return + min_mod2med_areacor = minval(mod2med_areacor) + max_mod2med_areacor = maxval(mod2med_areacor) + min_med2mod_areacor = minval(med2mod_areacor) + max_med2mod_areacor = maxval(med2mod_areacor) + call shr_mpi_max(max_mod2med_areacor, max_mod2med_areacor_glob, mpicom_rof) + call shr_mpi_min(min_mod2med_areacor, min_mod2med_areacor_glob, mpicom_rof) + call shr_mpi_max(max_med2mod_areacor, max_med2mod_areacor_glob, mpicom_rof) + call shr_mpi_min(min_med2mod_areacor, min_med2mod_areacor_glob, mpicom_rof) if (mainproc) then write(iulog,'(2A,2g23.15,A )') trim(subname),' : min_mod2med_areacor, max_mod2med_areacor ',& - min_mod2med_areacor_glob(1), max_mod2med_areacor_glob(1), 'MOSART' + min_mod2med_areacor_glob, max_mod2med_areacor_glob, 'MOSART' write(iulog,'(2A,2g23.15,A )') trim(subname),' : min_med2mod_areacor, max_med2mod_areacor ',& - min_med2mod_areacor_glob(1), max_med2mod_areacor_glob(1), 'MOSART' + min_med2mod_areacor_glob, max_med2mod_areacor_glob, 'MOSART' end if if (fldchk(importState, 'Fgrg_rofl') .and. fldchk(importState, 'Fgrg_rofl')) then @@ -240,6 +235,9 @@ subroutine realize_fields(gcomp, Emesh, flds_scalar_name, flds_scalar_num, rc) 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 @@ -309,7 +307,6 @@ subroutine import_fields( gcomp, begr, endr, rc ) ctl%qglc_liq(:) = 0._r8 ctl%qglc_ice(:) = 0._r8 end if - write(6,*)'DEBUG: ctl%rof_from_glc = ',ctl%rof_from_glc end subroutine import_fields From 4a80be551254cd98adbc71cda58e9832f66c3dee Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Sat, 15 Jun 2024 03:30:09 -0600 Subject: [PATCH 13/18] fixed mosart history for total runoff --- src/riverroute/mosart_driver.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index b161bf6..1a370b6 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -880,8 +880,8 @@ subroutine mosart_run(begr, endr, ntracers, rstwr, nlend, rdate, rc) ! final update from glc input do nr = begr,endr - ctl%runofftot(nr,nt_liq) = ctl%runoff(nr,nt_liq) + ctl%direct_glc(nr,nt_liq) - ctl%runofftot(nr,nt_ice) = ctl%runoff(nr,nt_ice) + ctl%direct_glc(nr,nt_ice) + 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') From e6fd7706b998f629c3360877f2fe0e7b2e9e7ccb Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Tue, 18 Jun 2024 03:45:28 -0600 Subject: [PATCH 14/18] made separate_glc2ocn_fluxes be true by default --- cime_config/namelist_definition_mosart.xml | 4 ++-- src/riverroute/mosart_driver.F90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/cime_config/namelist_definition_mosart.xml b/cime_config/namelist_definition_mosart.xml index 24916b8..adf10b2 100644 --- a/cime_config/namelist_definition_mosart.xml +++ b/cime_config/namelist_definition_mosart.xml @@ -299,10 +299,10 @@ mosart mosart_inparm - .false. + .true. - Default: .false. + Default: .true. If .true., glc2ocn fluxes that are passed through mosart will be sent as a separate fields to the mediator. diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index 1a370b6..bc14c52 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -107,7 +107,7 @@ subroutine mosart_read_namelist() use_halo_option = .false. mosart_tracers = 'LIQ:ICE' mosart_euler_calc = 'T:F' - separate_glc2ocn_fluxes = .false. + separate_glc2ocn_fluxes = .true. nlfilename_rof = "mosart_in" // trim(inst_suffix) inquire (file = trim(nlfilename_rof), exist = lexist) From 68995dd5ea90303ebcd0fa1ef5743326520479b8 Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 20 Jun 2024 03:02:32 -0600 Subject: [PATCH 15/18] removed namelist variable separate_glc2ocn_fluxes - mosart will only send glc fluxes separately to the mediator --- cime_config/namelist_definition_mosart.xml | 14 --------- src/cpl/nuopc/rof_import_export.F90 | 33 ++++++++-------------- src/riverroute/mosart_driver.F90 | 9 ++---- src/riverroute/mosart_physics.F90 | 2 +- src/riverroute/mosart_tctl_type.F90 | 2 +- src/riverroute/mosart_vars.F90 | 1 - 6 files changed, 15 insertions(+), 46 deletions(-) diff --git a/cime_config/namelist_definition_mosart.xml b/cime_config/namelist_definition_mosart.xml index adf10b2..86f5224 100644 --- a/cime_config/namelist_definition_mosart.xml +++ b/cime_config/namelist_definition_mosart.xml @@ -294,18 +294,4 @@ - - logical - mosart - mosart_inparm - - .true. - - - Default: .true. - If .true., glc2ocn fluxes that are passed through mosart will be sent - as a separate fields to the mediator. - - - diff --git a/src/cpl/nuopc/rof_import_export.F90 b/src/cpl/nuopc/rof_import_export.F90 index dc682a4..80cc7c7 100644 --- a/src/cpl/nuopc/rof_import_export.F90 +++ b/src/cpl/nuopc/rof_import_export.F90 @@ -9,7 +9,7 @@ module rof_import_export use NUOPC_Model , only : NUOPC_ModelGet use shr_kind_mod , only : r8 => shr_kind_r8 use shr_sys_mod , only : shr_sys_abort - use mosart_vars , only : iulog, mainproc, mpicom_rof, ice_runoff, separate_glc2ocn_fluxes + use mosart_vars , only : iulog, mainproc, mpicom_rof, ice_runoff use mosart_data , only : ctl, TRunoff, TUnit use mosart_timemanager , only : get_nstep use nuopc_shr_methods , only : chkerr @@ -87,10 +87,8 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) 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') - if (separate_glc2ocn_fluxes) then - call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofl_glc') - call fldlist_add(fldsFrRof_num, fldsFrRof, 'Forr_rofi_glc') - end if + 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') @@ -383,17 +381,10 @@ subroutine export_fields (gcomp, begr, endr, rc) end do end if - if (separate_glc2ocn_fluxes) then - 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 - else - do n = begr,endr - rofl(n) = rofl(n) + ctl%direct_glc(n,nliq) / (ctl%area(n)*0.001_r8) - rofi(n) = rofi(n) + ctl%direct_glc(n,nice) / (ctl%area(n)*0.001_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. @@ -417,13 +408,11 @@ 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 - if (separate_glc2ocn_fluxes) then - 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_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 - end if + 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 diff --git a/src/riverroute/mosart_driver.F90 b/src/riverroute/mosart_driver.F90 index bc14c52..652bf98 100644 --- a/src/riverroute/mosart_driver.F90 +++ b/src/riverroute/mosart_driver.F90 @@ -11,8 +11,7 @@ module mosart_driver frivinp, nsrContinue, nsrBranch, nsrStartup, nsrest, & inst_index, inst_suffix, inst_name, decomp_option, & bypass_routing_option, qgwl_runoff_option, barrier_timers, & - mainproc, npes, iam, mpicom_rof, budget_frq, isecspday, & - separate_glc2ocn_fluxes + mainproc, npes, iam, mpicom_rof, budget_frq, isecspday use mosart_data , only : ctl, Tctl, Tunit, TRunoff, Tpara use mosart_budget_type , only : budget_type use mosart_fileutils , only : getfil @@ -92,8 +91,7 @@ subroutine mosart_read_namelist() namelist /mosart_inparm / frivinp, finidat, nrevsn, coupling_period, ice_runoff, & ndens, mfilt, nhtfrq, fincl1, fincl2, fincl3, fexcl1, fexcl2, fexcl3, & avgflag_pertape, decomp_option, bypass_routing_option, qgwl_runoff_option, & - use_halo_option, delt_mosart, mosart_tracers, mosart_euler_calc, budget_frq, & - separate_glc2ocn_fluxes + use_halo_option, delt_mosart, mosart_tracers, mosart_euler_calc, budget_frq ! Preset values ice_runoff = .true. @@ -107,7 +105,6 @@ subroutine mosart_read_namelist() use_halo_option = .false. mosart_tracers = 'LIQ:ICE' mosart_euler_calc = 'T:F' - separate_glc2ocn_fluxes = .true. nlfilename_rof = "mosart_in" // trim(inst_suffix) inquire (file = trim(nlfilename_rof), exist = lexist) @@ -151,7 +148,6 @@ subroutine mosart_read_namelist() call mpi_bcast (mosart_tracers, CS, MPI_CHARACTER, 0, mpicom_rof, ier) 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) - call mpi_bcast (separate_glc2ocn_fluxes, 1, MPI_LOGICAL, 0, mpicom_rof, ier) ! Determine number of tracers and array of tracer names and initialize module variables call ctl%init_tracer_names(mosart_tracers) @@ -174,7 +170,6 @@ subroutine mosart_read_namelist() 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,l1)') ' separate_glc2ocn_fluxes = ',separate_glc2ocn_fluxes if (nsrest == nsrStartup .and. finidat /= ' ') then write(iulog,'(a)') ' mosart initial data = '//trim(finidat) end if diff --git a/src/riverroute/mosart_physics.F90 b/src/riverroute/mosart_physics.F90 index a03c61f..3700d42 100644 --- a/src/riverroute/mosart_physics.F90 +++ b/src/riverroute/mosart_physics.F90 @@ -280,7 +280,7 @@ subroutine mainchannelRouting(nr, nt, theDeltaT) if(Tctl%RoutingMethod == 1) then call Routing_KW(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 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_vars.F90 b/src/riverroute/mosart_vars.F90 index 101d5d7..6712c4d 100644 --- a/src/riverroute/mosart_vars.F90 +++ b/src/riverroute/mosart_vars.F90 @@ -39,7 +39,6 @@ module mosart_vars 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 - logical :: separate_glc2ocn_fluxes ! true => send fluxes from glc through mozart separately to mediator ! Metadata variables used in history and restart generation character(len=CL) :: caseid = ' ' ! case id From 01f40dcf589d4e6c5932c76128e3d42faea91e9f Mon Sep 17 00:00:00 2001 From: Mariana Vertenstein Date: Thu, 20 Jun 2024 03:10:41 -0600 Subject: [PATCH 16/18] updates for issues raised in PR review --- src/riverroute/mosart_control_type.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/riverroute/mosart_control_type.F90 b/src/riverroute/mosart_control_type.F90 index 04b297c..672613f 100644 --- a/src/riverroute/mosart_control_type.F90 +++ b/src/riverroute/mosart_control_type.F90 @@ -134,6 +134,9 @@ module mosart_control_type 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 From ca0acc0336365456bec653c4c233abb57d9bc62d Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Fri, 21 Jun 2024 16:27:38 -0600 Subject: [PATCH 17/18] Add I1850...G (ie cism) test to the mosart test-suite --- cime_config/testdefs/testlist_mosart.xml | 9 +++++++++ 1 file changed, 9 insertions(+) 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 @@ + + + + + + + + + From 19966c36fa79622d81d841033e223d5644363616 Mon Sep 17 00:00:00 2001 From: Samuel Levis Date: Fri, 21 Jun 2024 17:11:59 -0600 Subject: [PATCH 18/18] Updated ChangeLog --- docs/ChangeLog | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) 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