From 8e1f68d29f155b51b4b169a14bbeba65b8a913c2 Mon Sep 17 00:00:00 2001 From: Katetc Date: Thu, 29 Aug 2024 13:51:49 -0600 Subject: [PATCH 1/8] Fixed usrf bug and first draft simple hole filling for neg qice --- dglc/dglc_datamode_noevolve_mod.F90 | 58 ++++++++++++++++++++++++----- 1 file changed, 48 insertions(+), 10 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index d654e6aa..9d4d5cb7 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -49,8 +49,6 @@ module dglc_datamode_noevolve_mod ! type(icesheet_ptr_t), allocatable :: So_t(:) ! type(icesheet_ptr_t), allocatable :: So_q(:) - real(r8), allocatable :: usrf(:) ! upper surface elevation (m) on ice grid - ! Export Field names character(len=*), parameter :: field_out_area = 'Sg_area' character(len=*), parameter :: field_out_topo = 'Sg_topo' @@ -206,7 +204,7 @@ end subroutine dglc_datamode_noevolve_init_pointers !=============================================================================== subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & - model_meshes, model_internal_gridsize, model_datafiles, rc) + logunit, model_meshes, model_internal_gridsize, model_datafiles, rc) ! Assume that the model mesh is the same as the input data mesh @@ -214,6 +212,7 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & type(iosystem_desc_t) , pointer :: pio_subsystem ! pio info integer , intent(in) :: io_type ! pio info integer , intent(in) :: io_format ! pio info + integer , intent(in) :: logunit ! For writing logs type(ESMF_Mesh) , intent(in) :: model_meshes(:) ! ice sheets meshes real(r8) , intent(in) :: model_internal_gridsize(:) ! internal model gridsizes (m) character(len=*) , intent(in) :: model_datafiles(:) ! input file names @@ -240,7 +239,11 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & real(r8) :: rhoi ! density of ice ~ kg/m^3 real(r8) :: rhoo ! density of sea water ~ kg/m^3 real(r8) :: eus ! eustatic sea level + real(r8), allocatable :: Tot_smb(:) ! Sum smb on each ice sheet to reduce neg fluxes + real(r8) :: num_Tot ! Number of active grid cells in total calculation + real(r8) :: ice_out ! Scaled value of ice mass balance sent to rof real(r8), allocatable :: lsrf(:) ! lower surface elevation (m) on ice grid + real(r8), allocatable :: usrf(:) ! upper surface elevation (m) on ice grid character(len=*), parameter :: subname='(dglc_datamode_noevolve_advance): ' !------------------------------------------------------------------------------- @@ -361,7 +364,7 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & end do deallocate(lsrf) - + deallocate(usrf) call pio_closefile(pioid) call pio_freedecomp(pio_subsystem, pio_iodesc) @@ -369,23 +372,58 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & end if + ! Set initialized flag + initialized_noevolve = .true. + if (initialized_noevolve) then + ! Compute Fgrg_rofi do ns = 1,num_icesheets - lsize = size(Fgrg_rofi(ns)%ptr) + + ! Get mesh info + call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + + allocate(Tot_smb(lsize)) + Tot_smb(:) = 0.d0 + num_Tot = 0.d0 + ice_out = 0.d0 + + ! For No Evolve to reduce negative ice fluxes from DGLC, we will + ! Sum over all smb (Flgl_qice) for each ice sheet (not globally). + ! The output to the river model will be an equal division of the + ! final total (pos+Neg). do ng = 1,lsize - if (is_in_active_grid(usrf(ng))) then - Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + Tot_smb(ns) = Tot_smb(ns)+Flgl_qice(ns)%ptr(ng) + num_Tot = num_Tot+1.d0 + end if + end do + ! Now redistribute the scaled values + if (num_Tot.gt.0.d0) then + ice_out = Tot_smb(ns)/num_Tot + else + ice_out = 0.d0 + end if + do ng = 1,lsize + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + Fgrg_rofi(ns)%ptr(ng) = ice_out else Fgrg_rofi(ns)%ptr(ng) = 0._r8 end if end do + + deallocate(Tot_smb) + !if (is_in_active_grid(usrf(ng))) then + ! Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + !else + ! Fgrg_rofi(ns)%ptr(ng) = 0._r8 + !end if end do end if - ! Set initialized flag - initialized_noevolve = .true. - end subroutine dglc_datamode_noevolve_advance !=============================================================================== From da714ce8a7c7da64c9523818f848df967a2299f4 Mon Sep 17 00:00:00 2001 From: Katetc Date: Thu, 29 Aug 2024 14:44:00 -0600 Subject: [PATCH 2/8] Passing logunit to run comp subroutine --- dglc/glc_comp_nuopc.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dglc/glc_comp_nuopc.F90 b/dglc/glc_comp_nuopc.F90 index f5b6c2a3..11e66bfd 100644 --- a/dglc/glc_comp_nuopc.F90 +++ b/dglc/glc_comp_nuopc.F90 @@ -590,7 +590,7 @@ subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, valid_inp case('noevolve') if (valid_inputs) then call dglc_datamode_noevolve_advance(sdat(1)%pio_subsystem, sdat(1)%io_type, sdat(1)%io_format, & - model_meshes, model_internal_gridsize, model_datafiles, rc) + logunit, model_meshes, model_internal_gridsize, model_datafiles, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if end select From aedd7221d80c988adc8c61fb481520ad1f01fbb6 Mon Sep 17 00:00:00 2001 From: Katetc Date: Fri, 30 Aug 2024 10:07:25 -0600 Subject: [PATCH 3/8] Fix step 2, scaled removal by fraction of local positives --- dglc/dglc_datamode_noevolve_mod.F90 | 84 +++++++++++++++++++---------- 1 file changed, 56 insertions(+), 28 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 9d4d5cb7..8a8f96b9 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -239,9 +239,11 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & real(r8) :: rhoi ! density of ice ~ kg/m^3 real(r8) :: rhoo ! density of sea water ~ kg/m^3 real(r8) :: eus ! eustatic sea level - real(r8), allocatable :: Tot_smb(:) ! Sum smb on each ice sheet to reduce neg fluxes + real(r8) :: Tot_pos_smb ! Sum of positive smb values on each ice sheet for hole-filling + real(r8) :: Tot_neg_smb ! Sum of negative smb values on each ice sheet for hole-filling real(r8) :: num_Tot ! Number of active grid cells in total calculation - real(r8) :: ice_out ! Scaled value of ice mass balance sent to rof + real(r8) :: rat ! Ratio of hole-filling flux to apply + real(r8), allocatable :: ice_runoff_out(:) ! Scaled ice runoff output after holes filled real(r8), allocatable :: lsrf(:) ! lower surface elevation (m) on ice grid real(r8), allocatable :: usrf(:) ! upper surface elevation (m) on ice grid character(len=*), parameter :: subname='(dglc_datamode_noevolve_advance): ' @@ -386,41 +388,67 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - allocate(Tot_smb(lsize)) - Tot_smb(:) = 0.d0 + allocate(ice_runoff_out(lsize)) + ice_runoff_out(:) = 0.d0 + Tot_pos_smb = 0.d0 + Tot_neg_smb = 0.d0 num_Tot = 0.d0 - ice_out = 0.d0 + rat = 0.d0 ! For No Evolve to reduce negative ice fluxes from DGLC, we will - ! Sum over all smb (Flgl_qice) for each ice sheet (not globally). - ! The output to the river model will be an equal division of the - ! final total (pos+Neg). + ! Calculate the total positive and total negative fluxes on each + ! processor for each ice sheet. do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then - Tot_smb(ns) = Tot_smb(ns)+Flgl_qice(ns)%ptr(ng) - num_Tot = num_Tot+1.d0 + if(Flgl_qice(ns)%ptr(ng) > 0.d0) then + Tot_pos_smb = Tot_pos_smb+Flgl_qice(ns)%ptr(ng) + end if + ! Ignore places that are exactly 0.d0 + if(Flgl_qice(ns)%ptr(ng) < 0.d0) then + Tot_neg_smb = Tot_neg_smb+Flgl_qice(ns)%ptr(ng) + end if end if end do - ! Now redistribute the scaled values - if (num_Tot.gt.0.d0) then - ice_out = Tot_smb(ns)/num_Tot + ! If there's more positive than negative, then set all + ! negative to zero and destribute the negative flux amount + ! across the positive values, scaled by the size of the + ! positive value. This section also applies to any chunks + ! where there is no negative smb. In that case the ice + ! runoff is exactly equal to the input smb. + if(abs(Tot_pos_smb) >= abs(Tot_neg_smb)) then + do ng = 1,lsize + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + if(Flgl_qice(ns)%ptr(ng) > 0.d0) then + rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb + else if (Flgl_qice(ns)%ptr(ng) < 0.d0) then + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + else + Fgrg_rofi(ns)%ptr(ng) = 0._r8 + end if + end do else - ice_out = 0.d0 - end if - do ng = 1,lsize - if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then - Fgrg_rofi(ns)%ptr(ng) = ice_out - else - Fgrg_rofi(ns)%ptr(ng) = 0._r8 - end if - end do + ! If there's more negative than positive, set the positive to zero + ! and distribute the positive amount to the negative spaces to + ! reduce their negativity a bit. This shouldn't happen often. + ! This section of code also applies if Tot_pos_smb is zero. + do ng = 1,lsize + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + if(Flgl_qice(ns)%ptr(ng) < 0.d0) then + rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb + else if (Flgl_qice(ns)%ptr(ng) > 0.d0) then + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + else + Fgrg_rofi(ns)%ptr(ng) = 0._r8 + end if + end do + + end if ! More neg or pos smb - deallocate(Tot_smb) - !if (is_in_active_grid(usrf(ng))) then - ! Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) - !else - ! Fgrg_rofi(ns)%ptr(ng) = 0._r8 - !end if + deallocate(ice_runoff_out) end do end if From b73abc38b1e41e80a9e553be080f025fb69be66c Mon Sep 17 00:00:00 2001 From: Katetc Date: Fri, 30 Aug 2024 12:49:27 -0600 Subject: [PATCH 4/8] Third step complete, added global sums to spread calculation out over the entire ice sheet --- dglc/dglc_datamode_noevolve_mod.F90 | 49 +++++++++++++++++++---------- dglc/glc_comp_nuopc.F90 | 9 +++--- 2 files changed, 38 insertions(+), 20 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 8a8f96b9..068e1b81 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -4,6 +4,8 @@ module dglc_datamode_noevolve_mod use ESMF , only : ESMF_Mesh, ESMF_DistGrid, ESMF_FieldBundle, ESMF_Field use ESMF , only : ESMF_FieldBundleCreate, ESMF_FieldCreate, ESMF_MeshLoc_Element use ESMF , only : ESMF_FieldBundleAdd, ESMF_MeshGet, ESMF_DistGridGet, ESMF_Typekind_R8 + use ESMF , only : ESMF_GridComp, ESMF_GridCompGet + use ESMF , only : ESMF_VM, ESMF_VMAllreduce, ESMF_REDUCE_SUM use NUOPC , only : NUOPC_Advertise, NUOPC_IsConnected use shr_kind_mod , only : r8=>shr_kind_r8, i8=>shr_kind_i8, cl=>shr_kind_cl, cs=>shr_kind_cs use shr_sys_mod , only : shr_sys_abort @@ -203,12 +205,13 @@ subroutine dglc_datamode_noevolve_init_pointers(NStateExp, NstateImp, rc) end subroutine dglc_datamode_noevolve_init_pointers !=============================================================================== - subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & + subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_format, & logunit, model_meshes, model_internal_gridsize, model_datafiles, rc) ! Assume that the model mesh is the same as the input data mesh - ! input/output variables + ! input/output variables + type(ESMF_GridComp) :: gcomp type(iosystem_desc_t) , pointer :: pio_subsystem ! pio info integer , intent(in) :: io_type ! pio info integer , intent(in) :: io_format ! pio info @@ -222,6 +225,7 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & type(ESMF_FieldBundle) :: fldbun_noevolve type(ESMF_DistGrid) :: distgrid type(ESMF_Field) :: field_noevolve + type(ESMF_VM) :: vm type(file_desc_t) :: pioid type(io_desc_t) :: pio_iodesc integer :: ns ! ice sheet index @@ -239,8 +243,8 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & real(r8) :: rhoi ! density of ice ~ kg/m^3 real(r8) :: rhoo ! density of sea water ~ kg/m^3 real(r8) :: eus ! eustatic sea level - real(r8) :: Tot_pos_smb ! Sum of positive smb values on each ice sheet for hole-filling - real(r8) :: Tot_neg_smb ! Sum of negative smb values on each ice sheet for hole-filling + real(r8) :: loc_pos_smb(1), Tot_pos_smb(1) ! Sum of positive smb values on each ice sheet for hole-filling + real(r8) :: loc_neg_smb(1), Tot_neg_smb(1) ! Sum of negative smb values on each ice sheet for hole-filling real(r8) :: num_Tot ! Number of active grid cells in total calculation real(r8) :: rat ! Ratio of hole-filling flux to apply real(r8), allocatable :: ice_runoff_out(:) ! Scaled ice runoff output after holes filled @@ -390,42 +394,55 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & allocate(ice_runoff_out(lsize)) ice_runoff_out(:) = 0.d0 - Tot_pos_smb = 0.d0 - Tot_neg_smb = 0.d0 + loc_pos_smb(1) = 0.d0 + Tot_pos_smb(1) = 0.d0 + loc_neg_smb(1) = 0.d0 + Tot_neg_smb(1) = 0.d0 num_Tot = 0.d0 rat = 0.d0 ! For No Evolve to reduce negative ice fluxes from DGLC, we will ! Calculate the total positive and total negative fluxes on each - ! processor for each ice sheet. + ! processor first (local totals). do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) > 0.d0) then - Tot_pos_smb = Tot_pos_smb+Flgl_qice(ns)%ptr(ng) + loc_pos_smb(1) = loc_pos_smb(1)+Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng) end if ! Ignore places that are exactly 0.d0 if(Flgl_qice(ns)%ptr(ng) < 0.d0) then - Tot_neg_smb = Tot_neg_smb+Flgl_qice(ns)%ptr(ng) + loc_neg_smb(1) = loc_neg_smb(1)+Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng) end if end if end do + ! Now do two global sums to get the ice sheet total positive + ! and negative ice fluxes + call ESMF_GridCompGet(gcomp, vm=vm, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=loc_pos_smb, recvdata=Tot_pos_smb, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call ESMF_VMAllreduce(vm, senddata=loc_neg_smb, recvdata=Tot_neg_smb, count=1, & + reduceflag=ESMF_REDUCE_SUM, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! If there's more positive than negative, then set all ! negative to zero and destribute the negative flux amount ! across the positive values, scaled by the size of the ! positive value. This section also applies to any chunks ! where there is no negative smb. In that case the ice ! runoff is exactly equal to the input smb. - if(abs(Tot_pos_smb) >= abs(Tot_neg_smb)) then + if(abs(Tot_pos_smb(1)) >= abs(Tot_neg_smb(1))) then do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) > 0.d0) then - rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb - Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb + rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb(1) + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb(1)/Sg_area(ns)%ptr(ng) else if (Flgl_qice(ns)%ptr(ng) < 0.d0) then Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if else - Fgrg_rofi(ns)%ptr(ng) = 0._r8 + Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if end do else @@ -436,13 +453,13 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) < 0.d0) then - rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb - Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb + rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb(1) + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb(1)/Sg_area(ns)%ptr(ng) else if (Flgl_qice(ns)%ptr(ng) > 0.d0) then Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if else - Fgrg_rofi(ns)%ptr(ng) = 0._r8 + Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if end do diff --git a/dglc/glc_comp_nuopc.F90 b/dglc/glc_comp_nuopc.F90 index 11e66bfd..ab797bbd 100644 --- a/dglc/glc_comp_nuopc.F90 +++ b/dglc/glc_comp_nuopc.F90 @@ -430,7 +430,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) end do ! end loop over ice sheets ! Run dglc - call dglc_comp_run(clock, current_ymd, current_tod, restart_write=.false., valid_inputs=.true., rc=rc) + call dglc_comp_run(gcomp, clock, current_ymd, current_tod, restart_write=.false., valid_inputs=.true., rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_TraceRegionExit('dglc_strdata_init') @@ -498,19 +498,20 @@ subroutine ModelAdvance(gcomp, rc) restart_write = dshr_check_restart_alarm(clock, rc=rc) ! run dglc - call dglc_comp_run(clock, next_ymd, next_tod, restart_write, valid_inputs, rc=rc) + call dglc_comp_run(gcomp, clock, next_ymd, next_tod, restart_write, valid_inputs, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end subroutine ModelAdvance !=============================================================================== - subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, valid_inputs, rc) + subroutine dglc_comp_run(gcomp, clock, target_ymd, target_tod, restart_write, valid_inputs, rc) ! -------------------------- ! advance dglc ! -------------------------- ! input/output variables: + type(ESMF_GridComp) :: gcomp type(ESMF_Clock) , intent(in) :: clock integer , intent(in) :: target_ymd ! model date integer , intent(in) :: target_tod ! model sec into model date @@ -589,7 +590,7 @@ subroutine dglc_comp_run(clock, target_ymd, target_tod, restart_write, valid_inp select case (trim(datamode)) case('noevolve') if (valid_inputs) then - call dglc_datamode_noevolve_advance(sdat(1)%pio_subsystem, sdat(1)%io_type, sdat(1)%io_format, & + call dglc_datamode_noevolve_advance(gcomp, sdat(1)%pio_subsystem, sdat(1)%io_type, sdat(1)%io_format, & logunit, model_meshes, model_internal_gridsize, model_datafiles, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if From 9a8afbc4fc61d334f6c06890ea3ebabd293f9def Mon Sep 17 00:00:00 2001 From: Katetc Date: Fri, 30 Aug 2024 13:44:42 -0600 Subject: [PATCH 5/8] Fixed area weighting on application --- dglc/dglc_datamode_noevolve_mod.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 068e1b81..d544ac8a 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -436,8 +436,8 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) > 0.d0) then - rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb(1) - Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb(1)/Sg_area(ns)%ptr(ng) + rat = Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng)/Tot_pos_smb(1) + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb(1) else if (Flgl_qice(ns)%ptr(ng) < 0.d0) then Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if @@ -453,8 +453,8 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) < 0.d0) then - rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb(1) - Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb(1)/Sg_area(ns)%ptr(ng) + rat = Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng)/Tot_neg_smb(1) + Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb(1) else if (Flgl_qice(ns)%ptr(ng) > 0.d0) then Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if From 5251638b4d8b9d290b97b282b064f59f1fc7f42d Mon Sep 17 00:00:00 2001 From: Katetc Date: Tue, 3 Sep 2024 17:02:51 -0600 Subject: [PATCH 6/8] Removed second area weighting and second if clause to make sure the 0-0 case is handled correctly --- dglc/dglc_datamode_noevolve_mod.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index d544ac8a..3c18f36c 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -436,9 +436,9 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) > 0.d0) then - rat = Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng)/Tot_pos_smb(1) + rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb(1) Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb(1) - else if (Flgl_qice(ns)%ptr(ng) < 0.d0) then + else Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if else @@ -453,9 +453,9 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form do ng = 1,lsize if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then if(Flgl_qice(ns)%ptr(ng) < 0.d0) then - rat = Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng)/Tot_neg_smb(1) + rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb(1) Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb(1) - else if (Flgl_qice(ns)%ptr(ng) > 0.d0) then + else Fgrg_rofi(ns)%ptr(ng) = 0.d0 end if else From 34838d24dcd29ff00b74bdc0edf60db132c42254 Mon Sep 17 00:00:00 2001 From: Katetc Date: Thu, 5 Sep 2024 10:43:24 -0600 Subject: [PATCH 7/8] One more review comment addressed --- dglc/dglc_datamode_noevolve_mod.F90 | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index fef1201b..3b64c884 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -398,9 +398,8 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form ! Get mesh info call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - call ESMF_DistGridGet(distGrid, localDe=0, elementCount=lsize, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return - + lsize = size(Fgrg_rofi(ns)%ptr) + allocate(ice_runoff_out(lsize)) ice_runoff_out(:) = 0.d0 loc_pos_smb(1) = 0.d0 From 9b22a506d661f50f71d03e1d73170c2b2d4b3663 Mon Sep 17 00:00:00 2001 From: Katetc Date: Thu, 5 Sep 2024 15:17:48 -0600 Subject: [PATCH 8/8] Final code cleanup and review comments --- dglc/dglc_datamode_noevolve_mod.F90 | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 3b64c884..6a7caa58 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -261,9 +261,7 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form real(r8) :: usrf ! upper surface elevation (m) on ice grid real(r8) :: loc_pos_smb(1), Tot_pos_smb(1) ! Sum of positive smb values on each ice sheet for hole-filling real(r8) :: loc_neg_smb(1), Tot_neg_smb(1) ! Sum of negative smb values on each ice sheet for hole-filling - real(r8) :: num_Tot ! Number of active grid cells in total calculation real(r8) :: rat ! Ratio of hole-filling flux to apply - real(r8), allocatable :: ice_runoff_out(:) ! Scaled ice runoff output after holes filled character(len=*), parameter :: subname='(dglc_datamode_noevolve_advance): ' !------------------------------------------------------------------------------- @@ -395,25 +393,22 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form ! Compute Fgrg_rofi do ns = 1,num_icesheets - ! Get mesh info - call ESMF_MeshGet(model_meshes(ns), elementdistGrid=distGrid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! Get number of grid cells per ice sheet lsize = size(Fgrg_rofi(ns)%ptr) - - allocate(ice_runoff_out(lsize)) - ice_runoff_out(:) = 0.d0 + + ! reset output variables to zero + Fgrg_rofi(ns)%ptr(:) = 0.d0 loc_pos_smb(1) = 0.d0 Tot_pos_smb(1) = 0.d0 loc_neg_smb(1) = 0.d0 Tot_neg_smb(1) = 0.d0 - num_Tot = 0.d0 rat = 0.d0 ! For No Evolve to reduce negative ice fluxes from DGLC, we will ! Calculate the total positive and total negative fluxes on each ! processor first (local totals). do ng = 1,lsize - if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng) > 0.d0) then if(Flgl_qice(ns)%ptr(ng) > 0.d0) then loc_pos_smb(1) = loc_pos_smb(1)+Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng) end if @@ -433,7 +428,7 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form call ESMF_VMAllreduce(vm, senddata=loc_neg_smb, recvdata=Tot_neg_smb, count=1, & reduceflag=ESMF_REDUCE_SUM, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return - + ! If there's more positive than negative, then set all ! negative to zero and destribute the negative flux amount ! across the positive values, scaled by the size of the @@ -442,7 +437,7 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form ! runoff is exactly equal to the input smb. if(abs(Tot_pos_smb(1)) >= abs(Tot_neg_smb(1))) then do ng = 1,lsize - if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng) > 0.d0) then if(Flgl_qice(ns)%ptr(ng) > 0.d0) then rat = Flgl_qice(ns)%ptr(ng)/Tot_pos_smb(1) Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_neg_smb(1) @@ -459,7 +454,7 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form ! reduce their negativity a bit. This shouldn't happen often. ! This section of code also applies if Tot_pos_smb is zero. do ng = 1,lsize - if (Sg_icemask_coupled_fluxes(ns)%ptr(ng).gt.0.d0) then + if (Sg_icemask_coupled_fluxes(ns)%ptr(ng) > 0.d0) then if(Flgl_qice(ns)%ptr(ng) < 0.d0) then rat = Flgl_qice(ns)%ptr(ng)/Tot_neg_smb(1) Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) + rat*Tot_pos_smb(1) @@ -473,9 +468,8 @@ subroutine dglc_datamode_noevolve_advance(gcomp, pio_subsystem, io_type, io_form end if ! More neg or pos smb - deallocate(ice_runoff_out) - end do - end if + end do ! Each ice sheet + end if ! if initialized ! Set initialized flag initialized_noevolve = .true.