diff --git a/dglc/dglc_datamode_noevolve_mod.F90 b/dglc/dglc_datamode_noevolve_mod.F90 index 3393ea21..6a7caa58 100644 --- a/dglc/dglc_datamode_noevolve_mod.F90 +++ b/dglc/dglc_datamode_noevolve_mod.F90 @@ -4,7 +4,9 @@ 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_VMGetCurrent, ESMF_VMBroadCast, ESMF_VM + use ESMF , only : ESMF_GridComp, ESMF_GridCompGet + use ESMF , only : ESMF_VM, ESMF_VMAllreduce, ESMF_REDUCE_SUM + use ESMF , only : ESMF_VMGetCurrent, ESMF_VMBroadCast 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 @@ -217,15 +219,17 @@ 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, & - model_meshes, model_internal_gridsize, model_datafiles, rc) + 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 + 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 @@ -235,6 +239,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 @@ -254,6 +259,10 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & real(r8) :: eus ! eustatic sea level real(r8) :: lsrf ! lower surface elevation (m) on ice grid 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) :: rat ! Ratio of hole-filling flux to apply + character(len=*), parameter :: subname='(dglc_datamode_noevolve_advance): ' !------------------------------------------------------------------------------- @@ -380,25 +389,91 @@ subroutine dglc_datamode_noevolve_advance(pio_subsystem, io_type, io_format, & end if if (initialized_noevolve) then + ! Compute Fgrg_rofi do ns = 1,num_icesheets - 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 + + ! Get number of grid cells per ice sheet + lsize = size(Fgrg_rofi(ns)%ptr) + + ! 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 + 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) > 0.d0) then - Fgrg_rofi(ns)%ptr(ng) = Flgl_qice(ns)%ptr(ng) - else - Fgrg_rofi(ns)%ptr(ng) = 0._r8 + 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 + ! Ignore places that are exactly 0.d0 + if(Flgl_qice(ns)%ptr(ng) < 0.d0) then + loc_neg_smb(1) = loc_neg_smb(1)+Flgl_qice(ns)%ptr(ng)*Sg_area(ns)%ptr(ng) + end if end if end do - end do - end if + ! 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(1)) >= abs(Tot_neg_smb(1))) then + do ng = 1,lsize + 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) + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + end do + else + ! 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) > 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) + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + else + Fgrg_rofi(ns)%ptr(ng) = 0.d0 + end if + end do + + end if ! More neg or pos smb + + end do ! Each ice sheet + end if ! if initialized ! Set initialized flag initialized_noevolve = .true. - + end subroutine dglc_datamode_noevolve_advance !=============================================================================== diff --git a/dglc/glc_comp_nuopc.F90 b/dglc/glc_comp_nuopc.F90 index b74e56f7..68623456 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') @@ -503,7 +503,7 @@ subroutine ModelAdvance(gcomp, rc) end if ! 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 call ESMF_TraceRegionExit(subname) @@ -511,13 +511,14 @@ subroutine ModelAdvance(gcomp, rc) 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 @@ -598,8 +599,8 @@ 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, & - model_meshes, model_internal_gridsize, model_datafiles, rc) + 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 end select