From e7ce959e8b59a0c67e793706beab62e8765f30db Mon Sep 17 00:00:00 2001 From: Stephen Price Date: Tue, 2 Jul 2024 17:58:47 -0500 Subject: [PATCH 01/13] Improve and debug coupler budgets for GLC Update various bits of code needed for completing and testing of glc coupler budgets. Still very much in the testing and debugging phase. --- .../mpas-albany-landice/driver/glc_comp_mct.F | 2 +- driver-mct/main/cime_comp_mod.F90 | 9 +- driver-mct/main/seq_diag_mct.F90 | 124 +++++++++++++----- 3 files changed, 97 insertions(+), 38 deletions(-) diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F index ac33d17b1f31..0146fc74791c 100644 --- a/components/mpas-albany-landice/driver/glc_comp_mct.F +++ b/components/mpas-albany-landice/driver/glc_comp_mct.F @@ -1482,7 +1482,7 @@ subroutine glc_export_mct(g2x_g, errorCode) !call route_ice_runoff(0.0_RKIND, & !Recuperate runoff routing switch code (originally in glc_route_ice_runoff module in earlier code), and attach to ice calving flux once present... ! rofi_to_ocn=Fogg_rofi, & ! rofi_to_ice=Figg_rofi) - g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0 !...and remove these placeholders + g2x_g % rAttr(index_g2x_Fogg_rofi,n)=9999.0 !...and remove these placeholders g2x_g % rAttr(index_g2x_Figg_rofi,n)=0.0 !...and remove these placeholders g2x_g % rAttr(index_g2x_Fogg_rofl,n) = 0.0 !Attach to subglacial liquid flux once present diff --git a/driver-mct/main/cime_comp_mod.F90 b/driver-mct/main/cime_comp_mod.F90 index 2131d0c86844..77749415af8c 100644 --- a/driver-mct/main/cime_comp_mod.F90 +++ b/driver-mct/main/cime_comp_mod.F90 @@ -145,7 +145,8 @@ module cime_comp_mod ! diagnostic routines use seq_diag_mct, only : seq_diag_zero_mct , seq_diag_avect_mct, seq_diag_lnd_mct use seq_diag_mct, only : seq_diag_rof_mct , seq_diag_ocn_mct , seq_diag_atm_mct - use seq_diag_mct, only : seq_diag_ice_mct , seq_diag_accum_mct, seq_diag_print_mct + use seq_diag_mct, only : seq_diag_ice_mct , seq_diag_glc_mct + use seq_diag_mct, only : seq_diag_accum_mct, seq_diag_print_mct use seq_diagBGC_mct, only : seq_diagBGC_zero_mct , seq_diagBGC_avect_mct, seq_diagBGC_lnd_mct use seq_diagBGC_mct, only : seq_diagBGC_rof_mct , seq_diagBGC_ocn_mct , seq_diagBGC_atm_mct use seq_diagBGC_mct, only : seq_diagBGC_ice_mct , seq_diagBGC_accum_mct @@ -4743,6 +4744,9 @@ subroutine cime_run_calc_budgets1(in_cplrun) if (ice_present) then call seq_diag_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_x2i=.true.) endif + if (glc_present) then + call seq_diag_glc_mct(glc(ens1), fractions_ix(ens1), infodata, do_x2g=.true., do_g2x=.true.) + endif if (do_bgc_budgets) then if (rof_present) then call seq_diagBGC_rof_mct(rof(ens1), fractions_rx(ens1), infodata) @@ -4789,6 +4793,9 @@ subroutine cime_run_calc_budgets2(in_cplrun) if (ice_present) then call seq_diagBGC_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_i2x=.true., do_x2i=.true.) endif + if (glc_present) then + call seq_diag_glc_mct(glc(ens1), fractions_ix(ens1), infodata, do_x2g=.true., do_g2x=.true.) + endif if (lnd_present) then call seq_diagBGC_lnd_mct(lnd(ens1), fractions_lx(ens1), infodata, do_l2x=.true., do_x2l=.true.) endif diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 8534008f9be7..502373fc7a90 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -140,42 +140,45 @@ module seq_diag_mct integer(in),parameter :: f_hsen =10 ! heat : sensible integer(in),parameter :: f_hpolar =11 ! heat : AIS imbalance integer(in),parameter :: f_hh2ot =12 ! heat : water temperature - integer(in),parameter :: f_wfrz =13 ! water: freezing - integer(in),parameter :: f_wmelt =14 ! water: melting - integer(in),parameter :: f_wrain =15 ! water: precip, liquid - integer(in),parameter :: f_wsnow =16 ! water: precip, frozen - integer(in),parameter :: f_wpolar =17 ! water: AIS imbalance - integer(in),parameter :: f_wevap =18 ! water: evaporation - integer(in),parameter :: f_wroff =19 ! water: runoff/flood - integer(in),parameter :: f_wioff =20 ! water: frozen runoff - integer(in),parameter :: f_wirrig =21 ! water: irrigation - integer(in),parameter :: f_wfrz_16O =22 ! water: freezing - integer(in),parameter :: f_wmelt_16O =23 ! water: melting - integer(in),parameter :: f_wrain_16O =24 ! water: precip, liquid - integer(in),parameter :: f_wsnow_16O =25 ! water: precip, frozen - integer(in),parameter :: f_wevap_16O =26 ! water: evaporation - integer(in),parameter :: f_wroff_16O =27 ! water: runoff/flood - integer(in),parameter :: f_wioff_16O =28 ! water: frozen runoff - integer(in),parameter :: f_wfrz_18O =29 ! water: freezing - integer(in),parameter :: f_wmelt_18O =30 ! water: melting - integer(in),parameter :: f_wrain_18O =31 ! water: precip, liquid - integer(in),parameter :: f_wsnow_18O =32 ! water: precip, frozen - integer(in),parameter :: f_wevap_18O =33 ! water: evaporation - integer(in),parameter :: f_wroff_18O =34 ! water: runoff/flood - integer(in),parameter :: f_wioff_18O =35 ! water: frozen runoff - integer(in),parameter :: f_wfrz_HDO =36 ! water: freezing - integer(in),parameter :: f_wmelt_HDO =37 ! water: melting - integer(in),parameter :: f_wrain_HDO =38 ! water: precip, liquid - integer(in),parameter :: f_wsnow_HDO =39 ! water: precip, frozen - integer(in),parameter :: f_wevap_HDO =40 ! water: evaporation - integer(in),parameter :: f_wroff_HDO =41 ! water: runoff/flood - integer(in),parameter :: f_wioff_HDO =42 ! water: frozen runoff + integer(in),parameter :: f_hgsmb =13 ! heat : GIS SMB !SFP added + integer(in),parameter :: f_wfrz =14 ! water: freezing + integer(in),parameter :: f_wmelt =15 ! water: melting + integer(in),parameter :: f_wrain =16 ! water: precip, liquid + integer(in),parameter :: f_wsnow =17 ! water: precip, frozen + integer(in),parameter :: f_wpolar =18 ! water: AIS imbalance + integer(in),parameter :: f_wgsmb =19 ! water: GIS SMB !SFP added + integer(in),parameter :: f_wevap =20 ! water: evaporation + integer(in),parameter :: f_wroff =21 ! water: runoff/flood + integer(in),parameter :: f_wioff =22 ! water: frozen runoff + integer(in),parameter :: f_wirrig =23 ! water: irrigation + integer(in),parameter :: f_wfrz_16O =24 ! water: freezing + integer(in),parameter :: f_wmelt_16O =25 ! water: melting + integer(in),parameter :: f_wrain_16O =26 ! water: precip, liquid + integer(in),parameter :: f_wsnow_16O =27 ! water: precip, frozen + integer(in),parameter :: f_wevap_16O =28 ! water: evaporation + integer(in),parameter :: f_wroff_16O =29 ! water: runoff/flood + integer(in),parameter :: f_wioff_16O =30 ! water: frozen runoff + integer(in),parameter :: f_wfrz_18O =31 ! water: freezing + integer(in),parameter :: f_wmelt_18O =32 ! water: melting + integer(in),parameter :: f_wrain_18O =33 ! water: precip, liquid + integer(in),parameter :: f_wsnow_18O =34 ! water: precip, frozen + integer(in),parameter :: f_wevap_18O =35 ! water: evaporation + integer(in),parameter :: f_wroff_18O =36 ! water: runoff/flood + integer(in),parameter :: f_wioff_18O =37 ! water: frozen runoff + integer(in),parameter :: f_wfrz_HDO =38 ! water: freezing + integer(in),parameter :: f_wmelt_HDO =39 ! water: melting + integer(in),parameter :: f_wrain_HDO =40 ! water: precip, liquid + integer(in),parameter :: f_wsnow_HDO =41 ! water: precip, frozen + integer(in),parameter :: f_wevap_HDO =42 ! water: evaporation + integer(in),parameter :: f_wroff_HDO =43 ! water: runoff/flood + integer(in),parameter :: f_wioff_HDO =44 ! water: frozen runoff integer(in),parameter :: f_size = f_wioff_HDO ! Total array size of all elements integer(in),parameter :: f_a = f_area ! 1st index for area integer(in),parameter :: f_a_end = f_area ! last index for area integer(in),parameter :: f_h = f_hfrz ! 1st index for heat - integer(in),parameter :: f_h_end = f_hh2ot ! Last index for heat + !integer(in),parameter :: f_h_end = f_hh2ot ! Last index for heat + integer(in),parameter :: f_h_end = f_hgsmb ! Last index for heat integer(in),parameter :: f_w = f_wfrz ! 1st index for water integer(in),parameter :: f_w_end = f_wirrig ! Last index for water integer(in),parameter :: f_16O = f_wfrz_16O ! 1st index for 16O water isotope @@ -189,8 +192,10 @@ module seq_diag_mct (/' area',' hfreeze',' hmelt',' hnetsw',' hlwdn', & ' hlwup',' hlatvap',' hlatfus',' hiroff',' hsen', & - ' hpolar',' hh2otemp',' wfreeze',' wmelt',' wrain', & - ' wsnow',' wpolar',' wevap',' wrunoff',' wfrzrof', & +! ' hpolar',' hh2otemp',' wfreeze',' wmelt',' wrain', & + ' hpolar',' hh2otemp',' hgsmb',' wfreeze',' wmelt',' wrain', & +! ' wsnow',' wpolar',' wevap',' wrunoff',' wfrzrof', & + ' wsnow',' wpolar',' wgsmb',' wevap',' wrunoff',' wfrzrof', & ' wirrig', & ' wfreeze_16O',' wmelt_16O',' wrain_16O',' wsnow_16O', & ' wevap_16O',' wrunoff_16O',' wfrzrof_16O', & @@ -262,6 +267,7 @@ module seq_diag_mct integer :: index_l2x_Flrl_irrig integer :: index_l2x_Flrl_wslake + integer :: index_l2x_Flgl_qice(0:10) !SFP added integer :: index_x2l_Faxa_lwdn integer :: index_x2l_Faxa_rainc @@ -338,6 +344,8 @@ module seq_diag_mct integer :: index_g2x_Fogg_rofi integer :: index_g2x_Figg_rofi + integer :: index_x2g_Flgl_qice !SFP added + integer :: index_x2o_Foxx_rofl_16O integer :: index_x2o_Foxx_rofi_16O integer :: index_x2o_Foxx_rofl_18O @@ -874,6 +882,10 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) logical,save :: first_time = .true. logical,save :: flds_wiso_lnd = .false. + character(len=64) :: name !SFP: added this and next 2 + character(len= 2) :: cnum + integer(in) :: num + !----- formats ----- character(*),parameter :: subName = '(seq_diag_lnd_mct) ' @@ -909,6 +921,13 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_irrig = mct_aVect_indexRA(l2x_l,'Flrl_irrig', perrWith='quiet') index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') + do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) + write(cnum,'(i2.2)') num + name = 'Flgl_qice' // cnum + index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) !SFP added + !index_l2x_Flgl_qice = mct_aVect_indexRA(l2x_l,'Flgl_qice') + end do + index_l2x_Fall_evap_16O = mct_aVect_indexRA(l2x_l,'Fall_evap_16O',perrWith='quiet') if ( index_l2x_Fall_evap_16O /= 0 ) flds_wiso_lnd = .true. if ( flds_wiso_lnd )then @@ -944,6 +963,10 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) end if nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) + do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) + nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flgl_qice(num),n) !SFP added + end do + if ( flds_wiso_lnd )then nf = f_wevap_16O; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + & @@ -976,7 +999,11 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi_HDO,n) end if end do + budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice + + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP add + end if if (present(do_x2l)) then @@ -1256,11 +1283,13 @@ end subroutine seq_diag_rof_mct ! ! !INTERFACE: ------------------------------------------------------------------ - subroutine seq_diag_glc_mct( glc, frac_g, infodata) + subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) type(component_type) , intent(in) :: glc ! component type for instance1 type(mct_aVect) , intent(in) :: frac_g ! frac bundle type(seq_infodata_type) , intent(in) :: infodata + logical , intent(in), optional :: do_x2g + logical , intent(in), optional :: do_g2x !EOP @@ -1289,14 +1318,25 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata) g2x_g => component_get_c2x_cx(glc) x2g_g => component_get_x2c_cx(glc) - if (first_time) then +! eventually use the following if constructs to wrap relevant sections below? +! if (present(do_g2x)) then +! end if +! if (present(do_x2g)) then +! end if + +! if (first_time) then + index_g2x_Fogg_rofl = mct_aVect_indexRA(g2x_g,'Fogg_rofl') index_g2x_Fogg_rofi = mct_aVect_indexRA(g2x_g,'Fogg_rofi') index_g2x_Figg_rofi = mct_aVect_indexRA(g2x_g,'Figg_rofi') - end if + + index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') !SFP: might be cleaner to do "x2g" in its own section? + +! end if ip = p_inst ic = c_glc_gs + !ic = c_glc_gr !SFP: should this actually be used here? other sections of this code use"r" for c2x and "s" for x2c kArea = mct_aVect_indexRA(dom_g%data,afldname) lSize = mct_avect_lSize(g2x_g) do n=1,lSize @@ -1307,6 +1347,18 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata) end do budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice + ip = p_inst + !SFP: unclear if next should be 'send' or 'receive' index but for now we think send (because of x2g), + ! but budget results don't seem to be sensitive to this choice (same numbers appear using either). + ic = c_glc_gs ! cpl send + !ic = c_glc_gr ! cpl receive + lSize = mct_avect_lSize(x2g_g) + do n=1,lSize + ca_g = dom_g%data%rAttr(kArea,n) + nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) + end do + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice + first_time = .false. end subroutine seq_diag_glc_mct From 855468acd50031d50340771ab6e01b4cbd53dc9f Mon Sep 17 00:00:00 2001 From: Stephen Price Date: Tue, 27 Aug 2024 13:01:46 -0500 Subject: [PATCH 02/13] Continuation of glc coupler budget development Code additions and debugging for addition of Greenland surface mass balance terms to glc budget code, including temporary lines for debugging and validation. --- .../mpas-albany-landice/driver/glc_comp_mct.F | 3 +- driver-mct/main/cime_comp_mod.F90 | 8 +- driver-mct/main/seq_diag_mct.F90 | 156 +++++++++++++----- 3 files changed, 122 insertions(+), 45 deletions(-) diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F index 0146fc74791c..7bf72556603c 100644 --- a/components/mpas-albany-landice/driver/glc_comp_mct.F +++ b/components/mpas-albany-landice/driver/glc_comp_mct.F @@ -1482,7 +1482,8 @@ subroutine glc_export_mct(g2x_g, errorCode) !call route_ice_runoff(0.0_RKIND, & !Recuperate runoff routing switch code (originally in glc_route_ice_runoff module in earlier code), and attach to ice calving flux once present... ! rofi_to_ocn=Fogg_rofi, & ! rofi_to_ice=Figg_rofi) - g2x_g % rAttr(index_g2x_Fogg_rofi,n)=9999.0 !...and remove these placeholders + !g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0!...and remove these placeholders + g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0001d0 !SFP: dummy value to see if passes through coupler g2x_g % rAttr(index_g2x_Figg_rofi,n)=0.0 !...and remove these placeholders g2x_g % rAttr(index_g2x_Fogg_rofl,n) = 0.0 !Attach to subglacial liquid flux once present diff --git a/driver-mct/main/cime_comp_mod.F90 b/driver-mct/main/cime_comp_mod.F90 index 77749415af8c..f9678b3e7c9d 100644 --- a/driver-mct/main/cime_comp_mod.F90 +++ b/driver-mct/main/cime_comp_mod.F90 @@ -4745,7 +4745,7 @@ subroutine cime_run_calc_budgets1(in_cplrun) call seq_diag_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_x2i=.true.) endif if (glc_present) then - call seq_diag_glc_mct(glc(ens1), fractions_ix(ens1), infodata, do_x2g=.true., do_g2x=.true.) + !call seq_diag_glc_mct(glc(ens1), fractions_gx(ens1), infodata, do_x2g=.true., do_g2x=.true.) !SFP: comment out for now while debugging endif if (do_bgc_budgets) then if (rof_present) then @@ -4786,6 +4786,9 @@ subroutine cime_run_calc_budgets2(in_cplrun) if (ice_present) then call seq_diag_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_i2x=.true.) endif + if (glc_present) then + call seq_diag_glc_mct(glc(ens1), fractions_gx(ens1), infodata, do_x2g=.true., do_g2x=.true.) + endif if (do_bgc_budgets) then if (atm_present) then call seq_diagBGC_atm_mct(atm(ens1), fractions_ax(ens1), infodata, do_a2x=.true., do_x2a=.true.) @@ -4793,9 +4796,6 @@ subroutine cime_run_calc_budgets2(in_cplrun) if (ice_present) then call seq_diagBGC_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_i2x=.true., do_x2i=.true.) endif - if (glc_present) then - call seq_diag_glc_mct(glc(ens1), fractions_ix(ens1), infodata, do_x2g=.true., do_g2x=.true.) - endif if (lnd_present) then call seq_diagBGC_lnd_mct(lnd(ens1), fractions_lx(ens1), infodata, do_l2x=.true., do_x2l=.true.) endif diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 502373fc7a90..734b375070dc 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -963,9 +963,10 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) end if nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) - do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) - nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flgl_qice(num),n) !SFP added - end do +! do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) +! !SFP: this somehow needs to allow for each of the 11 vectors associate w/ each of the 11 elev classes +! nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flgl_qice(num),n) !SFP added +! end do if ( flds_wiso_lnd )then nf = f_wevap_16O; @@ -1000,9 +1001,9 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) end if end do - budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice +! budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice - budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP add +! budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP added end if @@ -1286,7 +1287,7 @@ end subroutine seq_diag_rof_mct subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) type(component_type) , intent(in) :: glc ! component type for instance1 - type(mct_aVect) , intent(in) :: frac_g ! frac bundle + type(mct_aVect) , intent(in) :: frac_g ! frac bundle !SFP: does not look like fractions are needed / used here? type(seq_infodata_type) , intent(in) :: infodata logical , intent(in), optional :: do_x2g logical , intent(in), optional :: do_g2x @@ -1303,6 +1304,9 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) real(r8) :: ca_g ! area of a grid cell logical,save :: first_time = .true. + integer,save :: counter,smb_counter,calving_counter !SFP: for debugging + integer,save :: smb_vector_length,calving_vector_length + !----- formats ----- character(*),parameter :: subName = '(seq_diag_glc_mct) ' @@ -1318,46 +1322,118 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) g2x_g => component_get_c2x_cx(glc) x2g_g => component_get_x2c_cx(glc) -! eventually use the following if constructs to wrap relevant sections below? -! if (present(do_g2x)) then -! end if -! if (present(do_x2g)) then -! end if + if( present(do_g2x))then !SPF: glc to coupler -! if (first_time) then + if (first_time) then - index_g2x_Fogg_rofl = mct_aVect_indexRA(g2x_g,'Fogg_rofl') - index_g2x_Fogg_rofi = mct_aVect_indexRA(g2x_g,'Fogg_rofi') - index_g2x_Figg_rofi = mct_aVect_indexRA(g2x_g,'Figg_rofi') + calving_counter=0 + calving_vector_length = 0 - index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') !SFP: might be cleaner to do "x2g" in its own section? + index_g2x_Fogg_rofl = mct_aVect_indexRA(g2x_g,'Fogg_rofl') + index_g2x_Fogg_rofi = mct_aVect_indexRA(g2x_g,'Fogg_rofi') + index_g2x_Figg_rofi = mct_aVect_indexRA(g2x_g,'Figg_rofi') -! end if + !SFP:debug + write(logunit,*) ' ' + write(logunit,*) ' index_g2x_Fogg_rofl = ', index_g2x_Fogg_rofl + write(logunit,*) ' index_g2x_Fogg_rofi = ', index_g2x_Fogg_rofi + write(logunit,*) ' index_g2x_Figg_rofi = ', index_g2x_Figg_rofi + write(logunit,*) ' ' - ip = p_inst - ic = c_glc_gs - !ic = c_glc_gr !SFP: should this actually be used here? other sections of this code use"r" for c2x and "s" for x2c - kArea = mct_aVect_indexRA(dom_g%data,afldname) - lSize = mct_avect_lSize(g2x_g) - do n=1,lSize - ca_g = dom_g%data%rAttr(kArea,n) - nf = f_wroff; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_g*g2x_g%rAttr(index_g2x_Fogg_rofl,n) - nf = f_wioff; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_g*g2x_g%rAttr(index_g2x_Fogg_rofi,n) & - - ca_g*g2x_g%rAttr(index_g2x_Figg_rofi,n) - end do - budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice + end if - ip = p_inst - !SFP: unclear if next should be 'send' or 'receive' index but for now we think send (because of x2g), - ! but budget results don't seem to be sensitive to this choice (same numbers appear using either). - ic = c_glc_gs ! cpl send - !ic = c_glc_gr ! cpl receive - lSize = mct_avect_lSize(x2g_g) - do n=1,lSize - ca_g = dom_g%data%rAttr(kArea,n) - nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) - end do - budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice + ip = p_inst + !ip = p_day + !ic = c_glc_gs + ic = c_glc_gr !SFP: use recieve here since this is coming from glc to coupler? + kArea = mct_aVect_indexRA(dom_g%data,afldname) + lSize = mct_avect_lSize(g2x_g) + + !SFP:debug + if(calving_counter==0)then !one day at 30 min land/atmos time steps + write(logunit,*) ' ' + write(logunit,*) ' calving vector length (should be 7425) = ', lSize + write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) + write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) + write(logunit,*) ' intial value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) + write(logunit,*) ' calving flux to ocean (Fogg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Fogg_rofi,1) + write(logunit,*) ' calving flux to ice (Figg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Figg_rofi,1) + write(logunit,*) ' calving flux X area to ocean (Fogg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Fogg_rofi,1) + write(logunit,*) ' calving flux X area to ice (Figg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Figg_rofi,1) + end if + + do n=1,lSize + ca_g = dom_g%data%rAttr(kArea,n) + nf = f_wroff; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_g*g2x_g%rAttr(index_g2x_Fogg_rofl,n) + nf = f_wioff; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_g*g2x_g%rAttr(index_g2x_Fogg_rofi,n) & + - ca_g*g2x_g%rAttr(index_g2x_Figg_rofi,n) + end do + + budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice + + calving_vector_length = calving_vector_length +lSize + calving_counter = calving_counter + 1 + + !SFP:debug + if(calving_counter==48)then !one day at 30 min land/atmos time steps + write(logunit,*) ' calving counter = ', calving_counter + write(logunit,*) ' final value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) + write(logunit,*) ' ' + end if + + endif !SFP: end 'do_g2x' + + if( present(do_x2g))then !SFP: coupler to glc + + if (first_time) then + + smb_counter=0 + smb_vector_length = 0 + + index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') + + !SFP:debug + write(logunit,*) ' ' + write(logunit,*) ' index_x2g_Flgl_qice = ', index_x2g_Flgl_qice + write(logunit,*) ' ' + + end if + + ip = p_inst + !ip = p_day + ic = c_glc_gs ! SFP: use send here since going from coupler to glc? + !ic = c_glc_gr + kArea = mct_aVect_indexRA(dom_g%data,afldname) + lSize = mct_avect_lSize(x2g_g) + + !SFP:debug + if(smb_counter==0)then !one day at 30 min land/atmos time steps + write(logunit,*) ' ' + write(logunit,*) ' smb vector length (should be 7425) = ', lSize + write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) + write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) + write(logunit,*) ' initial value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) + end if + + do n=1,lSize + ca_g = dom_g%data%rAttr(kArea,n) + nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) + end do + + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice + + smb_vector_length = smb_vector_length +lSize + smb_counter = smb_counter + 1 + + !SFP:debug + if(smb_counter==48)then !one day at 30 min land/atmos time steps + write(logunit,*) ' ' + write(logunit,*) ' smb_counter = ', smb_counter + write(logunit,*) ' final value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) + write(logunit,*) ' ' + end if + + end if !SPF: end do coupler to glc first_time = .false. From 564a80beb817a1089156ee9653df5184295c86d0 Mon Sep 17 00:00:00 2001 From: Stephen Price Date: Thu, 29 Aug 2024 18:30:50 -0500 Subject: [PATCH 03/13] Further work on adding support for and testing of glc budgets for GIS component Given known number of time steps smb flux is accumulated over, x2g_ smb flux term in budget table now agrees with value calculated from mali and cpl hist files (still need to add support for automatically determining number of averaging steps). Similar support on l2x_ side is close (summing product of accumulation fluxes in different MECs w/ area fractions and cell areas) but off in budget tables from cpl calculated values by ~5%. --- driver-mct/main/seq_diag_mct.F90 | 64 +++++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 18 deletions(-) diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 734b375070dc..af26aeacfbc2 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -48,6 +48,9 @@ module seq_diag_mct use shr_reprosum_mod, only : shr_reprosum_calc use seq_diagBGC_mct, only : seq_diagBGC_preprint_mct, seq_diagBGC_print_mct + use prep_glc_mod, only : prep_glc_get_x2gacc_gx_cnt !SFP: added this and next - unclear which is needed +! use prep_glc_mod, only : prep_glc_get_l2gacc_lx_cnt + implicit none save private @@ -219,6 +222,10 @@ module seq_diag_mct logical :: flds_wiso ! If water isotope fields are active + !--- temporary pointers --- + integer , pointer :: x2gacc_gx_cnt ! SFP added +! integer , pointer :: l2gacc_lx_cnt ! SFP: unclear if this or the above is needed / more relevant + ! !PUBLIC DATA MEMBERS !--- time-averaged (annual?) global budge diagnostics --- @@ -268,6 +275,7 @@ module seq_diag_mct integer :: index_l2x_Flrl_wslake integer :: index_l2x_Flgl_qice(0:10) !SFP added + integer :: index_x2l_Sg_ice_covered(0:10) !SFP added integer :: index_x2l_Faxa_lwdn integer :: index_x2l_Faxa_rainc @@ -882,6 +890,8 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) logical,save :: first_time = .true. logical,save :: flds_wiso_lnd = .false. + real(r8) :: l2x_Flgl_qice_col_sum !SFP: sum of fluxes over no. MECs (cols) + character(len=64) :: name !SFP: added this and next 2 character(len= 2) :: cnum integer(in) :: num @@ -921,11 +931,12 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_irrig = mct_aVect_indexRA(l2x_l,'Flrl_irrig', perrWith='quiet') index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') - do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) + do num=0,10 !SFP: change later to 0,glc_nec_max (no. of elev classes) write(cnum,'(i2.2)') num name = 'Flgl_qice' // cnum index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) !SFP added - !index_l2x_Flgl_qice = mct_aVect_indexRA(l2x_l,'Flgl_qice') + name = 'Sg_ice_covered' // cnum + index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) !SFP added end do index_l2x_Fall_evap_16O = mct_aVect_indexRA(l2x_l,'Fall_evap_16O',perrWith='quiet') @@ -963,10 +974,12 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) end if nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) -! do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) -! !SFP: this somehow needs to allow for each of the 11 vectors associate w/ each of the 11 elev classes -! nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flgl_qice(num),n) !SFP added -! end do + l2x_Flgl_qice_col_sum = 0.0d0 + do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) + !SFP: this somehow needs to allow for each of the 11 vectors associate w/ each of the 11 elev classes + l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) !SFP added + end do + nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_Flgl_qice_col_sum !SFP added if ( flds_wiso_lnd )then nf = f_wevap_16O; @@ -1003,7 +1016,7 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) ! budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice -! budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP added + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP added end if @@ -1296,7 +1309,8 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) !----- local ----- type(mct_aVect), pointer :: g2x_g - type(mct_aVect), pointer :: x2g_g +! type(mct_aVect), pointer :: x2g_g + type(mct_aVect), pointer :: x2gacc_g !SFP: replace above w/ vector for accumulated fluxes type(mct_ggrid), pointer :: dom_g integer(in) :: n,ic,nf,ip ! generic index integer(in) :: kArea ! index of area field in aVect @@ -1320,7 +1334,11 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) dom_g => component_get_dom_cx(glc) g2x_g => component_get_c2x_cx(glc) - x2g_g => component_get_x2c_cx(glc) +! x2g_g => component_get_x2c_cx(glc) + x2gacc_g => component_get_x2c_cx(glc) !SFP: use accum fluxes vector + + x2gacc_gx_cnt => prep_glc_get_x2gacc_gx_cnt() !SFP: counter for how many times SMB flux accumulation has occured +! l2gacc_lx_cnt => prep_glc_get l2gacc_lx_cnt() if( present(do_g2x))then !SPF: glc to coupler @@ -1342,8 +1360,8 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) end if - ip = p_inst - !ip = p_day + !ip = p_inst + ip = p_day !ic = c_glc_gs ic = c_glc_gr !SFP: use recieve here since this is coming from glc to coupler? kArea = mct_aVect_indexRA(dom_g%data,afldname) @@ -1352,7 +1370,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) !SFP:debug if(calving_counter==0)then !one day at 30 min land/atmos time steps write(logunit,*) ' ' - write(logunit,*) ' calving vector length (should be 7425) = ', lSize + write(logunit,*) ' calving vector length (7425 in coupler) = ', lSize write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) write(logunit,*) ' intial value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) @@ -1390,7 +1408,8 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) smb_counter=0 smb_vector_length = 0 - index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') + !index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') + index_x2g_Flgl_qice = mct_aVect_indexRA(x2gacc_g,'Flgl_qice') !SFP: use accum flux vector !SFP:debug write(logunit,*) ' ' @@ -1399,17 +1418,18 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) end if - ip = p_inst - !ip = p_day + !ip = p_inst + ip = p_day ic = c_glc_gs ! SFP: use send here since going from coupler to glc? !ic = c_glc_gr kArea = mct_aVect_indexRA(dom_g%data,afldname) - lSize = mct_avect_lSize(x2g_g) + !lSize = mct_avect_lSize(x2g_g) + lSize = mct_avect_lSize(x2gacc_g) !SFP: use accum flux vector !SFP:debug if(smb_counter==0)then !one day at 30 min land/atmos time steps write(logunit,*) ' ' - write(logunit,*) ' smb vector length (should be 7425) = ', lSize + write(logunit,*) ' smb vector length (7425 in coupler) = ', lSize write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) write(logunit,*) ' initial value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) @@ -1417,9 +1437,13 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) do n=1,lSize ca_g = dom_g%data%rAttr(kArea,n) - nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) + !nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) + nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2gacc_g%rAttr(index_x2g_Flgl_qice,n) !SFP: use accum flux vector end do + budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * 48.0d0 !SFP: hack to see if this recovers actual value from time averaged value + !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * x2gacc_gx_cnt !SFP: ideally use this or something like it to contain actual value + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice smb_vector_length = smb_vector_length +lSize @@ -1429,6 +1453,10 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) if(smb_counter==48)then !one day at 30 min land/atmos time steps write(logunit,*) ' ' write(logunit,*) ' smb_counter = ', smb_counter + write(logunit,*) ' x2gacc_gx_cnt = ', x2gacc_gx_cnt +! write(logunit,*) ' l2gacc_lx_cnt = ', l2gacc_lx_cnt +! write(logunit,*) ' current value of x2g_ vector = ', x2g_g%rAttr(index_x2g_Flgl_qice,:) +! write(logunit,*) ' current value of x2gacc_ vector = ', x2gacc_g%rAttr(index_x2g_Flgl_qice,:) write(logunit,*) ' final value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) write(logunit,*) ' ' end if From 0909a5af22474285c26e828a37f216d5f882336a Mon Sep 17 00:00:00 2001 From: Stephen Price Date: Thu, 19 Sep 2024 12:46:05 -0500 Subject: [PATCH 04/13] Checkpoint of working changes with debugging included. Checkpointing changes to coupler budget code that accounts for l2x_ and g2x_ Greenland surface mass balance fluxes, with various lines supporting debugging outputs included. Cleaned up code without debug lines to follow. --- .../mpas-albany-landice/driver/glc_comp_mct.F | 4 +- driver-mct/main/seq_diag_mct.F90 | 184 +++++++++--------- 2 files changed, 98 insertions(+), 90 deletions(-) diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F index 7bf72556603c..51e874c2e205 100644 --- a/components/mpas-albany-landice/driver/glc_comp_mct.F +++ b/components/mpas-albany-landice/driver/glc_comp_mct.F @@ -1482,8 +1482,8 @@ subroutine glc_export_mct(g2x_g, errorCode) !call route_ice_runoff(0.0_RKIND, & !Recuperate runoff routing switch code (originally in glc_route_ice_runoff module in earlier code), and attach to ice calving flux once present... ! rofi_to_ocn=Fogg_rofi, & ! rofi_to_ice=Figg_rofi) - !g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0!...and remove these placeholders - g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0001d0 !SFP: dummy value to see if passes through coupler + g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0!...and remove these placeholders + !g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0001d0 ! dummy value to allow tracking through coupler g2x_g % rAttr(index_g2x_Figg_rofi,n)=0.0 !...and remove these placeholders g2x_g % rAttr(index_g2x_Fogg_rofl,n) = 0.0 !Attach to subglacial liquid flux once present diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index af26aeacfbc2..556beae61752 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -48,9 +48,9 @@ module seq_diag_mct use shr_reprosum_mod, only : shr_reprosum_calc use seq_diagBGC_mct, only : seq_diagBGC_preprint_mct, seq_diagBGC_print_mct - use prep_glc_mod, only : prep_glc_get_x2gacc_gx_cnt !SFP: added this and next - unclear which is needed -! use prep_glc_mod, only : prep_glc_get_l2gacc_lx_cnt - + use prep_glc_mod, only : prep_glc_get_x2gacc_gx_cnt !SFP: added this and next + use glc_elevclass_mod, only: glc_get_num_elevation_classes + implicit none save private @@ -222,10 +222,6 @@ module seq_diag_mct logical :: flds_wiso ! If water isotope fields are active - !--- temporary pointers --- - integer , pointer :: x2gacc_gx_cnt ! SFP added -! integer , pointer :: l2gacc_lx_cnt ! SFP: unclear if this or the above is needed / more relevant - ! !PUBLIC DATA MEMBERS !--- time-averaged (annual?) global budge diagnostics --- @@ -274,8 +270,8 @@ module seq_diag_mct integer :: index_l2x_Flrl_irrig integer :: index_l2x_Flrl_wslake - integer :: index_l2x_Flgl_qice(0:10) !SFP added - integer :: index_x2l_Sg_ice_covered(0:10) !SFP added + integer, allocatable :: index_l2x_Flgl_qice(:) !SFP: added this and next; unclear if this is the best way to treat these + integer, allocatable :: index_x2l_Sg_ice_covered(:) integer :: index_x2l_Faxa_lwdn integer :: index_x2l_Faxa_rainc @@ -450,6 +446,9 @@ module seq_diag_mct integer :: index_x2i_Faxa_snow_18O integer :: index_x2i_Faxa_snow_HDO + integer :: glc_nec !SFP: added + integer :: x2gacc_gx_cnt ! SFP added (maybe not needed) + !=============================================================================== contains !=============================================================================== @@ -892,7 +891,7 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) real(r8) :: l2x_Flgl_qice_col_sum !SFP: sum of fluxes over no. MECs (cols) - character(len=64) :: name !SFP: added this and next 2 + character(len=64) :: name !SFP: added this and next 2 for support of working w/ data in MECs character(len= 2) :: cnum integer(in) :: num @@ -916,6 +915,13 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) kArea = mct_aVect_indexRA(dom_l%data,afldname) kl = mct_aVect_indexRA(frac_l,lfrinname) + ! get number of elevation classes and allocate relevant sets of indices + glc_nec = glc_get_num_elevation_classes() + if (first_time) then + allocate(index_l2x_Flgl_qice(0:glc_nec)) + allocate(index_x2l_Sg_ice_covered(0:glc_nec)) + end if + if (present(do_l2x)) then if (first_time) then index_l2x_Fall_swnet = mct_aVect_indexRA(l2x_l,'Fall_swnet') @@ -931,12 +937,13 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_irrig = mct_aVect_indexRA(l2x_l,'Flrl_irrig', perrWith='quiet') index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') - do num=0,10 !SFP: change later to 0,glc_nec_max (no. of elev classes) + !SFP: added this loop + do num=0,glc_nec write(cnum,'(i2.2)') num name = 'Flgl_qice' // cnum - index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) !SFP added + index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) name = 'Sg_ice_covered' // cnum - index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) !SFP added + index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) end do index_l2x_Fall_evap_16O = mct_aVect_indexRA(l2x_l,'Fall_evap_16O',perrWith='quiet') @@ -975,11 +982,12 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) l2x_Flgl_qice_col_sum = 0.0d0 - do num=0,10 !SFP: change later to 0,glc_nec_max (no elev classes) - !SFP: this somehow needs to allow for each of the 11 vectors associate w/ each of the 11 elev classes - l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) !SFP added + do num=0,glc_nec + !SFP: this should sum the contributions from each of the n vectors in the total no. of MECs + !SFP: product on RHS is the SMB flux times the fraction of area in that particular elevation class times the land cell area + l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * ca_l !SFP added end do - nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_Flgl_qice_col_sum !SFP added + nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - l2x_Flgl_qice_col_sum !SFP added if ( flds_wiso_lnd )then nf = f_wevap_16O; @@ -1014,10 +1022,14 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) end if end do -! budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice +! budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice !SFP: waiting for this to contain actual non-zero values budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP added + ! SFP: are these needed? Currently not sure how / if / when to deallocate these ... + !deallocate(index_l2x_Flgl_qice(0:glc_nec)) + !deallocate(index_x2l_Sg_ice_covered(0:glc_nec)) + end if if (present(do_x2l)) then @@ -1309,16 +1321,15 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) !----- local ----- type(mct_aVect), pointer :: g2x_g -! type(mct_aVect), pointer :: x2g_g - type(mct_aVect), pointer :: x2gacc_g !SFP: replace above w/ vector for accumulated fluxes + type(mct_aVect), pointer :: x2gacc_g type(mct_ggrid), pointer :: dom_g integer(in) :: n,ic,nf,ip ! generic index - integer(in) :: kArea ! index of area field in aVect - integer(in) :: lSize ! size of aVect - real(r8) :: ca_g ! area of a grid cell + integer(in) :: kArea ! index of area field in aVect + integer(in) :: lSize ! size of aVect + real(r8) :: ca_g ! area of a grid cell logical,save :: first_time = .true. - integer,save :: counter,smb_counter,calving_counter !SFP: for debugging + integer,save :: counter,smb_counter,calving_counter !SFP: added (mostly for debugging) integer,save :: smb_vector_length,calving_vector_length !----- formats ----- @@ -1334,11 +1345,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) dom_g => component_get_dom_cx(glc) g2x_g => component_get_c2x_cx(glc) -! x2g_g => component_get_x2c_cx(glc) - x2gacc_g => component_get_x2c_cx(glc) !SFP: use accum fluxes vector - - x2gacc_gx_cnt => prep_glc_get_x2gacc_gx_cnt() !SFP: counter for how many times SMB flux accumulation has occured -! l2gacc_lx_cnt => prep_glc_get l2gacc_lx_cnt() + x2gacc_g => component_get_x2c_cx(glc) if( present(do_g2x))then !SPF: glc to coupler @@ -1352,33 +1359,32 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) index_g2x_Figg_rofi = mct_aVect_indexRA(g2x_g,'Figg_rofi') !SFP:debug - write(logunit,*) ' ' - write(logunit,*) ' index_g2x_Fogg_rofl = ', index_g2x_Fogg_rofl - write(logunit,*) ' index_g2x_Fogg_rofi = ', index_g2x_Fogg_rofi - write(logunit,*) ' index_g2x_Figg_rofi = ', index_g2x_Figg_rofi - write(logunit,*) ' ' + !write(logunit,*) ' ' + !write(logunit,*) ' index_g2x_Fogg_rofl = ', index_g2x_Fogg_rofl + !write(logunit,*) ' index_g2x_Fogg_rofi = ', index_g2x_Fogg_rofi + !write(logunit,*) ' index_g2x_Figg_rofi = ', index_g2x_Figg_rofi + !write(logunit,*) ' ' end if - !ip = p_inst + !ip = p_inst !SFP: this value, day or inst, does not change anything here ip = p_day - !ic = c_glc_gs - ic = c_glc_gr !SFP: use recieve here since this is coming from glc to coupler? + ic = c_glc_gr !SFP: use recieve here ("_gr") since this is coming from glc to coupler? kArea = mct_aVect_indexRA(dom_g%data,afldname) lSize = mct_avect_lSize(g2x_g) !SFP:debug - if(calving_counter==0)then !one day at 30 min land/atmos time steps - write(logunit,*) ' ' - write(logunit,*) ' calving vector length (7425 in coupler) = ', lSize - write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) - write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) - write(logunit,*) ' intial value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) - write(logunit,*) ' calving flux to ocean (Fogg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Fogg_rofi,1) - write(logunit,*) ' calving flux to ice (Figg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Figg_rofi,1) - write(logunit,*) ' calving flux X area to ocean (Fogg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Fogg_rofi,1) - write(logunit,*) ' calving flux X area to ice (Figg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Figg_rofi,1) - end if + !if(calving_counter==0)then + !write(logunit,*) ' ' + !write(logunit,*) ' calving vector length (7425 in coupler) = ', lSize + !write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) + !write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) + !write(logunit,*) ' intial value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) + !write(logunit,*) ' calving flux to ocean (Fogg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Fogg_rofi,1) + !write(logunit,*) ' calving flux to ice (Figg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Figg_rofi,1) + !write(logunit,*) ' calving flux X area to ocean (Fogg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Fogg_rofi,1) + !write(logunit,*) ' calving flux X area to ice (Figg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Figg_rofi,1) + !end if do n=1,lSize ca_g = dom_g%data%rAttr(kArea,n) @@ -1389,60 +1395,63 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice - calving_vector_length = calving_vector_length +lSize - calving_counter = calving_counter + 1 + !SFP: only needed for debugging + !calving_vector_length = calving_vector_length +lSize + !calving_counter = calving_counter + 1 !SFP:debug - if(calving_counter==48)then !one day at 30 min land/atmos time steps - write(logunit,*) ' calving counter = ', calving_counter - write(logunit,*) ' final value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) - write(logunit,*) ' ' - end if + !if(calving_counter==48)then !one day at 30 min land/atmos time steps + !write(logunit,*) ' calving counter = ', calving_counter + !write(logunit,*) ' final value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) + !write(logunit,*) ' ' + !end if endif !SFP: end 'do_g2x' if( present(do_x2g))then !SFP: coupler to glc + x2gacc_gx_cnt = prep_glc_get_x2gacc_gx_cnt() !SFP: counter for how many times SMB flux accumulation has occured + ! note that this would be useful below but does not seem to work currently + ! (being reset to zero before being called here?) if (first_time) then - smb_counter=0 - smb_vector_length = 0 + smb_counter=0 !SFP: this may be needed in order to turn average flux into accumulated flux (by multiplying average by no of lnd coupling intervals) - !index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') - index_x2g_Flgl_qice = mct_aVect_indexRA(x2gacc_g,'Flgl_qice') !SFP: use accum flux vector + !smb_vector_length = 0 !SFP: debugging only + + index_x2g_Flgl_qice = mct_aVect_indexRA(x2gacc_g,'Flgl_qice') !SFP:debug - write(logunit,*) ' ' - write(logunit,*) ' index_x2g_Flgl_qice = ', index_x2g_Flgl_qice - write(logunit,*) ' ' + !write(logunit,*) ' ' + !write(logunit,*) ' index_x2g_Flgl_qice = ', index_x2g_Flgl_qice + !write(logunit,*) ' ' end if - !ip = p_inst + !ip = p_inst !SFP: as above, day vs. inst. does not seem to matter here ip = p_day - ic = c_glc_gs ! SFP: use send here since going from coupler to glc? - !ic = c_glc_gr + ic = c_glc_gs ! SFP: use send here ("_gs") since going from coupler to glc? kArea = mct_aVect_indexRA(dom_g%data,afldname) - !lSize = mct_avect_lSize(x2g_g) - lSize = mct_avect_lSize(x2gacc_g) !SFP: use accum flux vector + lSize = mct_avect_lSize(x2gacc_g) !SFP:debug - if(smb_counter==0)then !one day at 30 min land/atmos time steps - write(logunit,*) ' ' - write(logunit,*) ' smb vector length (7425 in coupler) = ', lSize - write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) - write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) - write(logunit,*) ' initial value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) - end if + !if(smb_counter==0)then !one day at 30 min land/atmos time steps + !write(logunit,*) ' ' + !write(logunit,*) ' smb vector length (7425 in coupler) = ', lSize + !write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) + !write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) + !write(logunit,*) ' initial value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) + !end if do n=1,lSize ca_g = dom_g%data%rAttr(kArea,n) - !nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) - nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2gacc_g%rAttr(index_x2g_Flgl_qice,n) !SFP: use accum flux vector + nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2gacc_g%rAttr(index_x2g_Flgl_qice,n) end do - budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * 48.0d0 !SFP: hack to see if this recovers actual value from time averaged value - !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * x2gacc_gx_cnt !SFP: ideally use this or something like it to contain actual value + !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * 48.0d0 !SFP: hack to see if this recovers actual value from time averaged value + !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * x2gacc_gx_cnt !SFP: ideally use this ... but always zero (zeroed before called?) + budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * smb_counter !SFP: this works but smb_counter seems like sloppy way to recover no. of lnd steps per glc coupling step + ! would be nicer to use value of x2gacc_gx_cnt (but always 0 as currently called) budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice @@ -1450,16 +1459,15 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) smb_counter = smb_counter + 1 !SFP:debug - if(smb_counter==48)then !one day at 30 min land/atmos time steps - write(logunit,*) ' ' - write(logunit,*) ' smb_counter = ', smb_counter - write(logunit,*) ' x2gacc_gx_cnt = ', x2gacc_gx_cnt -! write(logunit,*) ' l2gacc_lx_cnt = ', l2gacc_lx_cnt -! write(logunit,*) ' current value of x2g_ vector = ', x2g_g%rAttr(index_x2g_Flgl_qice,:) -! write(logunit,*) ' current value of x2gacc_ vector = ', x2gacc_g%rAttr(index_x2g_Flgl_qice,:) - write(logunit,*) ' final value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) - write(logunit,*) ' ' - end if + !if(smb_counter==48)then !one day at 30 min land/atmos time steps + !write(logunit,*) ' ' + !write(logunit,*) ' smb_counter = ', smb_counter + !write(logunit,*) ' x2gacc_gx_cnt = ', x2gacc_gx_cnt + !write(logunit,*) ' current value of x2g_ vector = ', x2g_g%rAttr(index_x2g_Flgl_qice,:) + !write(logunit,*) ' current value of x2gacc_ vector = ', x2gacc_g%rAttr(index_x2g_Flgl_qice,:) + !write(logunit,*) ' final value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) + !write(logunit,*) ' ' + !end if end if !SPF: end do coupler to glc From 8968cf34c28fde45a4f443a116f9c293c6f9c189 Mon Sep 17 00:00:00 2001 From: Stephen Price Date: Fri, 20 Sep 2024 13:58:43 -0500 Subject: [PATCH 05/13] Update glc copuler budgets for Grenland Clean up commenting and debugging lines, leaving bare minimum needed to make draft PR understandable. --- driver-mct/main/seq_diag_mct.F90 | 138 +++++++++---------------------- 1 file changed, 39 insertions(+), 99 deletions(-) diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 556beae61752..94446a3a3492 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -48,7 +48,7 @@ module seq_diag_mct use shr_reprosum_mod, only : shr_reprosum_calc use seq_diagBGC_mct, only : seq_diagBGC_preprint_mct, seq_diagBGC_print_mct - use prep_glc_mod, only : prep_glc_get_x2gacc_gx_cnt !SFP: added this and next + use prep_glc_mod, only : prep_glc_get_x2gacc_gx_cnt use glc_elevclass_mod, only: glc_get_num_elevation_classes implicit none @@ -143,13 +143,13 @@ module seq_diag_mct integer(in),parameter :: f_hsen =10 ! heat : sensible integer(in),parameter :: f_hpolar =11 ! heat : AIS imbalance integer(in),parameter :: f_hh2ot =12 ! heat : water temperature - integer(in),parameter :: f_hgsmb =13 ! heat : GIS SMB !SFP added + integer(in),parameter :: f_hgsmb =13 ! heat : Greenland ice sheet surface mass balance integer(in),parameter :: f_wfrz =14 ! water: freezing integer(in),parameter :: f_wmelt =15 ! water: melting integer(in),parameter :: f_wrain =16 ! water: precip, liquid integer(in),parameter :: f_wsnow =17 ! water: precip, frozen integer(in),parameter :: f_wpolar =18 ! water: AIS imbalance - integer(in),parameter :: f_wgsmb =19 ! water: GIS SMB !SFP added + integer(in),parameter :: f_wgsmb =19 ! water: Greenland ice sheet surface mass balance integer(in),parameter :: f_wevap =20 ! water: evaporation integer(in),parameter :: f_wroff =21 ! water: runoff/flood integer(in),parameter :: f_wioff =22 ! water: frozen runoff @@ -270,7 +270,7 @@ module seq_diag_mct integer :: index_l2x_Flrl_irrig integer :: index_l2x_Flrl_wslake - integer, allocatable :: index_l2x_Flgl_qice(:) !SFP: added this and next; unclear if this is the best way to treat these + integer, allocatable :: index_l2x_Flgl_qice(:) integer, allocatable :: index_x2l_Sg_ice_covered(:) integer :: index_x2l_Faxa_lwdn @@ -348,7 +348,7 @@ module seq_diag_mct integer :: index_g2x_Fogg_rofi integer :: index_g2x_Figg_rofi - integer :: index_x2g_Flgl_qice !SFP added + integer :: index_x2g_Flgl_qice integer :: index_x2o_Foxx_rofl_16O integer :: index_x2o_Foxx_rofi_16O @@ -446,8 +446,8 @@ module seq_diag_mct integer :: index_x2i_Faxa_snow_18O integer :: index_x2i_Faxa_snow_HDO - integer :: glc_nec !SFP: added - integer :: x2gacc_gx_cnt ! SFP added (maybe not needed) + integer :: glc_nec + integer :: x2gacc_gx_cnt !=============================================================================== contains @@ -881,7 +881,7 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) type(mct_aVect), pointer :: l2x_l ! model to drv bundle type(mct_aVect), pointer :: x2l_l ! drv to model bundle type(mct_ggrid), pointer :: dom_l - integer(in) :: n,ic,nf,ip ! generic index + integer(in) :: n,ic,nf,ip ! generic index integer(in) :: kArea ! index of area field in aVect integer(in) :: kl ! fraction indices integer(in) :: lSize ! size of aVect @@ -889,9 +889,9 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) logical,save :: first_time = .true. logical,save :: flds_wiso_lnd = .false. - real(r8) :: l2x_Flgl_qice_col_sum !SFP: sum of fluxes over no. MECs (cols) + real(r8) :: l2x_Flgl_qice_col_sum ! for summing fluxes over no. of elev. classes - character(len=64) :: name !SFP: added this and next 2 for support of working w/ data in MECs + character(len=64) :: name character(len= 2) :: cnum integer(in) :: num @@ -937,7 +937,6 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_irrig = mct_aVect_indexRA(l2x_l,'Flrl_irrig', perrWith='quiet') index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') - !SFP: added this loop do num=0,glc_nec write(cnum,'(i2.2)') num name = 'Flgl_qice' // cnum @@ -979,15 +978,15 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) if (index_l2x_Flrl_irrig /= 0) then nf = f_wroff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_irrig,n) end if - nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) + nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) ! contribution from land ice calving currently zero l2x_Flgl_qice_col_sum = 0.0d0 do num=0,glc_nec - !SFP: this should sum the contributions from each of the n vectors in the total no. of MECs - !SFP: product on RHS is the SMB flux times the fraction of area in that particular elevation class times the land cell area - l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * ca_l !SFP added + ! sums the contributions from fluxes in each set of elevation classes + ! RHS product is flux times fraction of area in specific elevation class times land cell area + l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * ca_l end do - nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - l2x_Flgl_qice_col_sum !SFP added + nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - l2x_Flgl_qice_col_sum if ( flds_wiso_lnd )then nf = f_wevap_16O; @@ -1022,11 +1021,10 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) end if end do -! budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice !SFP: waiting for this to contain actual non-zero values - - budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice !SFP added + budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice ! contribution from land ice calving currently zero + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice - ! SFP: are these needed? Currently not sure how / if / when to deallocate these ... + ! Nneeded? Not sure if / when these should be deallocated !deallocate(index_l2x_Flgl_qice(0:glc_nec)) !deallocate(index_x2l_Sg_ice_covered(0:glc_nec)) @@ -1305,14 +1303,14 @@ end subroutine seq_diag_rof_mct ! Compute global glc input/output flux diagnostics ! ! !REVISION HISTORY: - ! 2008-jul-10 - T. Craig - update + ! 2024-Sept. - S. Price - update ! ! !INTERFACE: ------------------------------------------------------------------ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) type(component_type) , intent(in) :: glc ! component type for instance1 - type(mct_aVect) , intent(in) :: frac_g ! frac bundle !SFP: does not look like fractions are needed / used here? + type(mct_aVect) , intent(in) :: frac_g ! frac bundle (may not be used / needed here) type(seq_infodata_type) , intent(in) :: infodata logical , intent(in), optional :: do_x2g logical , intent(in), optional :: do_g2x @@ -1329,7 +1327,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) real(r8) :: ca_g ! area of a grid cell logical,save :: first_time = .true. - integer,save :: counter,smb_counter,calving_counter !SFP: added (mostly for debugging) + integer,save :: counter,smb_counter,calving_counter ! SFP: Debugging integer,save :: smb_vector_length,calving_vector_length !----- formats ----- @@ -1347,7 +1345,9 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) g2x_g => component_get_c2x_cx(glc) x2gacc_g => component_get_x2c_cx(glc) - if( present(do_g2x))then !SPF: glc to coupler + ip = p_inst + + if( present(do_g2x))then ! do fields from glc to coupler (g2x_) if (first_time) then @@ -1358,34 +1358,12 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) index_g2x_Fogg_rofi = mct_aVect_indexRA(g2x_g,'Fogg_rofi') index_g2x_Figg_rofi = mct_aVect_indexRA(g2x_g,'Figg_rofi') - !SFP:debug - !write(logunit,*) ' ' - !write(logunit,*) ' index_g2x_Fogg_rofl = ', index_g2x_Fogg_rofl - !write(logunit,*) ' index_g2x_Fogg_rofi = ', index_g2x_Fogg_rofi - !write(logunit,*) ' index_g2x_Figg_rofi = ', index_g2x_Figg_rofi - !write(logunit,*) ' ' - end if - !ip = p_inst !SFP: this value, day or inst, does not change anything here - ip = p_day - ic = c_glc_gr !SFP: use recieve here ("_gr") since this is coming from glc to coupler? + ic = c_glc_gr kArea = mct_aVect_indexRA(dom_g%data,afldname) lSize = mct_avect_lSize(g2x_g) - !SFP:debug - !if(calving_counter==0)then - !write(logunit,*) ' ' - !write(logunit,*) ' calving vector length (7425 in coupler) = ', lSize - !write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) - !write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) - !write(logunit,*) ' intial value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) - !write(logunit,*) ' calving flux to ocean (Fogg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Fogg_rofi,1) - !write(logunit,*) ' calving flux to ice (Figg_rofi) in one cell = ', g2x_g%rAttr(index_g2x_Figg_rofi,1) - !write(logunit,*) ' calving flux X area to ocean (Fogg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Fogg_rofi,1) - !write(logunit,*) ' calving flux X area to ice (Figg_rofi) in one cell = ', dom_g%data%rAttr(kArea,1)*g2x_g%rAttr(index_g2x_Figg_rofi,1) - !end if - do n=1,lSize ca_g = dom_g%data%rAttr(kArea,n) nf = f_wroff; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_g*g2x_g%rAttr(index_g2x_Fogg_rofl,n) @@ -1395,81 +1373,43 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) budg_dataL(f_hioff,ic,ip) = -budg_dataL(f_wioff,ic,ip)*shr_const_latice - !SFP: only needed for debugging - !calving_vector_length = calving_vector_length +lSize - !calving_counter = calving_counter + 1 + endif ! end do fields from glc to coupler (g2x_) - !SFP:debug - !if(calving_counter==48)then !one day at 30 min land/atmos time steps - !write(logunit,*) ' calving counter = ', calving_counter - !write(logunit,*) ' final value of calving flux sum vector = ', budg_dataL(f_wioff,ic,ip) - !write(logunit,*) ' ' - !end if + if( present(do_x2g))then ! do fields from coupler to glc (x2g_) - endif !SFP: end 'do_g2x' - - if( present(do_x2g))then !SFP: coupler to glc - - x2gacc_gx_cnt = prep_glc_get_x2gacc_gx_cnt() !SFP: counter for how many times SMB flux accumulation has occured + x2gacc_gx_cnt = prep_glc_get_x2gacc_gx_cnt() ! counter for how many times SMB flux accumulation has occured ! note that this would be useful below but does not seem to work currently ! (being reset to zero before being called here?) if (first_time) then - smb_counter=0 !SFP: this may be needed in order to turn average flux into accumulated flux (by multiplying average by no of lnd coupling intervals) - - !smb_vector_length = 0 !SFP: debugging only - - index_x2g_Flgl_qice = mct_aVect_indexRA(x2gacc_g,'Flgl_qice') - - !SFP:debug - !write(logunit,*) ' ' - !write(logunit,*) ' index_x2g_Flgl_qice = ', index_x2g_Flgl_qice - !write(logunit,*) ' ' + smb_counter=0 ! something like this (or above) needed to turn average flux + ! into accumulated flux (i.e., multiply average flux by no. of lnd coupling intervals) + index_x2g_Flgl_qice = mct_aVect_indexRA(x2gacc_g,'Flgl_qice') ! While name suggests this holds accumulated flux, + ! it appears to actually be the average flux (e.g. see + ! subroutine 'prep_glc_accum_avg' in prep_glc_mod.f90. + ! (also note that this same value gets copied to x2g_) end if - !ip = p_inst !SFP: as above, day vs. inst. does not seem to matter here - ip = p_day - ic = c_glc_gs ! SFP: use send here ("_gs") since going from coupler to glc? + ic = c_glc_gs kArea = mct_aVect_indexRA(dom_g%data,afldname) lSize = mct_avect_lSize(x2gacc_g) - !SFP:debug - !if(smb_counter==0)then !one day at 30 min land/atmos time steps - !write(logunit,*) ' ' - !write(logunit,*) ' smb vector length (7425 in coupler) = ', lSize - !write(logunit,*) ' kArea(1) = ', dom_g%data%rAttr(kArea,1) - !write(logunit,*) ' kArea(50) = ', dom_g%data%rAttr(kArea,50) - !write(logunit,*) ' initial value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) - !end if - do n=1,lSize ca_g = dom_g%data%rAttr(kArea,n) nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2gacc_g%rAttr(index_x2g_Flgl_qice,n) end do - !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * 48.0d0 !SFP: hack to see if this recovers actual value from time averaged value - !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * x2gacc_gx_cnt !SFP: ideally use this ... but always zero (zeroed before called?) - budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * smb_counter !SFP: this works but smb_counter seems like sloppy way to recover no. of lnd steps per glc coupling step - ! would be nicer to use value of x2gacc_gx_cnt (but always 0 as currently called) + !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * x2gacc_gx_cnt ! ideally use something like this for multiplying average flux + ! to get accumulated flux (but currently always zero) + budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * smb_counter ! works for now, but sloppy and only works for a 1 day run budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice smb_vector_length = smb_vector_length +lSize smb_counter = smb_counter + 1 - !SFP:debug - !if(smb_counter==48)then !one day at 30 min land/atmos time steps - !write(logunit,*) ' ' - !write(logunit,*) ' smb_counter = ', smb_counter - !write(logunit,*) ' x2gacc_gx_cnt = ', x2gacc_gx_cnt - !write(logunit,*) ' current value of x2g_ vector = ', x2g_g%rAttr(index_x2g_Flgl_qice,:) - !write(logunit,*) ' current value of x2gacc_ vector = ', x2gacc_g%rAttr(index_x2g_Flgl_qice,:) - !write(logunit,*) ' final value of smb sum vector = ', budg_dataL(f_wgsmb,ic,ip) - !write(logunit,*) ' ' - !end if - - end if !SPF: end do coupler to glc + end if ! end do fields from coupler to glc (x2g_) first_time = .false. From 1e3f3b9e7c2fc6b2e7ca469ef74941ce8be3a065 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 25 Nov 2024 15:21:48 -0600 Subject: [PATCH 06/13] Add call to zero out x2g_gx if the accum counter is 1 --- driver-mct/main/cime_comp_mod.F90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/driver-mct/main/cime_comp_mod.F90 b/driver-mct/main/cime_comp_mod.F90 index f9678b3e7c9d..687050464a36 100644 --- a/driver-mct/main/cime_comp_mod.F90 +++ b/driver-mct/main/cime_comp_mod.F90 @@ -2523,6 +2523,7 @@ subroutine cime_run() logical :: lnd2glc_averaged_now ! Whether lnd2glc averages were taken this timestep logical :: prep_glc_accum_avg_called ! Whether prep_glc_accum_avg has been called this timestep integer :: i, nodeId + integer :: l2gacc_lx_cnt character(len=15) :: c_ymdtod character(len=18) :: c_mprof_file @@ -3047,6 +3048,14 @@ subroutine cime_run() !---------------------------------------------------------- !| GLC SETUP-SEND !---------------------------------------------------------- + ! zero out x2g_gx if this is the first call to prep_glc_accum_avg + if (glc_present) then + l2gacc_lx_cnt = prep_glc_get_l2gacc_lx_cnt() + if (l2gacc_lx_cnt.eq.1) then + call prep_glc_zero_fields() + endif + endif + if (glc_present .and. glcrun_alarm) then call cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_called) endif @@ -3095,6 +3104,7 @@ subroutine cime_run() endif endif + !---------------------------------------------------------- !| Budget with old fractions !---------------------------------------------------------- @@ -4745,7 +4755,7 @@ subroutine cime_run_calc_budgets1(in_cplrun) call seq_diag_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_x2i=.true.) endif if (glc_present) then - !call seq_diag_glc_mct(glc(ens1), fractions_gx(ens1), infodata, do_x2g=.true., do_g2x=.true.) !SFP: comment out for now while debugging + call seq_diag_glc_mct(glc(ens1), fractions_gx(ens1), infodata, do_x2g=.true.) endif if (do_bgc_budgets) then if (rof_present) then @@ -4787,7 +4797,7 @@ subroutine cime_run_calc_budgets2(in_cplrun) call seq_diag_ice_mct(ice(ens1), fractions_ix(ens1), infodata, do_i2x=.true.) endif if (glc_present) then - call seq_diag_glc_mct(glc(ens1), fractions_gx(ens1), infodata, do_x2g=.true., do_g2x=.true.) + call seq_diag_glc_mct(glc(ens1), fractions_gx(ens1), infodata, do_g2x=.true.) endif if (do_bgc_budgets) then if (atm_present) then @@ -5596,3 +5606,4 @@ function copy_and_trim_rpointer_file(src, dst) result(out) end function copy_and_trim_rpointer_file end module cime_comp_mod + From 61c49ba2408aeaa2187d3beed5686b330eab182d Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 25 Nov 2024 15:24:32 -0600 Subject: [PATCH 07/13] Add l2gacc_lx_cnt_avg to count number of times used in averaging --- driver-mct/main/prep_glc_mod.F90 | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/driver-mct/main/prep_glc_mod.F90 b/driver-mct/main/prep_glc_mod.F90 index 07aeb9890bda..d2b309e53740 100644 --- a/driver-mct/main/prep_glc_mod.F90 +++ b/driver-mct/main/prep_glc_mod.F90 @@ -45,6 +45,7 @@ module prep_glc_mod public :: prep_glc_get_l2gacc_lx public :: prep_glc_get_l2gacc_lx_one_instance public :: prep_glc_get_l2gacc_lx_cnt + public :: prep_glc_get_l2gacc_lx_cnt_avg public :: prep_glc_get_o2x_gx public :: prep_glc_get_x2gacc_gx @@ -91,6 +92,7 @@ module prep_glc_mod type(mct_aVect), pointer :: l2gacc_lx(:) ! Lnd export, lnd grid, cpl pes - allocated in driver integer , target :: l2gacc_lx_cnt ! l2gacc_lx: number of time samples accumulated + integer , target :: l2gacc_lx_cnt_avg ! l2gacc_lx: number of time samples averaged ! other module variables integer :: mpicom_CPLID ! MPI cpl communicator @@ -195,6 +197,7 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) call mct_aVect_zero(l2gacc_lx(eli)) end do l2gacc_lx_cnt = 0 + l2gacc_lx_cnt_avg = 0 end if if (glc_present .and. lnd_c2_glc) then @@ -502,6 +505,7 @@ subroutine prep_glc_accum_avg(timer, lnd2glc_averaged_now) call mct_avect_avg(l2gacc_lx(eli), l2gacc_lx_cnt) end do end if + l2gacc_lx_cnt_avg = l2gacc_lx_cnt l2gacc_lx_cnt = 0 ! Accumulation for OCN @@ -950,6 +954,7 @@ subroutine prep_glc_zero_fields() type(mct_avect), pointer :: x2g_gx !--------------------------------------------------------------- + do egi = 1,num_inst_glc x2g_gx => component_get_x2c_cx(glc(egi)) call mct_aVect_zero(x2g_gx) @@ -1424,6 +1429,11 @@ function prep_glc_get_l2gacc_lx_cnt() prep_glc_get_l2gacc_lx_cnt => l2gacc_lx_cnt end function prep_glc_get_l2gacc_lx_cnt + function prep_glc_get_l2gacc_lx_cnt_avg() + integer, pointer :: prep_glc_get_l2gacc_lx_cnt_avg + prep_glc_get_l2gacc_lx_cnt_avg => l2gacc_lx_cnt_avg + end function prep_glc_get_l2gacc_lx_cnt_avg + function prep_glc_get_o2x_gx() type(mct_aVect), pointer :: prep_glc_get_o2x_gx(:) prep_glc_get_o2x_gx => o2x_gx(:) From 8461481114698dfd3633aa1d4f32fae94ea30303 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 25 Nov 2024 15:27:48 -0600 Subject: [PATCH 08/13] Use l2gacc_lx_cnt_avg for calculating some glc budget terms --- driver-mct/main/seq_diag_mct.F90 | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 94446a3a3492..425ab1be5176 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -48,7 +48,7 @@ module seq_diag_mct use shr_reprosum_mod, only : shr_reprosum_calc use seq_diagBGC_mct, only : seq_diagBGC_preprint_mct, seq_diagBGC_print_mct - use prep_glc_mod, only : prep_glc_get_x2gacc_gx_cnt + use prep_glc_mod, only : prep_glc_get_l2gacc_lx_cnt_avg use glc_elevclass_mod, only: glc_get_num_elevation_classes implicit none @@ -447,7 +447,7 @@ module seq_diag_mct integer :: index_x2i_Faxa_snow_HDO integer :: glc_nec - integer :: x2gacc_gx_cnt + integer :: l2gacc_lx_cnt_avg !=============================================================================== contains @@ -1319,7 +1319,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) !----- local ----- type(mct_aVect), pointer :: g2x_g - type(mct_aVect), pointer :: x2gacc_g + type(mct_aVect), pointer :: x2g_g type(mct_ggrid), pointer :: dom_g integer(in) :: n,ic,nf,ip ! generic index integer(in) :: kArea ! index of area field in aVect @@ -1327,7 +1327,6 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) real(r8) :: ca_g ! area of a grid cell logical,save :: first_time = .true. - integer,save :: counter,smb_counter,calving_counter ! SFP: Debugging integer,save :: smb_vector_length,calving_vector_length !----- formats ----- @@ -1343,7 +1342,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) dom_g => component_get_dom_cx(glc) g2x_g => component_get_c2x_cx(glc) - x2gacc_g => component_get_x2c_cx(glc) + x2g_g => component_get_x2c_cx(glc) ip = p_inst @@ -1351,7 +1350,6 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) if (first_time) then - calving_counter=0 calving_vector_length = 0 index_g2x_Fogg_rofl = mct_aVect_indexRA(g2x_g,'Fogg_rofl') @@ -1377,37 +1375,27 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) if( present(do_x2g))then ! do fields from coupler to glc (x2g_) - x2gacc_gx_cnt = prep_glc_get_x2gacc_gx_cnt() ! counter for how many times SMB flux accumulation has occured - ! note that this would be useful below but does not seem to work currently - ! (being reset to zero before being called here?) if (first_time) then - smb_counter=0 ! something like this (or above) needed to turn average flux - ! into accumulated flux (i.e., multiply average flux by no. of lnd coupling intervals) + index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') - index_x2g_Flgl_qice = mct_aVect_indexRA(x2gacc_g,'Flgl_qice') ! While name suggests this holds accumulated flux, - ! it appears to actually be the average flux (e.g. see - ! subroutine 'prep_glc_accum_avg' in prep_glc_mod.f90. - ! (also note that this same value gets copied to x2g_) end if + l2gacc_lx_cnt_avg = prep_glc_get_l2gacc_lx_cnt_avg() ! counter for how many times SMB flux accumulation has occured ic = c_glc_gs kArea = mct_aVect_indexRA(dom_g%data,afldname) - lSize = mct_avect_lSize(x2gacc_g) + lSize = mct_avect_lSize(x2g_g) do n=1,lSize ca_g = dom_g%data%rAttr(kArea,n) - nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2gacc_g%rAttr(index_x2g_Flgl_qice,n) + nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) end do - !budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * x2gacc_gx_cnt ! ideally use something like this for multiplying average flux - ! to get accumulated flux (but currently always zero) - budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * smb_counter ! works for now, but sloppy and only works for a 1 day run + budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) * l2gacc_lx_cnt_avg budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice smb_vector_length = smb_vector_length +lSize - smb_counter = smb_counter + 1 end if ! end do fields from coupler to glc (x2g_) From 7bd74c698287050a42113790cae376fb7bf09b37 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 25 Nov 2024 20:35:18 -0600 Subject: [PATCH 09/13] Add logic to make sure glc_nec is at least one --- driver-mct/main/seq_diag_mct.F90 | 36 +++++++++++++++++++------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 425ab1be5176..00ef16dcd77f 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -917,9 +917,11 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) ! get number of elevation classes and allocate relevant sets of indices glc_nec = glc_get_num_elevation_classes() - if (first_time) then - allocate(index_l2x_Flgl_qice(0:glc_nec)) - allocate(index_x2l_Sg_ice_covered(0:glc_nec)) + if (glc_nec.ge.1) then + if (first_time) then + allocate(index_l2x_Flgl_qice(0:glc_nec)) + allocate(index_x2l_Sg_ice_covered(0:glc_nec)) + end if end if if (present(do_l2x)) then @@ -937,13 +939,15 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_irrig = mct_aVect_indexRA(l2x_l,'Flrl_irrig', perrWith='quiet') index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') - do num=0,glc_nec - write(cnum,'(i2.2)') num - name = 'Flgl_qice' // cnum - index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) - name = 'Sg_ice_covered' // cnum - index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) - end do + if (glc_nec.ge.1) then + do num=0,glc_nec + write(cnum,'(i2.2)') num + name = 'Flgl_qice' // cnum + index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) + name = 'Sg_ice_covered' // cnum + index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) + end do + end if index_l2x_Fall_evap_16O = mct_aVect_indexRA(l2x_l,'Fall_evap_16O',perrWith='quiet') if ( index_l2x_Fall_evap_16O /= 0 ) flds_wiso_lnd = .true. @@ -981,11 +985,13 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) nf = f_wioff ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_l*l2x_l%rAttr(index_l2x_Flrl_rofi,n) ! contribution from land ice calving currently zero l2x_Flgl_qice_col_sum = 0.0d0 - do num=0,glc_nec - ! sums the contributions from fluxes in each set of elevation classes - ! RHS product is flux times fraction of area in specific elevation class times land cell area - l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * ca_l - end do + if (glc_nec.ge.1) then + do num=0,glc_nec + ! sums the contributions from fluxes in each set of elevation classes + ! RHS product is flux times fraction of area in specific elevation class times land cell area + l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * ca_l + end do + end if nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - l2x_Flgl_qice_col_sum if ( flds_wiso_lnd )then From 31d46ed3414f9f0944977572ea7e2e17f95066e0 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Mon, 25 Nov 2024 20:39:19 -0600 Subject: [PATCH 10/13] Fix trailing white space --- driver-mct/main/seq_diag_mct.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 00ef16dcd77f..26d21a2879a6 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -940,12 +940,12 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') if (glc_nec.ge.1) then - do num=0,glc_nec + do num=0,glc_nec write(cnum,'(i2.2)') num name = 'Flgl_qice' // cnum index_l2x_Flgl_qice(num) = mct_avect_indexRA(l2x_l,trim(name)) name = 'Sg_ice_covered' // cnum - index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) + index_x2l_Sg_ice_covered(num) = mct_avect_indexRA(x2l_l,trim(name)) end do end if From 822cbef84ce423c6828fbeb61fe48d62d4e593fd Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Tue, 26 Nov 2024 15:39:30 -0600 Subject: [PATCH 11/13] Better way to handle call to prep_glc_zero_fields --- driver-mct/main/cime_comp_mod.F90 | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/driver-mct/main/cime_comp_mod.F90 b/driver-mct/main/cime_comp_mod.F90 index 687050464a36..5ec98de349ef 100644 --- a/driver-mct/main/cime_comp_mod.F90 +++ b/driver-mct/main/cime_comp_mod.F90 @@ -2523,7 +2523,6 @@ subroutine cime_run() logical :: lnd2glc_averaged_now ! Whether lnd2glc averages were taken this timestep logical :: prep_glc_accum_avg_called ! Whether prep_glc_accum_avg has been called this timestep integer :: i, nodeId - integer :: l2gacc_lx_cnt character(len=15) :: c_ymdtod character(len=18) :: c_mprof_file @@ -3048,18 +3047,14 @@ subroutine cime_run() !---------------------------------------------------------- !| GLC SETUP-SEND !---------------------------------------------------------- - ! zero out x2g_gx if this is the first call to prep_glc_accum_avg if (glc_present) then - l2gacc_lx_cnt = prep_glc_get_l2gacc_lx_cnt() - if (l2gacc_lx_cnt.eq.1) then + if (glcrun_alarm) then + call cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_called) + else call prep_glc_zero_fields() endif endif - if (glc_present .and. glcrun_alarm) then - call cime_run_glc_setup_send(lnd2glc_averaged_now, prep_glc_accum_avg_called) - endif - ! ------------------------------------------------------------------------ ! Also average lnd2glc fields if needed for requested l2x1yrg auxiliary history ! files, even if running with a stub glc model. From 7b42c2e8ee1c0230da2578875ccda882353d7fe1 Mon Sep 17 00:00:00 2001 From: Jon Wolfe Date: Thu, 5 Dec 2024 10:31:06 -0600 Subject: [PATCH 12/13] Updates for the glc budget --- driver-mct/main/prep_glc_mod.F90 | 5 ++++- driver-mct/main/seq_diag_mct.F90 | 13 ++++++++++--- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/driver-mct/main/prep_glc_mod.F90 b/driver-mct/main/prep_glc_mod.F90 index d2b309e53740..f818a9c2a14d 100644 --- a/driver-mct/main/prep_glc_mod.F90 +++ b/driver-mct/main/prep_glc_mod.F90 @@ -1200,8 +1200,9 @@ subroutine prep_glc_renormalize_smb(eli, fractions_lx, g2x_gx, mapper_Fg2l, area aream_l(:) = dom_l%data%rAttr(km,:) ! Export land fractions from fractions_lx to a local array + ! Note that for E3SM we are using lfrin instead of lfrac allocate(lfrac(lsize_l)) - call mct_aVect_exportRattr(fractions_lx, "lfrac", lfrac) + call mct_aVect_exportRattr(fractions_lx, "lfrin", lfrac) ! Map Sg_icemask from the glc grid to the land grid. ! This may not be necessary, if Sg_icemask_l has already been mapped from Sg_icemask_g. @@ -1384,6 +1385,8 @@ subroutine prep_glc_renormalize_smb(eli, fractions_lx, g2x_gx, mapper_Fg2l, area endif if (iamroot) then + write(logunit,*) 'global_accum_on_land_grid = ', global_accum_on_land_grid + write(logunit,*) 'global_accum_on_glc_grid = ', global_accum_on_glc_grid write(logunit,*) 'accum_renorm_factor = ', accum_renorm_factor write(logunit,*) 'ablat_renorm_factor = ', ablat_renorm_factor endif diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 26d21a2879a6..6177ce5db192 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -270,6 +270,7 @@ module seq_diag_mct integer :: index_l2x_Flrl_irrig integer :: index_l2x_Flrl_wslake + integer :: index_x2l_Sg_icemask integer, allocatable :: index_l2x_Flgl_qice(:) integer, allocatable :: index_x2l_Sg_ice_covered(:) @@ -349,6 +350,7 @@ module seq_diag_mct integer :: index_g2x_Figg_rofi integer :: index_x2g_Flgl_qice + integer :: index_g2x_Sg_icemask integer :: index_x2o_Foxx_rofl_16O integer :: index_x2o_Foxx_rofi_16O @@ -890,6 +892,7 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) logical,save :: flds_wiso_lnd = .false. real(r8) :: l2x_Flgl_qice_col_sum ! for summing fluxes over no. of elev. classes + real(r8) :: effective_area character(len=64) :: name character(len= 2) :: cnum @@ -940,6 +943,7 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) index_l2x_Flrl_wslake = mct_aVect_indexRA(l2x_l,'Flrl_wslake') if (glc_nec.ge.1) then + index_x2l_Sg_icemask = mct_avect_indexRA(x2l_l,'Sg_icemask') do num=0,glc_nec write(cnum,'(i2.2)') num name = 'Flgl_qice' // cnum @@ -986,10 +990,12 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) l2x_Flgl_qice_col_sum = 0.0d0 if (glc_nec.ge.1) then + effective_area = min(frac_l%rAttr(kl,n),x2l_l%rAttr(index_x2l_Sg_icemask,n)) * dom_l%data%rAttr(kArea,n) do num=0,glc_nec ! sums the contributions from fluxes in each set of elevation classes ! RHS product is flux times fraction of area in specific elevation class times land cell area - l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * ca_l + l2x_Flgl_qice_col_sum = l2x_Flgl_qice_col_sum + l2x_l%rAttr(index_l2x_Flgl_qice(num),n) * & + x2l_l%rAttr(index_x2l_Sg_ice_covered(num),n) * effective_area end do end if nf = f_wgsmb ; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - l2x_Flgl_qice_col_sum @@ -1383,7 +1389,8 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) if (first_time) then - index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') + index_x2g_Flgl_qice = mct_aVect_indexRA(x2g_g,'Flgl_qice') + index_g2x_Sg_icemask = mct_avect_indexRA(g2x_g,'Sg_icemask') end if @@ -1393,7 +1400,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) lSize = mct_avect_lSize(x2g_g) do n=1,lSize - ca_g = dom_g%data%rAttr(kArea,n) + ca_g = dom_g%data%rAttr(kArea,n)*g2x_g%rAttr(index_g2x_Sg_icemask,n) nf = f_wgsmb; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_g*x2g_g%rAttr(index_x2g_Flgl_qice,n) end do From 2a630b10ce3a0f27544d6cfffd5f49a5a139dc49 Mon Sep 17 00:00:00 2001 From: Stephen Price Date: Fri, 6 Dec 2024 17:39:33 -0600 Subject: [PATCH 13/13] Minor cleanup Remove placeholder / commmented out code in glc_comp_mct Correct revision history comment in glc section of seq_diag_mct --- components/mpas-albany-landice/driver/glc_comp_mct.F | 5 +++-- driver-mct/main/seq_diag_mct.F90 | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F index f6001a73e621..dd254de9807e 100644 --- a/components/mpas-albany-landice/driver/glc_comp_mct.F +++ b/components/mpas-albany-landice/driver/glc_comp_mct.F @@ -1490,11 +1490,12 @@ subroutine glc_export_mct(g2x_g, errorCode) do i = 1, nCellsSolve n = n + 1 - !call route_ice_runoff(0.0_RKIND, & !Recuperate runoff routing switch code (originally in glc_route_ice_runoff module in earlier code), and attach to ice calving flux once present... + ! Recuperate runoff routing switch code (originally in glc_route_ice_runoff module in earlier code), + ! and attach to ice calving flux once present... + !call route_ice_runoff(0.0_RKIND, & ! rofi_to_ocn=Fogg_rofi, & ! rofi_to_ice=Figg_rofi) g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0!...and remove these placeholders - !g2x_g % rAttr(index_g2x_Fogg_rofi,n)=0.0001d0 ! dummy value to allow tracking through coupler g2x_g % rAttr(index_g2x_Figg_rofi,n)=0.0 !...and remove these placeholders g2x_g % rAttr(index_g2x_Fogg_rofl,n) = 0.0 !Attach to subglacial liquid flux once present diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 6177ce5db192..be9c3a1426ce 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -1315,7 +1315,8 @@ end subroutine seq_diag_rof_mct ! Compute global glc input/output flux diagnostics ! ! !REVISION HISTORY: - ! 2024-Sept. - S. Price - update + ! 2008-Jul-10 - T. Craig - update + ! 2024-Dec-06 - S. Price, J. Wolfe - update ! ! !INTERFACE: ------------------------------------------------------------------