diff --git a/cime_config/tests.py b/cime_config/tests.py index c1f6cfaa1b70..12ab829c51dd 100644 --- a/cime_config/tests.py +++ b/cime_config/tests.py @@ -274,6 +274,7 @@ "ERS_P480_Ld5.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF.mpaso-jra_1958", "PEM_P480_Ld5.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF.mpaso-jra_1958", "SMS_P480_Ld5.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF-TMIX.mpaso-jra_1958", + "SMS_P480_Ld5.TL319_IcoswISC30E3r5.GMPAS-JRA1p5-DIB-PISMF-DSGR.mpaso-jra_1958", ) }, diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 1c55bcb7be1d..37c8b88f67e6 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -54,6 +54,9 @@ OPTIONS -ocn_ismf variable for defining how the ocn model will handle ice shelf melt fluxes Options are: none, data, internal, coupled + -ocn_sgr variable for defining how the ocn model will handle subglacial + runoff + Options are: none, data -decomp_prefix decomp_prefix variable -date_stamp date_stamp variable -cfg_grid Directory containing MPASO configuration scripts. @@ -111,6 +114,7 @@ my %opts = ( help => 0, ocn_forcing => undef, ocn_iceberg => undef, ocn_ismf => undef, + ocn_sgr => undef, decomp_prefix => undef, date_stamp => undef, ocn_bgc => undef, @@ -137,6 +141,7 @@ GetOptions( "ocn_forcing=s" => \$opts{'ocn_forcing'}, "ocn_iceberg=s" => \$opts{'ocn_iceberg'}, "ocn_ismf=s" => \$opts{'ocn_ismf'}, + "ocn_sgr=s" => \$opts{'ocn_sgr'}, "decomp_prefix=s" => \$opts{'decomp_prefix'}, "date_stamp=s" => \$opts{'date_stamp'}, "ocn_bgc=s" => \$opts{'ocn_bgc'}, @@ -180,6 +185,7 @@ my $OCN_GRID = $opts{'ocn_grid'}; my $OCN_FORCING = $opts{'ocn_forcing'}; my $OCN_ICEBERG = $opts{'ocn_iceberg'}; my $OCN_ISMF = $opts{'ocn_ismf'}; +my $OCN_SGR = $opts{'ocn_sgr'}; my $decomp_prefix = $opts{'decomp_prefix'}; my $date_stamp = $opts{'date_stamp'}; my $ocn_bgc = $opts{'ocn_bgc'}; @@ -702,6 +708,18 @@ add_default($nl, 'config_use_bulk_wind_stress'); add_default($nl, 'config_use_bulk_thickness_flux'); add_default($nl, 'config_flux_attenuation_coefficient'); add_default($nl, 'config_flux_attenuation_coefficient_runoff'); +if ($OCN_SGR eq 'data') { + add_default($nl, 'config_subglacial_runoff_mode', 'val'=>"data"); +} else { + add_default($nl, 'config_subglacial_runoff_mode'); +} +add_default($nl, 'config_flux_attenuation_coefficient_subglacial_runoff'); +add_default($nl, 'config_sgr_flux_vertical_location'); +add_default($nl, 'config_use_sgr_opt_kpp'); +add_default($nl, 'config_use_sgr_opt_temp_prescribed'); +add_default($nl, 'config_use_sgr_opt_salt_prescribed'); +add_default($nl, 'config_sgr_temperature_prescribed'); +add_default($nl, 'config_sgr_salinity_prescribed'); ############################ # Namelist group: coupling # diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index cc1695243faa..9e646beb6b5e 100644 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -224,6 +224,14 @@ add_default($nl, 'config_use_bulk_wind_stress'); add_default($nl, 'config_use_bulk_thickness_flux'); add_default($nl, 'config_flux_attenuation_coefficient'); add_default($nl, 'config_flux_attenuation_coefficient_runoff'); +add_default($nl, 'config_subglacial_runoff_mode'); +add_default($nl, 'config_flux_attenuation_coefficient_subglacial_runoff'); +add_default($nl, 'config_sgr_flux_vertical_location'); +add_default($nl, 'config_use_sgr_opt_kpp'); +add_default($nl, 'config_use_sgr_opt_temp_prescribed'); +add_default($nl, 'config_use_sgr_opt_salt_prescribed'); +add_default($nl, 'config_sgr_temperature_prescribed'); +add_default($nl, 'config_sgr_salinity_prescribed'); ############################ # Namelist group: coupling # diff --git a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml index c7441c6b2214..b70f7942686e 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -342,6 +342,14 @@ .true. 0.001 10.0 +'off' +0.001 +'bottom' +.true. +.false. +.false. +0.0 +0.0 .false. diff --git a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml index f6c2f9c5b765..601b121e5768 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -1188,6 +1188,70 @@ Valid values: Any positive real number. Default: Defined in namelist_defaults.xml + +Selects the mode in which subglacial runoff fluxes are implemented. + +Valid values: 'off', 'data' +Default: Defined in namelist_defaults.xml + + + +The length scale of exponential decay of subglacial runoff, used when config_sgr_flux_vertical_location is 'top' or 'bottom'. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$. + +Valid values: Any positive real number. +Default: Defined in namelist_defaults.xml + + + +Selects the vertical location where subglacial runoff is fluxed. + +Valid values: 'top','uniform', 'bottom' +Default: Defined in namelist_defaults.xml + + + +If true, include subglacial runoff (sgr) contribution in kpp calculation. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +If true, subglacial runoff temperature is set to config_sgr_temperature_prescribed. If false, the temperature is set to local freezing point. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +If true, subglacial runoff salinity is set to config_sgr_salinity_prescribed. If false, the salinity is set to 0 PSU. + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + + +Prescribed subglacial runoff temperature value, applied when config_use_sgr_opt_temp_prescribed = .true. + +Valid values: Any real number. +Default: Defined in namelist_defaults.xml + + + +Prescribed subglacial runoff salinity value, applied when config_use_sgr_opt_salt_prescribed = .true. + +Valid values: Any positive real number +Default: Defined in namelist_defaults.xml + + diff --git a/components/mpas-ocean/cime_config/buildnml b/components/mpas-ocean/cime_config/buildnml index 4c7ed2c0a311..8f5193c334be 100755 --- a/components/mpas-ocean/cime_config/buildnml +++ b/components/mpas-ocean/cime_config/buildnml @@ -35,6 +35,7 @@ def buildnml(case, caseroot, compname): ocn_forcing = case.get_value("MPASO_FORCING") ocn_iceberg = case.get_value("MPASO_ICEBERG") ocn_ismf = case.get_value("MPASO_ISMF") + ocn_sgr = case.get_value("MPASO_SGR") ocn_bgc = case.get_value("MPASO_BGC") ocn_wave = case.get_value("MPASO_WAVE") ocn_tidal_mixing = case.get_value("MPASO_TIDAL_MIXING") @@ -76,6 +77,7 @@ def buildnml(case, caseroot, compname): analysis_mask_file = '' eco_forcing_file = '' u_tidal_rms_file = '' + data_sgr_file = '' if ocn_grid == 'oEC60to30v3': decomp_date = '230424' @@ -140,6 +142,8 @@ def buildnml(case, caseroot, compname): data_ismf_file = 'prescribed_ismf_paolo2023.oQU240wLI.20240404.nc' if ocn_tidal_mixing == 'true': u_tidal_rms_file = 'velocityTidalRMS_CATS2008.oQU240wLI.20240221.nc' + if ocn_sgr == 'data': + data_sgr_file = 'DSGR.massFlux.MALI.out2055.oQU240wLI.20240328.nc' elif ocn_grid == 'oQU120': decomp_date = '230424' @@ -277,6 +281,8 @@ def buildnml(case, caseroot, compname): data_ismf_file = 'prescribed_ismf_adusumilli2020.SOwISC12to60E2r4.230516.nc' if ocn_tidal_mixing == 'true': u_tidal_rms_file = 'velocityTidalRMS_CATS2008.SOwISC12to60E2r4.20210114.nc' + if ocn_sgr == 'data': + data_sgr_file = 'DSGR.massFlux.MALI.out2055.SOwISC12to60E2r4.20240328.nc' elif ocn_grid == 'FRISwISC08to60E3r1': decomp_date = '20230913' # changed to date of partiotions in ../files_for_e3sm/assembled_files/inputdata/ocn/mpas-o/FRISwISC08to60E3r1/partitions @@ -352,6 +358,8 @@ def buildnml(case, caseroot, compname): data_ismf_file = 'prescribed_ismf_adusumilli2020.ECwISC30to60E2r1.230429.nc' if ocn_tidal_mixing == 'true': u_tidal_rms_file = 'velocityTidalRMS_CATS2008.ECwISC30to60E2r1.20240221.nc' + if ocn_sgr == 'data': + data_sgr_file = 'DSGR.massFlux.MALI.out2055.ECwISC30to60E2r1.20240328.nc' elif ocn_grid == 'IcoswISC30E3r5': decomp_date = '20231120' @@ -371,6 +379,8 @@ def buildnml(case, caseroot, compname): data_ismf_file = 'prescribed_ismf_paolo2023.IcoswISC30E3r5.20240227.nc' if ocn_tidal_mixing == 'true': u_tidal_rms_file = 'velocityTidalRMS_CATS2008.IcoswISC30E3r5.20231120.nc' + if ocn_sgr == 'data': + data_sgr_file = 'DSGR.massFlux.MALI.out2055.IcoswISC30E3r5.20240328.nc' elif ocn_grid == 'IcosXISC30E3r7': decomp_date = '20240314' @@ -396,7 +406,6 @@ def buildnml(case, caseroot, compname): if ocn_ismf == 'data': data_ismf_file = 'prescribed_ismf_paolo2023.RRSwISC6to18E3r5.20240327.nc' - #-------------------------------------------------------------------- # Set OCN_FORCING = datm_forced_restoring if restoring file is available #-------------------------------------------------------------------- @@ -436,6 +445,9 @@ def buildnml(case, caseroot, compname): if u_tidal_rms_file != '': input_list.write("u_tidal_rms = {}/ocn/mpas-o/{}/{}\n".format(din_loc_root, ocn_mask, u_tidal_rms_file)) + if data_sgr_file != '': + input_list.write("subglacial_runoff = {}/ocn/mpas-o/{}/{}\n".format(din_loc_root, ocn_mask, data_sgr_file)) + #-------------------------------------------------------------------- # Invoke mpas build-namelist - output will go in $CASEBUILD/mpasoconf #-------------------------------------------------------------------- @@ -495,6 +507,7 @@ def buildnml(case, caseroot, compname): sysmod += " -ocn_forcing '{}'".format(ocn_forcing) sysmod += " -ocn_iceberg '{}'".format(ocn_iceberg) sysmod += " -ocn_ismf '{}'".format(ocn_ismf) + sysmod += " -ocn_sgr '{}'".format(ocn_sgr) sysmod += " -ocn_bgc '{}'".format(ocn_bgc) sysmod += " -ocn_wave '{}'".format(ocn_wave) sysmod += " -ocn_tidal_mixing '{}'".format(ocn_tidal_mixing) @@ -702,6 +715,19 @@ def buildnml(case, caseroot, compname): lines.append('') lines.append('') + if data_sgr_file != '': + lines.append('') + lines.append('') + lines.append(' ') + lines.append('') + lines.append('') + if analysis_mask_file != '': lines.append('') lines.append('') lines.append('') + lines.append('') lines.append('') lines.append('') lines.append('') lines.append('') lines.append('') lines.append('') + lines.append('') lines.append('') lines.append('') lines.append('') @@ -1183,6 +1211,7 @@ def buildnml(case, caseroot, compname): lines.append('') lines.append('') lines.append('') + lines.append('') lines.append('') lines.append('') lines.append('') @@ -1331,6 +1360,7 @@ def buildnml(case, caseroot, compname): lines.append(' ') lines.append(' ') lines.append(' ') + lines.append(' ') lines.append(' ') lines.append(' ') lines.append(' ') @@ -1533,6 +1563,7 @@ def buildnml(case, caseroot, compname): lines.append(' ') lines.append(' ') + lines.append(' ') lines.append(' ') lines.append(' ') lines.append(' ') @@ -1598,6 +1629,7 @@ def buildnml(case, caseroot, compname): lines.append(' ') lines.append(' ') + lines.append(' ') lines.append(' ') lines.append(' ') lines.append(' ') @@ -1812,6 +1844,7 @@ def buildnml(case, caseroot, compname): lines.append(' ') lines.append(' ') lines.append(' ') + lines.append(' ') lines.append(' ') lines.append(' ') lines.append(' ') diff --git a/components/mpas-ocean/cime_config/config_component.xml b/components/mpas-ocean/cime_config/config_component.xml index e95b63682d95..83b4501ee5e5 100644 --- a/components/mpas-ocean/cime_config/config_component.xml +++ b/components/mpas-ocean/cime_config/config_component.xml @@ -39,6 +39,19 @@ Option to describe the MPASO prescribed tidal mixing + + char + none,data + none + + none + data + + case_comp + env_case.xml + Option to describe how MPASO will handle subglacial runoff fluxes + + char false,true diff --git a/components/mpas-ocean/cime_config/config_compsets.xml b/components/mpas-ocean/cime_config/config_compsets.xml index 448aa7a82b57..be36f0103472 100644 --- a/components/mpas-ocean/cime_config/config_compsets.xml +++ b/components/mpas-ocean/cime_config/config_compsets.xml @@ -37,6 +37,11 @@ 2000_DATM%NYF_SLND_MPASSI_MPASO%PISMFDATMFORCED_DROF%NYFAIS45_SGLC_SWAV + + GMPAS-NYF-PISMF-DSGR + 2000_DATM%NYF_SLND_MPASSI_MPASO%PISMFDATMFORCEDDSGR_DROF%NYFAIS45_SGLC_SWAV + + GMPAS-NYF-DISMF 2000_DATM%NYF_SLND_MPASSI_MPASO%DISMFDATMFORCED_DROF%NYFAIS45_SGLC_SWAV @@ -87,6 +92,16 @@ 2000_DATM%JRA-1p5_SLND_MPASSI%DIB_MPASO%IBPISMFDATMFORCED_DROF%JRA-1p5-AIS0ROF_SGLC_SWAV + + GMPAS-JRA1p5-DIB-PISMF-DSGR + 2000_DATM%JRA-1p5_SLND_MPASSI%DIB_MPASO%IBPISMFDATMFORCEDDSGR_DROF%JRA-1p5-AIS0ROF_SGLC_SWAV + + + + GMPAS-JRA1p5-DIB-PISMF-DSGR-TMIX + 2000_DATM%JRA-1p5_SLND_MPASSI%DIB_MPASO%IBPISMFDATMFORCEDDSGRTMIX_DROF%JRA-1p5-AIS0ROF_SGLC_SWAV + + GMPAS-JRA1p5-DIB-DISMF 2000_DATM%JRA-1p5_SLND_MPASSI%DIB_MPASO%IBDISMFDATMFORCED_DROF%JRA-1p5-AIS0ROF_SGLC_SWAV diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 825fe607a5b4..bf02edf96b51 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -702,6 +702,39 @@ description="The length scale of exponential decay of river runoff. Fluxes are multiplied by $e^{z/\gamma}$, where this coefficient is $\gamma$." possible_values="Any positive real number." /> + + + + + + + + + + @@ -2229,6 +2263,7 @@ + @@ -3347,6 +3382,10 @@ + + + @@ -3760,6 +3807,10 @@ description="Fresh water flux from river runoff at cell centers from coupler. Positive into the ocean." packages="thicknessBulkPKG" /> + + @@ -145,6 +149,10 @@ description="Fresh water flux from river runoff from coupler. Positive into the ocean." packages="thicknessBulkPKG" /> + @@ -194,6 +202,11 @@ + + + + @@ -302,6 +317,7 @@ + @@ -343,12 +359,14 @@ + + @@ -358,6 +376,7 @@ + diff --git a/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F b/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F index bdd07e799b5d..beed937d6767 100644 --- a/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F +++ b/components/mpas-ocean/src/analysis_members/mpas_ocn_conservation_check.F @@ -28,6 +28,7 @@ module ocn_conservation_check use ocn_constants use ocn_config + use ocn_mesh implicit none private @@ -427,6 +428,7 @@ subroutine energy_conservation(domain, err) accumulatedEvapTemperatureFlux, & accumulatedSeaIceTemperatureFlux, & accumulatedRiverRunoffTemperatureFlux, & + accumulatedSubglacialRunoffTemperatureFlux, & accumulatedIcebergTemperatureFlux real(kind=RKIND), dimension(:), allocatable :: & @@ -463,7 +465,8 @@ subroutine energy_conservation(domain, err) real(kind=RKIND), dimension(:,:), pointer :: & - activeTracersSurfaceFluxRunoff + activeTracersSurfaceFluxRunoff, & + activeTracersSurfaceFluxSubglacialRunoff type (MPAS_timeInterval_type) :: & timeStepESMF @@ -480,7 +483,7 @@ subroutine energy_conservation(domain, err) ierr integer, parameter :: & - nSums = 19 + nSums = 20 character(len=160) :: & m @@ -513,6 +516,7 @@ subroutine energy_conservation(domain, err) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedEvapTemperatureFlux", accumulatedEvapTemperatureFlux) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedSeaIceTemperatureFlux", accumulatedSeaIceTemperatureFlux) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedRiverRunoffTemperatureFlux", accumulatedRiverRunoffTemperatureFlux) + call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedSubglacialRunoffTemperatureFlux", accumulatedSubglacialRunoffTemperatureFlux) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedIcebergTemperatureFlux", accumulatedIcebergTemperatureFlux) !------------------------------------------------------------- @@ -556,6 +560,7 @@ subroutine energy_conservation(domain, err) call mpas_pool_get_subpool(forcingPool, 'tracersSurfaceFlux',tracersSurfaceFluxPool) call mpas_pool_get_dimension(tracersSurfaceFluxPool, 'index_temperatureSurfaceFlux', index_temperature_flux) call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFluxRunoff', activeTracersSurfaceFluxRunoff) + call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFluxSubglacialRunoff', activeTracersSurfaceFluxSubglacialRunoff) call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassNew, 2) call mpas_pool_get_array(statePool, 'accumulatedFrazilIceMass', accumulatedFrazilIceMassOld, 1) @@ -581,6 +586,11 @@ subroutine energy_conservation(domain, err) sumArray(14) = sumArray(14) + areaCell(iCell) * activeTracersSurfaceFluxRunoff(index_temperature_flux,iCell) sumArray(15) = sumArray(15) + areaCell(iCell) * icebergTemperatureFlux(iCell) + ! subglacial river runoff temperature flux + if (trim(config_subglacial_runoff_mode) == 'data') then + sumArray(20) = sumArray(20) + areaCell(iCell) * activeTracersSurfaceFluxSubglacialRunoff(index_temperature_flux,iCell) + end if + enddo ! iCell if (config_use_frazil_ice_formation) then @@ -634,6 +644,9 @@ subroutine energy_conservation(domain, err) accumulatedLandIceHeatFlux = accumulatedLandIceHeatFlux + sumArrayOut(17) accumulatedLandIceFrazilHeatFlux = accumulatedLandIceFrazilHeatFlux + sumArrayOut(18) accumulatedRemovedIceRunoffHeatFlux = accumulatedRemovedIceRunoffHeatFlux + sumArrayOut(19) + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffTemperatureFlux = accumulatedSubglacialRunoffTemperatureFlux + sumArrayOut(20) + end if ! cleanup deallocate(sumArray) @@ -663,6 +676,9 @@ subroutine energy_conservation(domain, err) accumulatedEvapTemperatureFlux = accumulatedEvapTemperatureFlux /accumulatedFluxCounter accumulatedSeaIceTemperatureFlux = accumulatedSeaIceTemperatureFlux /accumulatedFluxCounter accumulatedRiverRunoffTemperatureFlux = accumulatedRiverRunoffTemperatureFlux /accumulatedFluxCounter + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffTemperatureFlux = accumulatedSubglacialRunoffTemperatureFlux /accumulatedFluxCounter + end if accumulatedIcebergTemperatureFlux = accumulatedIcebergTemperatureFlux /accumulatedFluxCounter accumulatedLandIceFrazilHeatFlux = accumulatedLandIceFrazilHeatFlux /accumulatedFluxCounter accumulatedRemovedIceRunoffHeatFlux = accumulatedRemovedIceRunoffHeatFlux /accumulatedFluxCounter @@ -699,6 +715,9 @@ subroutine energy_conservation(domain, err) + accumulatedRiverRunoffTemperatureFlux *rho_sw*cp_sw & + accumulatedIcebergTemperatureFlux*rho_sw*cp_sw ! note, accumulatedLandIceFrazilHeatFlux not added because already in accumulatedFrazilHeatFlux + if (trim(config_subglacial_runoff_mode) == 'data') then + netEnergyFlux = netEnergyFlux + accumulatedSubglacialRunoffTemperatureFlux * rho_sw*cp_sw + end if ! compute the final energy error call MPAS_pool_get_array(conservationCheckEnergyAMPool, "absoluteEnergyError", absoluteEnergyError) @@ -760,6 +779,9 @@ subroutine energy_conservation(domain, err) v=accumulatedEvapTemperatureFlux *rho_sw*cp_sw; write(m,"('EvapTemperatureFlux ',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v v=accumulatedSeaIceTemperatureFlux *rho_sw*cp_sw; write(m,"('SeaIceTemperatureFlux ',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v v=accumulatedRiverRunoffTemperatureFlux*rho_sw*cp_sw; write(m,"('RiverRunoffTempFlux ',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v +if (trim(config_subglacial_runoff_mode) == 'data') then + v=accumulatedSubglacialRunoffTemperatureFlux*rho_sw*cp_sw; write(m,"('SubglacialRunoffTempFlux ',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v +end if v=accumulatedIcebergTemperatureFlux*rho_sw*cp_sw; write(m,"('IcebergTemperatureFlux ',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v write(m,"('SUM IMPLICIT HEAT FLUXES ',es16.8,' hh20temp ',f16.8,es16.8)") s, s/A; call mpas_log_write(m) s = s + explicitHeatFluxSum @@ -828,6 +850,7 @@ subroutine mass_conservation(domain, err) accumulatedEvaporationFlux, & accumulatedSeaIceFlux, & accumulatedRiverRunoffFlux, & + accumulatedSubglacialRunoffFlux, & accumulatedIceRunoffFlux, & accumulatedRemovedRiverRunoffFlux, & accumulatedRemovedIceRunoffFlux, & @@ -852,6 +875,7 @@ subroutine mass_conservation(domain, err) evaporationFlux, & seaIceFreshwaterFlux, & riverRunoffFlux, & + subglacialRunoffFlux, & iceRunoffFlux, & removedRiverRunoffFlux, & removedIceRunoffFlux, & @@ -877,7 +901,7 @@ subroutine mass_conservation(domain, err) iCell, ierr, k integer, parameter :: & - nSums = 12 + nSums = 13 integer, dimension(:), pointer :: minLevelCell, maxLevelCell @@ -900,6 +924,7 @@ subroutine mass_conservation(domain, err) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedEvaporationFlux", accumulatedEvaporationFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedSeaIceFlux", accumulatedSeaIceFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedRiverRunoffFlux", accumulatedRiverRunoffFlux) + call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedSubglacialRunoffFlux", accumulatedSubglacialRunoffFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedIceRunoffFlux", accumulatedIceRunoffFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedRemovedRiverRunoffFlux",accumulatedRemovedRiverRunoffFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedRemovedIceRunoffFlux", accumulatedRemovedIceRunoffFlux) @@ -937,6 +962,7 @@ subroutine mass_conservation(domain, err) call mpas_pool_get_array(forcingPool, 'evaporationFlux', evaporationFlux) call mpas_pool_get_array(forcingPool, 'seaIceFreshWaterFlux', seaIceFreshwaterFlux) call mpas_pool_get_array(forcingPool, 'riverRunoffFlux', riverRunoffFlux) + call mpas_pool_get_array(forcingPool, 'subglacialRunoffFlux', subglacialRunoffFlux) call mpas_pool_get_array(forcingPool, 'iceRunoffFlux', iceRunoffFlux) call mpas_pool_get_array(forcingPool, 'removedRiverRunoffFlux', removedRiverRunoffFlux) call mpas_pool_get_array(forcingPool, 'removedIceRunoffFlux', removedIceRunoffFlux) @@ -956,6 +982,9 @@ subroutine mass_conservation(domain, err) sumArray( 7) = sumArray( 7) + areaCell(iCell) * removedRiverRunoffFlux(iCell) sumArray( 8) = sumArray( 8) + areaCell(iCell) * removedIceRunoffFlux(iCell) sumArray( 9) = sumArray( 9) + areaCell(iCell) * icebergFreshwaterFlux(iCell) + if (trim(config_subglacial_runoff_mode) == 'data') then + sumArray(13) = sumArray(13) + areaCell(iCell) * subglacialRunoffFlux(iCell) + end if enddo if (config_use_frazil_ice_formation) then @@ -1002,6 +1031,9 @@ subroutine mass_conservation(domain, err) accumulatedFrazilFlux = accumulatedFrazilFlux + sumArrayOut(10) accumulatedLandIceFlux = accumulatedLandIceFlux + sumArrayOut(11) accumulatedLandIceFrazilFlux = accumulatedLandIceFrazilFlux + sumArrayOut(12) + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffFlux = accumulatedSubglacialRunoffFlux + sumArrayOut(13) + end if ! cleanup deallocate(sumArray) @@ -1021,6 +1053,9 @@ subroutine mass_conservation(domain, err) accumulatedEvaporationFlux = accumulatedEvaporationFlux /accumulatedFluxCounter accumulatedSeaIceFlux = accumulatedSeaIceFlux /accumulatedFluxCounter accumulatedRiverRunoffFlux = accumulatedRiverRunoffFlux /accumulatedFluxCounter + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffFlux = accumulatedSubglacialRunoffFlux /accumulatedFluxCounter + end if accumulatedIceRunoffFlux = accumulatedIceRunoffFlux /accumulatedFluxCounter accumulatedRemovedRiverRunoffFlux = accumulatedRemovedRiverRunoffFlux /accumulatedFluxCounter accumulatedRemovedIceRunoffFlux = accumulatedRemovedIceRunoffFlux /accumulatedFluxCounter @@ -1054,6 +1089,9 @@ subroutine mass_conservation(domain, err) + accumulatedFrazilFlux & + accumulatedLandIceFlux ! note, accumulatedLandIceFrazilFlux not added because already in accumulatedFrazilFlux + if (trim(config_subglacial_runoff_mode) == 'data') then + netMassFlux = netMassFlux + accumulatedSubglacialRunoffFlux + end if ! compute the final mass error call MPAS_pool_get_array(conservationCheckMassAMPool, "absoluteMassError", absoluteMassError) @@ -1083,6 +1121,10 @@ subroutine mass_conservation(domain, err) v=accumulatedIcebergFlux ; write(m,"('icebergFreshwaterFlux ',es16.8,' x2o_Fioi_bergw wberg ',f16.8)") v,v*c; call mpas_log_write(m); s=s+v v=accumulatedEvaporationFlux ; write(m,"('evaporationFlux ',es16.8,' x2o_Foxx_evap wevap ',f16.8)") v,v*c; call mpas_log_write(m); s=s+v v=accumulatedRiverRunoffFlux ; write(m,"('riverRunoffFlux ',es16.8,' x2o_Foxx_rofl wrunoff ',f16.8)") v,v*c; call mpas_log_write(m); s=s+v +if (trim(config_subglacial_runoff_mode) == 'data') then + v=accumulatedSubglacialRunoffFlux ; write(m,"('subglacialRunoffFlux ',es16.8,' x2o_Foxx_rofl wsgr ',f16.8)") v,v*c; call mpas_log_write(m); s=s+v +end if + if (landIceFreshwaterFluxesOn) then v=accumulatedRemovedRiverRunoffFlux; write(m,"('removedRiverRunoffFlux ',es16.8,' wrunoff ',f16.8)") v,v*c; call mpas_log_write(m); v=accumulatedRiverRunoffFlux+accumulatedRemovedRiverRunoffFlux; @@ -1150,7 +1192,8 @@ subroutine salt_conservation(domain, err) conservationCheckSaltAMPool, & meshPool, & statePool, & - forcingPool + forcingPool, & + tracersSurfaceFluxPool real(kind=RKIND), pointer :: & initialSalt, & @@ -1160,9 +1203,13 @@ subroutine salt_conservation(domain, err) absoluteSaltError, & relativeSaltError + real(kind=RKIND), dimension(:,:), pointer :: & + activeTracersSurfaceFluxSubglacialRunoff + real(kind=RKIND), pointer :: & accumulatedSeaIceSalinityFlux, & - accumulatedFrazilSalinityFlux + accumulatedFrazilSalinityFlux, & + accumulatedSubglacialRunoffSalinityFlux ! accumulatedLandIceFrazilSalinityFlux is not present because it is always 0 real(kind=RKIND), dimension(:), allocatable :: & @@ -1184,14 +1231,15 @@ subroutine salt_conservation(domain, err) dt, dtAvg, v, A, s, c integer, pointer :: & - nCellsSolve + nCellsSolve, & + index_salinity_flux integer :: & iCell, & ierr integer, parameter :: & - nSums = 3 + nSums = 4 logical, pointer :: & activeTracersBulkRestoringPKG @@ -1209,6 +1257,7 @@ subroutine salt_conservation(domain, err) call MPAS_pool_get_array(conservationCheckSaltAMPool, "accumulatedSeaIceSalinityFlux", accumulatedSeaIceSalinityFlux) call MPAS_pool_get_array(conservationCheckSaltAMPool, "accumulatedFrazilSalinityFlux", accumulatedFrazilSalinityFlux) + call MPAS_pool_get_array(conservationCheckSaltAMPool, "accumulatedSubglacialRunoffSalinityFlux", accumulatedSubglacialRunoffSalinityFlux) !------------------------------------------------------------- ! Net salt flux to ice @@ -1234,12 +1283,23 @@ subroutine salt_conservation(domain, err) call mpas_pool_get_array(forcingPool, 'seaIceSalinityFlux', seaIceSalinityFlux) call mpas_pool_get_array(statePool, 'accumulatedFrazilIceSalinity', accumulatedFrazilIceSalinityOld, 1) call mpas_pool_get_array(statePool, 'accumulatedFrazilIceSalinity', accumulatedFrazilIceSalinityNew, 2) + call mpas_pool_get_subpool(forcingPool, 'tracersSurfaceFlux',tracersSurfaceFluxPool) + call mpas_pool_get_array(tracersSurfaceFluxPool, 'activeTracersSurfaceFluxSubglacialRunoff', activeTracersSurfaceFluxSubglacialRunoff) + call mpas_pool_get_dimension(tracersSurfaceFluxPool, 'index_salinitySurfaceFlux', index_salinity_flux) + + do iCell = 1, nCellsSolve ! salt flux to ocean sumArray(1) = sumArray(1) + areaCell(iCell) * seaIceSalinityFlux(iCell) enddo ! iCell + ! subglacial runoff + if (trim(config_subglacial_runoff_mode) == 'data') then + do iCell = 1, nCellsSolve + sumArray(4) = sumArray(4) + areaCell(iCell) * activeTracersSurfaceFluxSubglacialRunoff(index_salinity_flux,iCell) + enddo ! iCell + end if if (config_use_frazil_ice_formation) then do iCell = 1, nCellsSolve @@ -1265,6 +1325,9 @@ subroutine salt_conservation(domain, err) ! accumulate fluxes accumulatedSeaIceSalinityFlux = accumulatedSeaIceSalinityFlux + sumArrayOut(1) accumulatedFrazilSalinityFlux = accumulatedFrazilSalinityFlux + sumArrayOut(2) + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffSalinityFlux = accumulatedSubglacialRunoffSalinityFlux + sumArrayOut(4) + end if ! cleanup deallocate(sumArray) @@ -1281,6 +1344,9 @@ subroutine salt_conservation(domain, err) ! Average the fluxes accumulatedSeaIceSalinityFlux = accumulatedSeaIceSalinityFlux /accumulatedFluxCounter accumulatedFrazilSalinityFlux = accumulatedFrazilSalinityFlux /accumulatedFluxCounter + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffSalinityFlux = accumulatedSubglacialRunoffSalinityFlux /accumulatedFluxCounter + end if ! get initial salt content call MPAS_pool_get_array(conservationCheckSaltAMPool, "initialSalt", initialSalt) @@ -1299,6 +1365,10 @@ subroutine salt_conservation(domain, err) netSaltFlux = accumulatedSeaIceSalinityFlux & + accumulatedFrazilSalinityFlux + if (trim(config_subglacial_runoff_mode) == 'data') then + netSaltFlux = netSaltFlux + accumulatedSubglacialRunoffSalinityFlux + end if + ! compute the final salt error call MPAS_pool_get_array(conservationCheckSaltAMPool, "absoluteSaltError", absoluteSaltError) call MPAS_pool_get_array(conservationCheckSaltAMPool, "relativeSaltError", relativeSaltError) @@ -1327,6 +1397,10 @@ subroutine salt_conservation(domain, err) .and.config_frazil_under_land_ice) then v=0; write(m,"('LandIceFrazilSalinityFlux',es16.8,' (already in wmelt, do not sum) ',f16.8)") v,v*c; call mpas_log_write(m); !no sum: s=s+v end if +if (trim(config_subglacial_runoff_mode) == 'data') then + v=accumulatedSubglacialRunoffSalinityFlux ; write(m,"('SubglacialRunoffSalinityFlux ',es16.8,' x2o_Fioi_salt salt ',f16.8)") v,v*c; call mpas_log_write(m); s=s+v +end if + write(m,"('SUM VOLUME FLUXES ',es16.8,' ',f16.8,es16.8)") s, s*c; call mpas_log_write(m) call mpas_log_write(' ') @@ -1523,11 +1597,11 @@ subroutine carbon_conservation(domain, err) block => domain % blocklist do while (associated(block)) call MPAS_pool_get_dimension(block % dimensions, "nCellsSolve", nCellsSolve) - + call MPAS_pool_get_subpool(block % structs, "mesh", meshPool) call MPAS_pool_get_subpool(block % structs, "forcing", forcingPool) call MPAS_pool_get_subpool(block % structs, "state", statePool) - + call MPAS_pool_get_array(meshPool, "areaCell", areaCell) call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) @@ -1594,7 +1668,7 @@ subroutine carbon_conservation(domain, err) + ecosysTracersSurfaceFlux(index_spCFlux, iCell) & + ecosysTracersSurfaceFlux(index_diatCFlux, iCell) & + ecosysTracersSurfaceFlux(index_diazCFlux, iCell)) - + do k = 1, maxLevelCell(iCell) sumArray(4) = sumArray(4) + areaCell(iCell) * ( & ecosysTracersTend(index_DICTend, k, iCell) & @@ -1645,7 +1719,7 @@ subroutine carbon_conservation(domain, err) !------------------------------------------------------------- if (MPAS_stream_mgr_ringing_alarms(domain % streamManager, "conservationCheckOutput", ierr=ierr)) then - + ! Average the fluxes accumulatedCarbonSourceSink = accumulatedCarbonSourceSink /accumulatedFluxCounter accumulatedCarbonSedimentFlux = accumulatedCarbonSedimentFlux /accumulatedFluxCounter @@ -1657,7 +1731,7 @@ subroutine carbon_conservation(domain, err) ! get initial carbon call MPAS_pool_get_array(conservationCheckCarbonAMPool, "initialCarbon", initialCarbon) - + ! get final carbon call MPAS_pool_get_array(conservationCheckCarbonAMPool, "finalCarbon", finalCarbon) call compute_total_carbon(domain, finalCarbon) @@ -1694,13 +1768,13 @@ subroutine carbon_conservation(domain, err) relativeCarbonErrorBounds = & relativeCarbonErrorStepBounds*relativeCarbonErrorBoundsFac*totalTimeSteps relativeCarbonErrorPerTimeStep = relativeCarbonError/accumulatedFluxCounter - + !------------------------------------------------------------- ! Output to log file !------------------------------------------------------------- if (config_AM_conservationCheck_write_to_logfile) then - + call mpas_log_write('') call mpas_log_write('----------------------------------------------------------') call mpas_log_write('CARBON CONSERVATION CHECK') @@ -1731,7 +1805,7 @@ subroutine carbon_conservation(domain, err) v=accumulatedIceOceanInorganicCarbonFlux*mmol_to_kg_C write(m,"('Ice-Ocean Inorganic Flux ',es16.8,' ',f16.8)") v,v*c call mpas_log_write(m) - write(m,"('SUM FLUXES (surf + sed) ',es16.8,' ',f16.8,es16.8)") s, s*c; + write(m,"('SUM FLUXES (surf + sed) ',es16.8,' ',f16.8,es16.8)") s, s*c; call mpas_log_write(m) call mpas_log_write(' ') @@ -1751,24 +1825,24 @@ subroutine carbon_conservation(domain, err) write(m,"('Carbon change ', 3es16.8)") & carbonChange*mmol_to_kg_C, & carbonChange*mmol_to_kg_C/dtAvg, & - carbonChange*mmol_to_kg_C/dtAvg*c + carbonChange*mmol_to_kg_C/dtAvg*c call mpas_log_write(m) write(m,"('Net carbon flux ', 3es16.8)") & netCarbonFlux*mmol_to_kg_C*dtAvg, & netCarbonFlux*mmol_to_kg_C, & - netCarbonFlux*mmol_to_kg_C*c + netCarbonFlux*mmol_to_kg_C*c call mpas_log_write(m) write(m,"('Absolute carbon error ', 3es16.8)") & absoluteCarbonError*mmol_to_kg_C, & absoluteCarbonError*mmol_to_kg_C/dtAvg, & - absoluteCarbonError*mmol_to_kg_C/dtAvg*c + absoluteCarbonError*mmol_to_kg_C/dtAvg*c call mpas_log_write(m) call mpas_log_write(' ') write(m,"('RELATIVE CARBON ERROR =', es16.8)") & relativeCarbonError call mpas_log_write(m) call mpas_log_write(' ') - + write(m,"('Relative carbon error per timestep = ', es16.8)") & relativeCarbonErrorPerTimeStep call mpas_log_write(m) @@ -1814,7 +1888,7 @@ subroutine carbon_conservation(domain, err) endif end subroutine carbon_conservation - + !*********************************************************************** ! ! routine compute_total_energy @@ -2190,7 +2264,7 @@ subroutine compute_total_carbon(domain, totalCarbon) call MPAS_dmpar_sum_real(domain % dminfo, carbon, totalCarbon) end subroutine compute_total_carbon - + !*********************************************************************** ! ! routine reset_accumulated_variables @@ -2232,6 +2306,7 @@ subroutine reset_accumulated_variables(domain) accumulatedEvapTemperatureFlux, & accumulatedSeaIceTemperatureFlux, & accumulatedRiverRunoffTemperatureFlux, & + accumulatedSubglacialRunoffTemperatureFlux, & accumulatedIcebergTemperatureFlux real(kind=RKIND), pointer :: & @@ -2240,6 +2315,7 @@ subroutine reset_accumulated_variables(domain) accumulatedEvaporationFlux, & accumulatedSeaIceFlux, & accumulatedRiverRunoffFlux, & + accumulatedSubglacialRunoffFlux, & accumulatedIceRunoffFlux, & accumulatedRemovedRiverRunoffFlux, & accumulatedRemovedIceRunoffFlux, & @@ -2250,7 +2326,8 @@ subroutine reset_accumulated_variables(domain) real(kind=RKIND), pointer :: & accumulatedSeaIceSalinityFlux, & - accumulatedFrazilSalinityFlux + accumulatedFrazilSalinityFlux, & + accumulatedSubglacialRunoffSalinityFlux real(kind=RKIND), pointer :: & accumulatedCarbonSourceSink, & @@ -2280,6 +2357,7 @@ subroutine reset_accumulated_variables(domain) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedEvapTemperatureFlux", accumulatedEvapTemperatureFlux) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedSeaIceTemperatureFlux", accumulatedSeaIceTemperatureFlux) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedRiverRunoffTemperatureFlux", accumulatedRiverRunoffTemperatureFlux) + call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedSubglacialRunoffTemperatureFlux", accumulatedSubglacialRunoffTemperatureFlux) call MPAS_pool_get_array(conservationCheckEnergyAMPool, "accumulatedIcebergTemperatureFlux", accumulatedIcebergTemperatureFlux) accumulatedFluxCounter = 0 @@ -2298,6 +2376,9 @@ subroutine reset_accumulated_variables(domain) accumulatedEvapTemperatureFlux = 0.0_RKIND accumulatedSeaIceTemperatureFlux = 0.0_RKIND accumulatedRiverRunoffTemperatureFlux = 0.0_RKIND + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffTemperatureFlux = 0.0_RKIND + end if accumulatedIcebergTemperatureFlux = 0.0_RKIND accumulatedLandIceFrazilHeatFlux = 0.0_RKIND accumulatedRemovedIceRunoffHeatFlux = 0.0_RKIND @@ -2310,6 +2391,7 @@ subroutine reset_accumulated_variables(domain) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedEvaporationFlux", accumulatedEvaporationFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedSeaIceFlux", accumulatedSeaIceFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedRiverRunoffFlux", accumulatedRiverRunoffFlux) + call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedSubglacialRunoffFlux", accumulatedSubglacialRunoffFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedIceRunoffFlux", accumulatedIceRunoffFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedRemovedRiverRunoffFlux",accumulatedRemovedRiverRunoffFlux) call MPAS_pool_get_array(conservationCheckMassAMPool, "accumulatedRemovedIceRunoffFlux", accumulatedRemovedIceRunoffFlux) @@ -2323,6 +2405,9 @@ subroutine reset_accumulated_variables(domain) accumulatedEvaporationFlux = 0.0_RKIND accumulatedSeaIceFlux = 0.0_RKIND accumulatedRiverRunoffFlux = 0.0_RKIND + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffFlux = 0.0_RKIND + end if accumulatedIceRunoffFlux = 0.0_RKIND accumulatedRemovedRiverRunoffFlux = 0.0_RKIND accumulatedRemovedIceRunoffFlux = 0.0_RKIND @@ -2336,9 +2421,13 @@ subroutine reset_accumulated_variables(domain) call MPAS_pool_get_array(conservationCheckSaltAMPool, "accumulatedSeaIceSalinityFlux", accumulatedSeaIceSalinityFlux) call MPAS_pool_get_array(conservationCheckSaltAMPool, "accumulatedFrazilSalinityFlux", accumulatedFrazilSalinityFlux) + call MPAS_pool_get_array(conservationCheckSaltAMPool, "accumulatedSubglacialRunoffSalinityFlux", accumulatedSubglacialRunoffSalinityFlux) accumulatedSeaIceSalinityFlux = 0.0_RKIND accumulatedFrazilSalinityFlux = 0.0_RKIND + if (trim(config_subglacial_runoff_mode) == 'data') then + accumulatedSubglacialRunoffSalinityFlux = 0.0_RKIND + end if call MPAS_pool_get_subpool(domain % blocklist % structs, "conservationCheckCarbonAM", conservationCheckCarbonAMPool) diff --git a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F index 7d9eb083f907..6d542fb0c3f3 100644 --- a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F +++ b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F @@ -123,6 +123,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: dataLandIceFluxesPKGActive logical, pointer :: landIceFluxesPKGActive logical, pointer :: landIceCouplingPKGActive + logical, pointer :: dataSubglacialRunoffFluxPKGActive logical, pointer :: thicknessBulkPKGActive logical, pointer :: frazilIceActive logical, pointer :: tidalForcingActive @@ -186,6 +187,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: config_use_bulk_thickness_flux logical, pointer :: config_compute_active_tracer_budgets character (len=StrKIND), pointer :: config_land_ice_flux_mode + character (len=StrKIND), pointer :: config_subglacial_runoff_mode type (mpas_pool_iterator_type) :: groupItr character (len=StrKIND) :: tracerGroupName, configName, packageName @@ -319,6 +321,15 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ landIceCouplingPKGActive = .true. end if + ! + ! test for use of subglacial runoff flux, dataSubglacialRunoffFluxPKGActive + ! + call mpas_pool_get_package(packagePool, 'dataSubglacialRunoffFluxPKGActive', dataSubglacialRunoffFluxPKGActive) + call mpas_pool_get_config(configPool, 'config_subglacial_runoff_mode', config_subglacial_runoff_mode) + if ( trim(config_subglacial_runoff_mode) == 'data' ) then + dataSubglacialRunoffFluxPKGActive = .true. + end if + ! ! test for use of frazil ice formation, frazilIceActive ! diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index 52ef08267fc1..a2d6586a7d63 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -218,6 +218,7 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe !! tracersSurfaceValue calc :: minLevelCell (mesh array) !! normalVelocitySurfaceLayer calc :: minLevelEdgeBot (mesh array) !! surfaceFluxAttenuationCoefficientRunoff :: layerThickEdgeFlux (diagnostics array) + !! surfaceFluxAttenuationCoefficientSubglacialRunoff :: layerThickEdgeFlux (diagnostics array) !! ocn_diagnostic_solve_ssh :: ssh, (not on device) !! landIceDraft (not on device) @@ -457,6 +458,11 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe indexSurfaceLayerDepth, normalVelocitySurfaceLayer, & sfcFlxAttCoeff, surfaceFluxAttenuationCoefficientRunoff) + if (trim(config_subglacial_runoff_mode) == 'data') then + ! outputs: surfaceFluxAttenuationCoefficientSubglacialRunoff + call ocn_diagnostic_solve_surfaceLayer_subglacialRunoff(surfaceFluxAttenuationCoefficientSubglacialRunoff) + end if + end if ! full_compute ! @@ -522,6 +528,11 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe !! normalVelocitySurfaceLayer, (diagnostics array) !! surfaceFluxAttenuationCoefficient, ???? !! surfaceFluxAttenuationCoefficientRunoff (diagnostics array) + !! surfaceFluxAttenuationCoefficientSubglacialRunoff :: tracersSurfaceLayerValue, (diagnostics array) + !! indexSurfaceLayerDepth, (diagnostics array) + !! normalVelocitySurfaceLayer, (diagnostics array) + !! surfaceFluxAttenuationCoefficient, ???? + !! surfaceFluxAttenuationCoefficientSubglacialRunoff (diagnostics array) !! ocn_diagnostic_solve_ssh :: pressureAdjustedSSH, (diagnostics array) !! gradSSH (diagnostics array) @@ -558,6 +569,10 @@ subroutine ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMe ! !$acc normalVelocitySurfaceLayer, & ! !$acc sfcFlxAttCoeff, & ! !$acc surfaceFluxAttenuationCoefficientRunoff) +! if (trim(config_subglacial_runoff_mode) == 'data') then +! !$acc update host(surfaceFluxAttenuationCoefficientSubglacialRunoff) +! end if + ! end if ! full_compute ! !$acc exit data delete (atmosphericPressure, seaIcePressure) @@ -1902,6 +1917,69 @@ subroutine ocn_diagnostic_solve_surfaceLayer(layerThickness, & end subroutine ocn_diagnostic_solve_surfaceLayer !}}} +!*********************************************************************** +! +! routine ocn_diagnostic_solve_surfaceLayer_subglacialRunoff +! +!> \brief Computes diagnostic subglacialRunoff variables for surface layer +!> \author Irena Vankova +!> \date July 2024 +!> \details +! +!----------------------------------------------------------------------- + + subroutine ocn_diagnostic_solve_surfaceLayer_subglacialRunoff(surfaceFluxAttenuationCoefficientSubglacialRunoff)!{{{ + + !----------------------------------------------------------------- + ! input variables + !----------------------------------------------------------------- + + !----------------------------------------------------------------- + ! output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:), intent(out) :: & + surfaceFluxAttenuationCoefficientSubglacialRunoff !< [out] sfc flx attenuation + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + nCells, &! num of cells + iCell ! loop indices for cell + + ! End preamble + !----------------------------------------------------------------- + ! Begin code + + ! compute the attenuation coefficient for subglacial runoff surface fluxes + + nCells = nCellsHalo( 1 ) + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(surfaceFluxAttenuationCoefficientSubglacialRunoff) +#else + !$omp parallel + !$omp do schedule(runtime) +#endif + do iCell = 1, nCells + surfaceFluxAttenuationCoefficientSubglacialRunoff(iCell) = & + config_flux_attenuation_coefficient_subglacial_runoff + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + !-------------------------------------------------------------------- + + end subroutine ocn_diagnostic_solve_surfaceLayer_subglacialRunoff !}}} + +!*********************************************************************** + + !*********************************************************************** ! ! routine ocn_diagnostic_solve_GMvel @@ -3275,15 +3353,18 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & !----------------------------------------------------------------- integer :: & - iCell, iEdge, i, &! loop indices for cell, edge, neighbors - kmin, &! topmost active cell index - nCells, &! number of cells - err, &! local error code - timeLevel ! time level for state variables (default 1) + iCell, iEdge, i, k, &! loop indices for cell, edge, neighbors + kmin, &! topmost active cell index + nCells, &! number of cells + err, &! local error code + timeLevel ! time level for state variables (default 1) real (kind=RKIND) :: & fracAbsorbed, &! fraction of sfc flux absorbed fracAbsorbedRunoff, &! same for runoff + fracAbsorbedSubglacialRunoff, &! same for subglacial runoff + zTop,zBot, &! temporary variables + transmissionCoeffTop,transmissionCoeffBot, &! temporary variables sumSurfaceStressSquared ! sum of sfc stress squared ! pointers for variable/pool retrievals @@ -3299,6 +3380,7 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & penetrativeTemperatureFlux, &! various sfc flux components surfaceThicknessFlux, & surfaceThicknessFluxRunoff, & + surfaceThicknessFluxSubglacialRunoff, & rainTemperatureFlux, & evapTemperatureFlux, & icebergTemperatureFlux, & @@ -3311,6 +3393,7 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & normalVelocity, &! normal velocity activeTracersSurfaceFlux, &! sfc flux of active tracers (T,S) activeTracersSurfaceFluxRunoff, &! flx of tracers in runoff + activeTracersSurfaceFluxSubglacialRunoff, &! flx of tracers in subglacial runoff nonLocalSurfaceTracerFlux ! non-local flux of tracers real (kind=RKIND), dimension(:,:,:), pointer :: & @@ -3374,11 +3457,16 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & call mpas_pool_get_array(tracersSurfaceFluxPool, & 'activeTracersSurfaceFluxRunoff', & activeTracersSurfaceFluxRunoff) + call mpas_pool_get_array(tracersSurfaceFluxPool, & + 'activeTracersSurfaceFluxSubglacialRunoff', & + activeTracersSurfaceFluxSubglacialRunoff) call mpas_pool_get_array(forcingPool, 'surfaceThicknessFlux', & surfaceThicknessFlux) call mpas_pool_get_array(forcingPool, 'surfaceThicknessFluxRunoff', & surfaceThicknessFluxRunoff) + call mpas_pool_get_array(forcingPool, 'surfaceThicknessFluxSubglacialRunoff', & + surfaceThicknessFluxSubglacialRunoff) call mpas_pool_get_array(forcingPool, 'penetrativeTemperatureFlux', & penetrativeTemperatureFlux) call mpas_pool_get_array(forcingPool, 'surfaceStress', & @@ -3451,7 +3539,7 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & !$omp parallel !$omp do schedule(runtime) & - !$omp private(kmin, fracAbsorbed, fracAbsorbedRunoff, & + !$omp private(kmin, fracAbsorbed, fracAbsorbedRunoff, fracAbsorbedSubglacialRunoff, & !$omp sumSurfaceStressSquared, i, iEdge) do iCell = 1, nCells @@ -3465,6 +3553,37 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & - exp( max(-100.0_RKIND, & -layerThickness(kmin, iCell)/ & config_flux_attenuation_coefficient_runoff)) + if (trim(config_subglacial_runoff_mode) == 'data') then + if (config_use_sgr_opt_kpp) then + if ( trim(config_sgr_flux_vertical_location) == 'top' ) then + fracAbsorbedSubglacialRunoff = 1.0_RKIND & + - exp( max(-100.0_RKIND, & + -layerThickness(kmin, iCell)/ & + config_flux_attenuation_coefficient_subglacial_runoff)) + else if ( trim(config_sgr_flux_vertical_location) == 'uniform' ) then + ! calculate total thickness into variable zTop + zTop = 0.0_RKIND + do k = minLevelCell(iCell), maxLevelCell(iCell) + zTop = zTop + layerThickness(k,iCell) + end do + ! distribute flux evenly throughout water column + fracAbsorbedSubglacialRunoff = layerThickness(kmin, iCell) / zTop + else if ( trim(config_sgr_flux_vertical_location) == 'bottom' ) then + zTop = 0.0_RKIND + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + zBot = zTop - layerThickness(k,iCell) + if (k == minLevelCell(iCell)) then + transmissionCoeffTop = exp( max(zTop / config_flux_attenuation_coefficient_runoff, -100.0_RKIND) ) + transmissionCoeffBot = exp( max(zBot / config_flux_attenuation_coefficient_runoff, -100.0_RKIND) ) + fracAbsorbedSubglacialRunoff = transmissionCoeffTop - transmissionCoeffBot + end if + zTop = zBot + end do + end if + else + fracAbsorbedSubglacialRunoff = 0.0_RKIND + end if + end if ! Store the total tracer flux below in ! nonLocalSurfaceTemperatureFlux for use in the CVMix nonlocal @@ -3480,6 +3599,11 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & icebergTemperatureFlux(iCell)) & - fracAbsorbedRunoff* & activeTracersSurfaceFluxRunoff(indexTempFlux,iCell) + if (trim(config_subglacial_runoff_mode) == 'data') then + nonLocalSurfaceTracerFlux(indexTempFlux, iCell) = nonLocalSurfaceTracerFlux(indexTempFlux, iCell) & + - fracAbsorbedSubglacialRunoff* & + activeTracersSurfaceFluxSubglacialRunoff(indexTempFlux,iCell) + end if nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) = & activeTracersSurfaceFlux(indexSaltFlux,iCell) & @@ -3487,6 +3611,11 @@ subroutine ocn_compute_KPP_input_fields(statePool, forcingPool, & activeTracers(indexSaltFlux,kmin,iCell) & - fracAbsorbedRunoff*surfaceThicknessFluxRunoff(iCell)* & activeTracers(indexSaltFlux,kmin,iCell) + if (trim(config_subglacial_runoff_mode) == 'data') then + nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) = nonLocalSurfaceTracerFlux(indexSaltFlux,iCell) & + - fracAbsorbedSubglacialRunoff*surfaceThicknessFluxSubglacialRunoff(iCell)* & + activeTracers(indexSaltFlux,kmin,iCell) + end if surfaceBuoyancyForcing(iCell) = & thermalExpansionCoeff(kmin,iCell)* & @@ -4428,6 +4557,21 @@ subroutine ocn_validate_state(domain, timeLevel)!{{{ call ocn_write_field_statistics(debugUnit, fieldName, minValue, maxValue) end if + if (trim(config_subglacial_runoff_mode) == 'data') then + ! Test subglacialRunoffFlux + fieldName = 'subglacialRunoffFlux' + minValue = HUGE(minValue) + maxValue = -HUGE(maxValue) + call mpas_pool_get_array(forcingPool, fieldName, real1DArr) + if ( associated(real1DArr) ) then + do iCell = 1, nCellsSolve + minValue = min( minValue, real1DArr(iCell) ) + maxValue = max( maxValue, real1DArr(iCell) ) + end do + call ocn_write_field_statistics(debugUnit, fieldName, minValue, maxValue) + end if + end if + ! Test seaIceSalinityFlux fieldName = 'seaIceSalinityFlux' minValue = HUGE(minValue) @@ -4598,6 +4742,12 @@ subroutine ocn_diagnostics_init(domain, err)!{{{ call ocn_diagnostics_variables_init(domain, jenkinsOn, & hollandJenkinsOn, err) + if ( trim(config_sgr_flux_vertical_location) /= 'top' .and. & + trim(config_sgr_flux_vertical_location) /= 'uniform' .and. & + trim(config_sgr_flux_vertical_location) /= 'bottom' ) then + call mpas_log_write("config_sgr_flux_vertical_location not one of 'bottom', 'uniform', 'top'.", MPAS_LOG_CRIT) + end if + end subroutine ocn_diagnostics_init!}}} !*********************************************************************** diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F index 30a62f68c8b9..38f32f88e17c 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -85,6 +85,7 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:), pointer :: sfcFlxAttCoeff real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientRunoff + real (kind=RKIND), dimension(:), pointer :: surfaceFluxAttenuationCoefficientSubglacialRunoff real (kind=RKIND), dimension(:), pointer :: landIceFrictionVelocity real (kind=RKIND), dimension(:), pointer :: velocityTidalRMS @@ -680,6 +681,9 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e call mpas_pool_get_array(diagnosticsPool, & 'surfaceFluxAttenuationCoefficientRunoff', & surfaceFluxAttenuationCoefficientRunoff) + call mpas_pool_get_array(diagnosticsPool, & + 'surfaceFluxAttenuationCoefficientSubglacialRunoff', & + surfaceFluxAttenuationCoefficientSubglacialRunoff) call mpas_pool_get_array(diagnosticsPool, & 'boundaryLayerDepth', & boundaryLayerDepth) @@ -964,7 +968,8 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc Time_bnds, & !$acc simulationStartTime, & !$acc boundaryLayerDepthSmooth, & - !$acc surfaceFluxAttenuationCoefficientRunoff & + !$acc surfaceFluxAttenuationCoefficientRunoff, & + !$acc surfaceFluxAttenuationCoefficientSubglacialRunoff & !$acc ) end subroutine ocn_diagnostics_variables_init!}}} @@ -1223,7 +1228,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc Time_bnds, & !$acc simulationStartTime, & !$acc boundaryLayerDepthSmooth, & - !$acc surfaceFluxAttenuationCoefficientRunoff & + !$acc surfaceFluxAttenuationCoefficientRunoff, & + !$acc surfaceFluxAttenuationCoefficientSubglacialRunoff & !$acc ) ! Nullify pointers @@ -1425,7 +1431,8 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ Time_bnds, & simulationStartTime, & boundaryLayerDepthSmooth, & - surfaceFluxAttenuationCoefficientRunoff) + surfaceFluxAttenuationCoefficientRunoff, & + surfaceFluxAttenuationCoefficientSubglacialRunoff) end subroutine ocn_diagnostics_variables_destroy!}}} diff --git a/components/mpas-ocean/src/shared/mpas_ocn_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_forcing.F index 99a674cd45fa..5cfce4c4e690 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_forcing.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_forcing.F @@ -27,6 +27,7 @@ module ocn_forcing use mpas_dmpar use ocn_constants use ocn_diagnostics_variables + use ocn_config implicit none private @@ -122,7 +123,7 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, forcin real (kind=RKIND) :: zTop, zBot, transmissionCoeffTop, transmissionCoeffBot - real (kind=RKIND), dimension(:,:), pointer :: layerThickness, fractionAbsorbed, fractionAbsorbedRunoff + real (kind=RKIND), dimension(:,:), pointer :: layerThickness, fractionAbsorbed, fractionAbsorbedRunoff, fractionAbsorbedSubglacialRunoff integer :: iCell, k, timeLevel, nCells @@ -147,6 +148,7 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, forcin call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', fractionAbsorbed) call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', fractionAbsorbedRunoff) + call mpas_pool_get_array(forcingPool, 'fractionAbsorbedSubglacialRunoff', fractionAbsorbedSubglacialRunoff) nCells = nCellsArray( 2 ) @@ -180,6 +182,48 @@ subroutine ocn_forcing_build_fraction_absorbed_array(meshPool, statePool, forcin end do end do +! now do subglacial runoff separately + if ( trim(config_subglacial_runoff_mode) == 'data' ) then + if ( trim(config_sgr_flux_vertical_location) == 'top' ) then + do iCell = 1, nCells + zTop = 0.0_RKIND + transmissionCoeffTop = ocn_forcing_transmission(zTop, surfaceFluxAttenuationCoefficientSubglacialRunoff(iCell)) + do k = minLevelCell(iCell), maxLevelCell(iCell) + zBot = zTop - layerThickness(k,iCell) + transmissionCoeffBot = ocn_forcing_transmission(zBot, surfaceFluxAttenuationCoefficientSubglacialRunoff(iCell)) + fractionAbsorbedSubglacialRunoff(k, iCell) = transmissionCoeffTop - transmissionCoeffBot + zTop = zBot + transmissionCoeffTop = transmissionCoeffBot + end do + end do + else if ( trim(config_sgr_flux_vertical_location) == 'uniform' ) then + do iCell = 1, nCells + ! calculate total thickness + zTop = 0.0_RKIND + do k = minLevelCell(iCell), maxLevelCell(iCell) + zTop = zTop + layerThickness(k,iCell) + end do + ! distribute flux evenly throughout water column + zBot = 0.0_RKIND + do k = minLevelCell(iCell), maxLevelCell(iCell) + fractionAbsorbedSubglacialRunoff(k, iCell) = layerThickness(k,iCell) / zTop + end do + end do + else if ( trim(config_sgr_flux_vertical_location) == 'bottom' ) then + do iCell = 1, nCells + zTop = 0.0_RKIND + transmissionCoeffTop = ocn_forcing_transmission(zTop, surfaceFluxAttenuationCoefficientSubglacialRunoff(iCell)) + do k = maxLevelCell(iCell), minLevelCell(iCell), -1 + zBot = zTop - layerThickness(k,iCell) + transmissionCoeffBot = ocn_forcing_transmission(zBot, surfaceFluxAttenuationCoefficientSubglacialRunoff(iCell)) + fractionAbsorbedSubglacialRunoff(k, iCell) = transmissionCoeffTop - transmissionCoeffBot + zTop = zBot + transmissionCoeffTop = transmissionCoeffBot + end do + end do + end if + end if + end subroutine ocn_forcing_build_fraction_absorbed_array!}}} !*********************************************************************** diff --git a/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F b/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F index e28c8f1e5819..11d70ae9c40d 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_surface_bulk_forcing.F @@ -49,7 +49,10 @@ module ocn_surface_bulk_forcing public :: ocn_surface_bulk_forcing_tracers, & ocn_surface_bulk_forcing_vel, & ocn_surface_bulk_forcing_thick, & - ocn_surface_bulk_forcing_init + ocn_surface_bulk_forcing_init, & + ocn_surface_bulk_forcing_tracers_subglacial_runoff, & + ocn_surface_bulk_forcing_thick_subglacial_runoff + !-------------------------------------------------------------------- ! @@ -128,6 +131,63 @@ subroutine ocn_surface_bulk_forcing_tracers(meshPool, groupName, forcingPool, tr end subroutine ocn_surface_bulk_forcing_tracers!}}} +!*********************************************************************** +! +! routine ocn_surface_bulk_forcing_tracers_subglacial_runoff +! +!> \brief Determines the tracers forcing array used for the bulk forcing. +!> \author Irena Vankova +!> \date July 2024 +!> \details +!> This routine computes the tracers forcing arrays used later in MPAS. +! +!----------------------------------------------------------------------- + + subroutine ocn_surface_bulk_forcing_tracers_subglacial_runoff(meshPool, groupName, forcingPool, & + tracersSurfaceFluxSubglacialRunoff, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information + character (len=*) :: groupName !< Input: Name of tracer group + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + type (mpas_pool_type), intent(inout) :: forcingPool !< Input: Forcing information + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tracersSurfaceFluxSubglacialRunoff !< Input/Output: Surface flux for tracer group due to subglacial runoff + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: Error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + err = 0 + + call mpas_timer_start("bulk_" // trim(groupName)) + if ( trim(groupName) == 'activeTracers' ) then + call ocn_surface_bulk_forcing_active_tracers_subglacial_runoff(meshPool, forcingPool, & + tracersSurfaceFluxSubglacialRunoff, err) + end if + call mpas_timer_stop("bulk_" // trim(groupName)) + + end subroutine ocn_surface_bulk_forcing_tracers_subglacial_runoff!}}} + !*********************************************************************** ! ! routine ocn_surface_bulk_forcing_vel @@ -333,8 +393,10 @@ subroutine ocn_surface_bulk_forcing_thick(forcingPool, surfaceThicknessFlux, sur !$acc riverRunoffFlux, iceRunoffFlux, rainFlux) !$acc parallel loop & - !$acc present(surfaceThicknessFlux, surfaceThicknessFluxRunoff, evaporationFlux, snowFlux, & - !$acc seaIceFreshWaterFlux, icebergFreshWaterFlux, riverRunoffFlux, iceRunoffFlux, rainFlux) + !$acc present(surfaceThicknessFlux, surfaceThicknessFluxRunoff, & + !$acc evaporationFlux, snowFlux, & + !$acc seaIceFreshWaterFlux, icebergFreshWaterFlux, riverRunoffFlux, & + !$acc iceRunoffFlux, rainFlux) #else !$omp parallel !$omp do schedule(runtime) @@ -358,6 +420,85 @@ subroutine ocn_surface_bulk_forcing_thick(forcingPool, surfaceThicknessFlux, sur end subroutine ocn_surface_bulk_forcing_thick!}}} +!*********************************************************************** +! +! routine ocn_surface_bulk_forcing_thick_subglacial_runoff +! +!> \brief Determines the thickness forcing array used for the bulk forcing. +!> \author Irena Vankova +!> \date July 2024 +!> \details +!> This routine computes the thickness forcing arrays used later in MPAS. +! +!----------------------------------------------------------------------- + + subroutine ocn_surface_bulk_forcing_thick_subglacial_runoff(forcingPool, surfaceThicknessFluxSubglacialRunoff, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + type (mpas_pool_type), intent(inout) :: forcingPool !< Input: Forcing information + real (kind=RKIND), dimension(:), intent(inout) :: & + surfaceThicknessFluxSubglacialRunoff !< Input/Output: Array for surface thickness flux due to subglacial runoff + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: Error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, nCells + + real (kind=RKIND), dimension(:), pointer :: subglacialRunoffFlux + + err = 0 + + if (bulkThicknessFluxOff) return + + call mpas_timer_start("bulk_thick", .false.) + + call mpas_pool_get_array(forcingPool, 'subglacialRunoffFlux', subglacialRunoffFlux) + + nCells = nCellsHalo( 2 ) + + ! Build surface fluxes at cell centers + +#ifdef MPAS_OPENACC + !$acc enter data copyin(subglacialRunoffFlux) + + !$acc parallel loop & + !$acc present(surfaceThicknessFluxSubglacialRunoff, subglacialRunoffFlux) +#else + !$omp parallel + !$omp do schedule(runtime) +#endif + do iCell = 1, nCells + surfaceThicknessFluxSubglacialRunoff(iCell) = subglacialRunoffFlux(iCell) / rho_sw + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef MPAS_OPENACC + !$acc exit data delete(subglacialRunoffFlux) +#endif + + call mpas_timer_stop("bulk_thick") + + end subroutine ocn_surface_bulk_forcing_thick_subglacial_runoff!}}} + + !*********************************************************************** ! ! routine ocn_surface_bulk_forcing_init @@ -470,6 +611,7 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer real (kind=RKIND), dimension(:), pointer :: rainTemperatureFlux, evapTemperatureFlux, & seaIceTemperatureFlux, icebergTemperatureFlux, & totalFreshWaterTemperatureFlux + real (kind=RKIND), dimension(:), pointer :: landIcePressure real (kind=RKIND) :: requiredSalt, allowedSalt, surfaceTemperatureFluxWithoutRunoff err = 0 @@ -505,6 +647,8 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer call mpas_pool_get_array(forcingPool, 'riverRunoffFlux', riverRunoffFlux) call mpas_pool_get_array(forcingPool, 'penetrativeTemperatureFlux', penetrativeTemperatureFlux) + call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure) + call mpas_pool_get_array(forcingPool, 'rainTemperatureFlux', rainTemperatureFlux) call mpas_pool_get_array(forcingPool, 'evapTemperatureFlux', evapTemperatureFlux) call mpas_pool_get_array(forcingPool, 'seaIceTemperatureFlux', seaIceTemperatureFlux) @@ -565,7 +709,7 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer * max(tracerGroup(index_temperature_flux,minLevelCell(iCell),iCell), 0.0_RKIND) / rho_sw ! Accumulate fluxes that use the freezing point -! mrp performance note: should call ocn_freezing_temperature just once here +! mrp performance note: should call ocn_freezing_temperature just once here seaIceTemperatureFlux(iCell) = seaIceFreshWaterFlux(iCell) * & ocn_freezing_temperature( tracerGroup(index_salinity_flux, minLevelCell(iCell), iCell), pressure=0.0_RKIND, & inLandIceCavity=.false.) / rho_sw @@ -579,7 +723,7 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer tracersSurfaceFlux(index_temperature_flux, iCell) = tracersSurfaceFlux(index_temperature_flux, iCell) & + surfaceTemperatureFluxWithoutRunoff - ! add runoff contribution for sending through coupler + ! add river runoff contribution for sending through coupler totalFreshWaterTemperatureFlux(iCell) = surfaceTemperatureFluxWithoutRunoff & + tracersSurfaceFluxRunoff(index_temperature_flux,iCell) @@ -603,6 +747,125 @@ subroutine ocn_surface_bulk_forcing_active_tracers(meshPool, forcingPool, tracer end subroutine ocn_surface_bulk_forcing_active_tracers!}}} +!*********************************************************************** +! +! routine ocn_surface_bulk_forcing_active_tracers_subglacial_runoff +! +!> \brief Determines the active tracers forcing array used for the bulk forcing. +!> \author Irena Vankova +!> \date July 2024 +!> \details +!> This routine computes the active tracers forcing arrays used later in MPAS. +! +!----------------------------------------------------------------------- + + subroutine ocn_surface_bulk_forcing_active_tracers_subglacial_runoff(meshPool, forcingPool, & + tracersSurfaceFluxSubglacialRunoff, err)!{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + type (mpas_pool_type), intent(in) :: meshPool !< Input: mesh information + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + type (mpas_pool_type), intent(inout) :: forcingPool !< Input: Forcing information + real (kind=RKIND), dimension(:,:), intent(inout) :: tracersSurfaceFluxSubglacialRunoff + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: Error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, nCells + integer, pointer :: index_temperature_fluxPtr, index_salinity_fluxPtr + integer :: index_temperature_flux, index_salinity_flux + integer, dimension(:), pointer :: nCellsArray + + type(mpas_pool_type),pointer :: tracersSurfaceFluxPool + + real (kind=RKIND), dimension(:), pointer :: subglacialRunoffFlux + real (kind=RKIND), dimension(:), pointer :: totalFreshWaterTemperatureFlux + real (kind=RKIND), dimension(:), pointer :: landIcePressure + + err = 0 + + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + + call mpas_pool_get_subpool(forcingPool, 'tracersSurfaceFlux',tracersSurfaceFluxPool) + + call mpas_pool_get_dimension(tracersSurfaceFluxPool, & + 'index_temperatureSurfaceFlux', & + index_temperature_fluxPtr) + call mpas_pool_get_dimension(tracersSurfaceFluxPool, & + 'index_salinitySurfaceFlux', & + index_salinity_fluxPtr) + index_temperature_flux = index_temperature_fluxPtr + index_salinity_flux = index_salinity_fluxPtr + + call mpas_pool_get_array(forcingPool, 'subglacialRunoffFlux', subglacialRunoffFlux) + + call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure) + + call mpas_pool_get_array(forcingPool, 'totalFreshWaterTemperatureFlux', totalFreshWaterTemperatureFlux) + + nCells = nCellsArray( 3 ) + + ! Surface fluxes of water have an associated heat content, but the coupled system does not account for this + ! Assume surface fluxes of water have a temperature dependent on the incoming mass flux. + ! Assume surface fluxes of water have zero salinity. So the RHS forcing is zero for salinity. + ! Only include this heat forcing when bulk thickness is turned on + ! indices on tracerGroup are (iTracer, iLevel, iCell) + if (.not. bulkThicknessFluxOff) then + !$omp parallel + !$omp do schedule(runtime) + do iCell = 1, nCells + + if ( config_use_sgr_opt_temp_prescribed ) then + !sgr with fixed prescribed temperature + tracersSurfaceFluxSubglacialRunoff(index_temperature_flux,iCell) = subglacialRunoffFlux(iCell) * & + config_sgr_temperature_prescribed / rho_sw + else + !sgr with temperature equal to the local freezing point of freshwater + tracersSurfaceFluxSubglacialRunoff(index_temperature_flux,iCell) = subglacialRunoffFlux(iCell) * & + ocn_freezing_temperature(salinity=0.0_RKIND, pressure=landIcePressure(iCell), & + inLandIceCavity=.true.) / rho_sw + end if + + if ( config_use_sgr_opt_salt_prescribed ) then + !sgr with fixed prescribed temperature + tracersSurfaceFluxSubglacialRunoff(index_salinity_flux,iCell) = subglacialRunoffFlux(iCell) * & + config_sgr_salinity_prescribed / rho_sw + else + !sgr with temperature equal to the local freezing point of freshwater + tracersSurfaceFluxSubglacialRunoff(index_salinity_flux,iCell) = 0.0_RKIND + end if + + ! add subglacial runoff contribution for sending through coupler + totalFreshWaterTemperatureFlux(iCell) = totalFreshWaterTemperatureFlux(iCell) & + + tracersSurfaceFluxSubglacialRunoff(index_temperature_flux,iCell) + + end do + !$omp end do + !$omp end parallel + endif ! bulkThicknessFluxOn + + end subroutine ocn_surface_bulk_forcing_active_tracers_subglacial_runoff!}}} + end module ocn_surface_bulk_forcing diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F index a16e52f14499..601d134ceec8 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F @@ -128,12 +128,14 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ ! pointers for the actual output variables within the pools above real (kind=RKIND), dimension(:), pointer, contiguous :: & surfaceThicknessFlux, &! surface thickness flux - surfaceThicknessFluxRunoff ! surface thickness flux from runoff + surfaceThicknessFluxRunoff, &! surface thickness flux from runoff + surfaceThicknessFluxSubglacialRunoff ! surface thickness flux from runoff real (kind=RKIND), dimension(:,:), pointer, contiguous :: & tendThick, &! accumulated layer thickness tendency fractionAbsorbed, &! fraction of sfc flux absorbed - fractionAbsorbedRunoff ! fraction of runoff flux absorbed + fractionAbsorbedRunoff, &! fraction of sfc flux absorbed + fractionAbsorbedSubglacialRunoff ! fraction of subglacial runoff flux absorbed !----------------------------------------------------------------- ! local variables @@ -154,10 +156,16 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ surfaceThicknessFlux) call mpas_pool_get_array(forcingPool, 'surfaceThicknessFluxRunoff', & surfaceThicknessFluxRunoff) + call mpas_pool_get_array(forcingPool, 'surfaceThicknessFluxSubglacialRunoff', & + surfaceThicknessFluxSubglacialRunoff) call mpas_pool_get_array(forcingPool, 'fractionAbsorbed', & fractionAbsorbed) call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', & fractionAbsorbedRunoff) + call mpas_pool_get_array(forcingPool, 'fractionAbsorbedSubglacialRunoff', & + fractionAbsorbedSubglacialRunoff) + + ! ! layer thickness tendency: @@ -185,6 +193,25 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ !$omp end parallel #endif + if (trim(config_subglacial_runoff_mode) == 'data') then +#ifdef MPAS_OPENACC + !$acc enter data create(surfaceThicknessFluxSubglacialRunoff) + + !$acc parallel loop & + !$acc present(surfaceThicknessFluxSubglacialRunoff) +#else + !$omp parallel + !$omp do schedule(runtime) +#endif + do iCell = 1, nCellsAll + surfaceThicknessFluxSubglacialRunoff(iCell) = 0.0_RKIND + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + end if + ! If turned off, return with zero fluxes, tendencies ! Otherwise, start time and call routines to accumulate if (config_disable_thick_all_tend) return @@ -195,6 +222,11 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ surfaceThicknessFlux, & surfaceThicknessFluxRunoff, err) + if (trim(config_subglacial_runoff_mode) == 'data') then + call ocn_surface_bulk_forcing_thick_subglacial_runoff(forcingPool, & + surfaceThicknessFluxSubglacialRunoff, err) + end if + ! Compute surface thickness flux from land ice call ocn_surface_land_ice_fluxes_thick(forcingPool, & surfaceThicknessFlux, err) @@ -217,6 +249,11 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ surfaceThicknessFlux, & surfaceThicknessFluxRunoff, & tendThick, err) + if (trim(config_subglacial_runoff_mode) == 'data') then + call ocn_thick_surface_flux_tend_subglacial_runoff(fractionAbsorbedSubglacialRunoff, & + surfaceThicknessFluxSubglacialRunoff, & + tendThick, err) + end if ! Compute contribution from frazil ice formation call ocn_frazil_forcing_layer_thickness(forcingPool, & @@ -230,7 +267,8 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{ call ocn_manufactured_solution_tend_thick(tendThick, err) #ifdef MPAS_OPENACC - !$acc exit data copyout(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff) + !$acc exit data copyout(tendThick, surfaceThicknessFlux, surfaceThicknessFluxRunoff, & + !$acc surfaceThicknessFluxSubglacialRunoff) #endif call mpas_timer_stop("ocn_tend_thick") @@ -600,9 +638,11 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & real (kind=RKIND), dimension(:,:), pointer, contiguous :: & layerThickness, &! layer thickness tracerGroupSurfaceFlux, &! tracer flux at surface - fractionAbsorbed, &! frac sfc flux aborbed - fractionAbsorbedRunoff, &! frac runoff flux aborbed + fractionAbsorbed, &! frac sfc flux absorbed + fractionAbsorbedRunoff, &! frac runoff flux absorbed + fractionAbsorbedSubglacialRunoff, &! frac subglacial runoff flux absorbed tracerGroupSurfaceFluxRunoff, &! runoff flux + tracerGroupSurfaceFluxSubglacialRunoff, &! subglacial runoff flux tracerGroupSurfaceFluxRemoved,&! total sfc flux absorbed nonLocalSurfaceTracerFlux ! non-local fluxes (eg KPP) @@ -664,6 +704,8 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & fractionAbsorbed) call mpas_pool_get_array(forcingPool, 'fractionAbsorbedRunoff', & fractionAbsorbedRunoff) + call mpas_pool_get_array(forcingPool, 'fractionAbsorbedSubglacialRunoff', & + fractionAbsorbedSubglacialRunoff) ! allocate and transfer data not specific to tracer groups allocate(normalThicknessFlux(nVertLevels, nEdgesAll+1)) @@ -788,6 +830,16 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & modifiedGroupName, & tracerGroupSurfaceFluxRunoff) + if (trim(config_subglacial_runoff_mode) == 'data') then + ! Get surface flux due to subglacial runoff array + ! only active tracers have subglacial runoff flux for now, + ! but we still need to associate for ALL tracers + modifiedGroupName = groupName // "SurfaceFluxSubglacialRunoff" + call mpas_pool_get_array(tracersSurfaceFluxPool, & + modifiedGroupName, & + tracerGroupSurfaceFluxSubglacialRunoff) + end if + ! Get surface flux removed array to keep track of how much ! flux is ignored modifiedGroupName = groupName // "SurfaceFluxRemoved" @@ -820,6 +872,18 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & !$omp end do !$omp end parallel + if (trim(config_subglacial_runoff_mode) == 'data') then + !$omp parallel + !$omp do schedule(runtime) private(n) + do iCell = 1, nCellsAll + do n=1,nTracersGroup + tracerGroupSurfaceFluxSubglacialRunoff (n,iCell) = 0.0_RKIND + end do + end do + !$omp end do + !$omp end parallel + end if + ! ! compute surface tracer flux from bulk forcing ! @@ -831,6 +895,11 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & tracerGroupSurfaceFluxRunoff, & tracerGroupSurfaceFluxRemoved, & dt, layerThickness, err) + if (trim(config_subglacial_runoff_mode) == 'data') then + call ocn_surface_bulk_forcing_tracers_subglacial_runoff(meshPool, groupName, & + forcingPool, & + tracerGroupSurfaceFluxSubglacialRunoff, err) + end if end if ! sfc bulk forcing @@ -1188,10 +1257,15 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, & endif ! compute budgets call ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, & - fractionAbsorbedRunoff, layerThickness, & - tracerGroupSurfaceFlux, & + fractionAbsorbedRunoff, & + layerThickness, tracerGroupSurfaceFlux, & tracerGroupSurfaceFluxRunoff, & tracerGroupTend, err) + if (trim(config_subglacial_runoff_mode) == 'data') then + call ocn_tracer_surface_flux_tend_subglacial_runoff(meshPool, fractionAbsorbedSubglacialRunoff, & + layerThickness, tracerGroupSurfaceFluxSubglacialRunoff, & + tracerGroupTend, err) + end if ! ! Compute shortwave absorption diff --git a/components/mpas-ocean/src/shared/mpas_ocn_thick_surface_flux.F b/components/mpas-ocean/src/shared/mpas_ocn_thick_surface_flux.F index 438aabafee17..1c535532cea2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_thick_surface_flux.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_thick_surface_flux.F @@ -46,7 +46,8 @@ module ocn_thick_surface_flux !-------------------------------------------------------------------- public :: ocn_thick_surface_flux_tend, & - ocn_thick_surface_flux_init + ocn_thick_surface_flux_init, & + ocn_thick_surface_flux_tend_subglacial_runoff !-------------------------------------------------------------------- ! @@ -69,7 +70,7 @@ module ocn_thick_surface_flux !> \date 15 September 2011 !> \details !> This routine computes the horizontal advection tendency for -!> thicknes based on current state and user choices of forcings. +!> thickness based on current state and user choices of forcings. ! !----------------------------------------------------------------------- @@ -86,8 +87,8 @@ subroutine ocn_thick_surface_flux_tend(transmissionCoefficients, transmissionCoe transmissionCoefficientsRunoff !< Input: Coefficients for the transmission of surface fluxes due to river runoff real (kind=RKIND), dimension(:), intent(in) :: & - surfaceThicknessFlux, &!< Input: surface flux of thickness - surfaceThicknessFluxRunoff !< Input: surface flux of thickness due to river runoff + surfaceThicknessFlux, &!< Input: surface flux of thickness + surfaceThicknessFluxRunoff !< Input: surface flux of thickness due to river runoff !----------------------------------------------------------------- @@ -140,7 +141,6 @@ subroutine ocn_thick_surface_flux_tend(transmissionCoefficients, transmissionCoe do k = minLevelCell(iCell), maxLevelCell(iCell) remainingFlux = remainingFlux - transmissionCoefficients(k, iCell) remainingFluxRunoff = remainingFluxRunoff - transmissionCoefficientsRunoff(k, iCell) - tend(k, iCell) = tend(k, iCell) + surfaceThicknessFlux(iCell) * transmissionCoefficients(k, iCell) & + surfaceThicknessFluxRunoff(iCell) * transmissionCoefficientsRunoff(k, iCell) end do @@ -153,6 +153,7 @@ subroutine ocn_thick_surface_flux_tend(transmissionCoefficients, transmissionCoe tend(maxLevelCell(iCell), iCell) = tend(maxLevelCell(iCell), iCell) & + remainingFluxRunoff * surfaceThicknessFluxRunoff(iCell) end if + end do #ifndef MPAS_OPENACC !$omp end do @@ -169,6 +170,107 @@ subroutine ocn_thick_surface_flux_tend(transmissionCoefficients, transmissionCoe end subroutine ocn_thick_surface_flux_tend!}}} +!*********************************************************************** +! +! routine ocn_thick_surface_flux_tend_subglacial_runoff +! +!> \brief Computes tendency term from horizontal advection of thickness +!> \author Irena Vankova +!> \date July 2024 +!> \details +!> This routine computes the horizontal advection tendency for +!> thickness based on current state and user choices of forcings. +! +!----------------------------------------------------------------------- + + subroutine ocn_thick_surface_flux_tend_subglacial_runoff(transmissionCoefficientsSubglacialRunoff, & + surfaceThicknessFluxSubglacialRunoff, tend, err)!{{{ + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(in) :: & + transmissionCoefficientsSubglacialRunoff !< Input: Coefficients for the transmission of surface fluxes due to subglacial runoff + + real (kind=RKIND), dimension(:), intent(in) :: & + surfaceThicknessFluxSubglacialRunoff !< Input: surface flux of thickness due to subglacial runoff + + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:), intent(inout) :: & + tend !< Input/Output: thickness tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, k + + real (kind=RKIND) :: remainingFlux, remainingFluxRunoff, remainingFluxSubglacialRunoff + + err = 0 + + if (.not. surfaceThicknessFluxOn) return + + call mpas_timer_start("thick surface flux") + +#ifdef MPAS_OPENACC + !$acc enter data copyin(transmissionCoefficientsSubglacialRunoff) + + !$acc parallel loop & + !$acc present(tend, & + !$acc surfaceThicknessFluxSubglacialRunoff, transmissionCoefficientsSubglacialRunoff, & + !$acc minLevelCell, maxLevelCell) & + !$acc private(k, remainingFluxSubglacialRunoff) +#else + !$omp parallel + !$omp do schedule(runtime) private(remainingFluxSubglacialRunoff, k) +#endif + do iCell = 1, nCellsOwned + remainingFluxSubglacialRunoff = 1.0_RKIND + do k = minLevelCell(iCell), maxLevelCell(iCell) + remainingFluxSubglacialRunoff = remainingFluxSubglacialRunoff - transmissionCoefficientsSubglacialRunoff(k, iCell) + tend(k, iCell) = tend(k, iCell) + surfaceThicknessFluxSubglacialRunoff(iCell) * transmissionCoefficientsSubglacialRunoff(k, iCell) + end do + + if(maxLevelCell(iCell) > 0 .and. remainingFluxSubglacialRunoff > 0.0_RKIND) then + tend(maxLevelCell(iCell), iCell) = tend(maxLevelCell(iCell), iCell) & + + remainingFluxSubglacialRunoff * surfaceThicknessFluxSubglacialRunoff(iCell) + end if + + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef MPAS_OPENACC + !$acc exit data delete(transmissionCoefficientsSubglacialRunoff) +#endif + + call mpas_timer_stop("thick surface flux") + + !-------------------------------------------------------------------- + + end subroutine ocn_thick_surface_flux_tend_subglacial_runoff!}}} + !*********************************************************************** ! ! routine ocn_thick_surface_flux_init diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tracer_surface_flux_to_tend.F b/components/mpas-ocean/src/shared/mpas_ocn_tracer_surface_flux_to_tend.F index 57e2949bac9b..ba97203c4ec3 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_tracer_surface_flux_to_tend.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_tracer_surface_flux_to_tend.F @@ -45,7 +45,8 @@ module ocn_tracer_surface_flux_to_tend !-------------------------------------------------------------------- public :: ocn_tracer_surface_flux_tend, & - ocn_tracer_surface_flux_init + ocn_tracer_surface_flux_init, & + ocn_tracer_surface_flux_tend_subglacial_runoff !-------------------------------------------------------------------- ! @@ -71,8 +72,8 @@ module ocn_tracer_surface_flux_to_tend ! !----------------------------------------------------------------------- - subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbsorbedRunoff, layerThickness, & - surfaceTracerFlux, surfaceTracerFluxRunoff, tend, err)!{{{ + subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbsorbedRunoff, & + layerThickness, surfaceTracerFlux, surfaceTracerFluxRunoff, tend, err)!{{{ !----------------------------------------------------------------- ! ! input variables @@ -197,11 +198,120 @@ subroutine ocn_tracer_surface_flux_tend(meshPool, fractionAbsorbed, fractionAbso call mpas_timer_stop("surface_tracer_runoff_flux") end if - !-------------------------------------------------------------------- end subroutine ocn_tracer_surface_flux_tend!}}} + +!*********************************************************************** +! +! routine ocn_tracer_surface_flux_tend_subglacial_runoff +! +!> \brief Computes tendency term for surface fluxes +!> \author Irena Vankova +!> \date July 2024 +!> \details +!> This routine computes the tendency for tracers based on surface fluxes. +! +!----------------------------------------------------------------------- + + subroutine ocn_tracer_surface_flux_tend_subglacial_runoff(meshPool, fractionAbsorbedSubglacialRunoff, & + layerThickness, surfaceTracerFluxSubglacialRunoff, tend, err)!{{{ + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: & + meshPool !< Input: mesh information + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness !< Input: Layer thickness + + real (kind=RKIND), dimension(:,:), intent(in), pointer :: & + surfaceTracerFluxSubglacialRunoff !< Input: surface tracer fluxes from subglacial runoff + + real (kind=RKIND), dimension(:,:), intent(in) :: & + fractionAbsorbedSubglacialRunoff !< Input: Coefficients for the application of surface fluxes due to subglacial runoff + + + !----------------------------------------------------------------- + ! + ! input/output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + tend !< Input/Output: velocity tendency + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + integer, intent(out) :: err !< Output: error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + integer :: iCell, k, iTracer, nTracers, nCells + integer, pointer :: nVertLevels + integer, dimension(:), pointer :: nCellsArray + integer, dimension(:), pointer :: minLevelCell, maxLevelCell + + real (kind=RKIND) :: remainingFlux + + err = 0 + + if (.not. surfaceTracerFluxOn) return + + call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + nTracers = size(tend, dim=1) + + call mpas_pool_get_array(meshPool, 'minLevelCell', minLevelCell) + call mpas_pool_get_array(meshPool, 'maxLevelCell', maxLevelCell) + + nCells = nCellsArray( 1 ) + + ! now do subglacial runoff component + + call mpas_timer_start("surface_tracer_subglacial_runoff_flux") + + !$omp parallel + !$omp do schedule(runtime) private(remainingFlux, k, iTracer) + do iCell = 1, nCells + remainingFlux = 1.0_RKIND + do k = minLevelCell(iCell), maxLevelCell(iCell) + remainingFlux = remainingFlux - fractionAbsorbedSubglacialRunoff(k, iCell) + + do iTracer = 1, nTracers + tend(iTracer, k, iCell) = tend(iTracer, k, iCell) + & + surfaceTracerFluxSubglacialRunoff(iTracer, iCell) * fractionAbsorbedSubglacialRunoff(k, iCell) + end do + end do + + if(maxLevelCell(iCell) > 0 .and. remainingFlux > 0.0_RKIND) then + do iTracer = 1, nTracers + tend(iTracer, maxLevelCell(iCell), iCell) = tend(iTracer, maxLevelCell(iCell), iCell) & + + surfaceTracerFluxSubglacialRunoff(iTracer, iCell) * remainingFlux + end do + end if + end do + !$omp end do + !$omp end parallel + + call mpas_timer_stop("surface_tracer_subglacial_runoff_flux") + + !-------------------------------------------------------------------- + + end subroutine ocn_tracer_surface_flux_tend_subglacial_runoff!}}} + !*********************************************************************** ! ! routine ocn_tracer_surface_flux_init diff --git a/components/mpas-ocean/src/tracer_groups/Registry_CFC.xml b/components/mpas-ocean/src/tracer_groups/Registry_CFC.xml index d77859dcef86..b0adb7d595f5 100644 --- a/components/mpas-ocean/src/tracer_groups/Registry_CFC.xml +++ b/components/mpas-ocean/src/tracer_groups/Registry_CFC.xml @@ -95,6 +95,15 @@ description="CFC12 Surface Flux Due to Runoff" /> + + + + + + + + + + + + + + + + + + + + + + + + + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/components/mpas-ocean/src/tracer_groups/Registry_idealAge.xml b/components/mpas-ocean/src/tracer_groups/Registry_idealAge.xml index c6adca3ea636..3a45bc10b794 100644 --- a/components/mpas-ocean/src/tracer_groups/Registry_idealAge.xml +++ b/components/mpas-ocean/src/tracer_groups/Registry_idealAge.xml @@ -69,6 +69,11 @@ description="Flux of iAge through the ocean surface due to river runoff. Positive into ocean." /> + + +