diff --git a/components/mpas-albany-landice/driver/glc_comp_mct.F b/components/mpas-albany-landice/driver/glc_comp_mct.F index fb4a6fdae2a3..fe0749056831 100644 --- a/components/mpas-albany-landice/driver/glc_comp_mct.F +++ b/components/mpas-albany-landice/driver/glc_comp_mct.F @@ -914,9 +914,9 @@ subroutine glc_run_mct( EClock, cdata_g, x2g_g, g2x_g)!{{{ ! Finally, set whence to latest_before so we have piecewise-constant forcing. ! Could add, e.g., linear interpolation later. - ! NOTE: If any input ("forcing") files here contain fields that the coupler also passes, - ! the input file versions here will clobber the fields passed from the coupler. - + ! NOTE: If any input ("forcing") files here contain fields that the coupler also passes, + ! the input file versions here will clobber the fields passed from the coupler. + call mpas_stream_mgr_read(domain % streamManager, whence=MPAS_STREAM_LATEST_BEFORE, saveActualWhen = .true., ierr=err_tmp) err = ior(err, err_tmp) call mpas_stream_mgr_reset_alarms(domain % streamManager, direction=MPAS_STREAM_INPUT, ierr=err_tmp) @@ -1410,11 +1410,13 @@ subroutine glc_import_mct(x2g_g, errorCode) real (kind=RKIND), dimension(:), pointer :: sfcMassBal,& floatingBasalMassBal,& surfaceTemperature,& - basalOceanHeatflx,& - OceanDensity + basalOceanHeatFlux,& + iceOceanInterfaceTemperature + real (kind=RKIND), dimension(:,:), pointer :: ismip6shelfMelt_3dThermalForcing integer, dimension(:,:), pointer :: orig3dOceanMask - real(kind=RKIND), pointer :: config_invalid_value_TF + real(kind=RKIND), pointer :: config_invalid_value_TF, config_min_floating_frac, & + config_invalid_value_ocn_intr_temp character(len=StrKIND), pointer :: config_smb_source real(kind=RKIND) :: fractionalMaskVal @@ -1422,6 +1424,10 @@ subroutine glc_import_mct(x2g_g, errorCode) call mpas_pool_get_config(domain % configs, 'config_invalid_value_TF', config_invalid_value_TF) call mpas_pool_get_config(domain % configs, 'config_smb_source', config_smb_source) + call mpas_pool_get_config(domain % configs, 'config_min_floating_frac', config_min_floating_frac) + call mpas_pool_get_config(domain % configs, 'config_min_floating_frac', config_min_floating_frac) + call mpas_pool_get_config(domain % configs, 'config_invalid_value_ocn_intr_temp', & + config_invalid_value_ocn_intr_temp) n = 0 block => domain % blocklist @@ -1440,10 +1446,10 @@ subroutine glc_import_mct(x2g_g, errorCode) call mpas_pool_get_array(geometryPool, 'sfcMassBal', sfcMassBal) call mpas_pool_get_array(geometryPool, 'floatingBasalMassBal',floatingBasalMassBal) call mpas_pool_get_array(thermalPool, 'surfaceTemperature',surfaceTemperature) + call mpas_pool_get_array(thermalPool, 'basalOceanHeatFlux', basalOceanHeatFlux) call mpas_pool_get_array(geometryPool, 'ismip6shelfMelt_3dThermalForcing', ismip6shelfMelt_3dThermalForcing) call mpas_pool_get_array(extrapOceanDataPool, 'orig3dOceanMask', orig3dOceanMask) -! call mpas_pool_get_array(thermalPool, 'basalOceanHeatflx',basalOceanHeatflx) - !call mpas_pool_get_array(geometryPool, 'OceanDensity',OceanDensity) + call mpas_pool_get_array(thermalPool, 'iceOceanInterfaceTemperature', iceOceanInterfaceTemperature) if (nISMIP6OceanLayers > 0) then orig3dOceanMask(:,:) = 0 @@ -1459,7 +1465,8 @@ subroutine glc_import_mct(x2g_g, errorCode) call mpas_log_write("Unknown value for 'config_smb_source'", MPAS_LOG_ERR) errorCode = ior(errorCode, 1) endif - floatingBasalMassBal(i) = x2g_g % rAttr(index_x2g_Fogx_qiceli, n) + floatingBasalMassBal(i) = -(x2g_g % rAttr(index_x2g_Foxo_ismw, n) & + + x2g_g % rAttr(index_x2g_Foxo_frazil_li, n)) if (nISMIP6OceanLayers > 0) then do iLev = 1, nISMIP6OceanLayers fractionalMaskVal = x2g_g % rAttr(index_x2g_So_tf3d_mask(iLev), n) @@ -1475,10 +1482,14 @@ subroutine glc_import_mct(x2g_g, errorCode) endif enddo endif -! surfaceTemperature(i) = x2g_g % rAttr(index_x2g_Sl_tsrf, n) -!JW basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogo_qiceh, n) -! basalOceanHeatflx(i) = x2g_g % rAttr(index_x2g_Fogx_qicehi, n) - !OceanDensity (i) = x2g_g % rAttr(, n) + basalOceanHeatFlux(i) = -x2g_g % rAttr(index_x2g_Fogx_qicehi, n) + fractionalMaskVal = x2g_g % rAttr(index_x2g_So_liflfrac, n) + if (fractionalMaskVal > config_min_floating_frac) then + iceOceanInterfaceTemperature(i) = x2g_g % rAttr(index_x2g_So_tfrz_isf, n) & + / fractionalMaskVal + else + iceOceanInterfaceTemperature(i) = config_invalid_value_ocn_intr_temp + endif end do block => block % next @@ -1516,13 +1527,12 @@ subroutine glc_export_mct(g2x_g, errorCode) integer, pointer :: nCellsSolve, nVertLevels real (kind=RKIND), dimension(:), pointer :: upperSurface - real (kind=RKIND), dimension(:), pointer :: layerThicknessFractions real (kind=RKIND), dimension(:), pointer :: thickness + real (kind=RKIND), dimension(:), pointer :: basalConductiveFlux real (kind=RKIND), dimension(:), pointer :: avgBareIceAblationApplied real (kind=RKIND), dimension(:), pointer :: avgCalvingFlux real (kind=RKIND), dimension(:), pointer :: avgFaceMeltFlux real (kind=RKIND), dimension(:), pointer :: avgFloatingBMBFlux - real (kind=RKIND), dimension(:,:), pointer :: temperature integer, dimension(:), pointer :: cellMask !------------------------------------------------------------------- @@ -1545,9 +1555,8 @@ subroutine glc_export_mct(g2x_g, errorCode) call mpas_pool_get_subpool(block % structs, 'timeAveraging', timeAveragingPool) call mpas_pool_get_array(geometryPool, 'upperSurface', upperSurface) - call mpas_pool_get_array(meshPool, 'layerThicknessFractions', layerThicknessFractions) call mpas_pool_get_array(geometryPool, 'thickness', Thickness, timeLevel = 1) - call mpas_pool_get_array(thermalPool, 'temperature', temperature) + call mpas_pool_get_array(thermalPool, 'basalConductiveFlux', basalConductiveFlux) call mpas_pool_get_array(timeAveragingPool, 'avgBareIceAblationApplied', avgBareIceAblationApplied) call mpas_pool_get_array(timeAveragingPool, 'avgCalvingFlux', avgCalvingFlux) @@ -1560,11 +1569,11 @@ subroutine glc_export_mct(g2x_g, errorCode) n = n + 1 ! Fogg_rofl - g2x_g % rAttr(index_g2x_Fogg_rofl,n) = avgBareIceAblationApplied(i) + g2x_g % rAttr(index_g2x_Fogg_rofl,n) = avgBareIceAblationApplied(i) ! Figg_rofi g2x_g % rAttr(index_g2x_Figg_rofi,n) = 0.0 ! placeholder ! Fogg_rofi - g2x_g % rAttr(index_g2x_Fogg_rofi,n) = avgCalvingFlux(i) + g2x_g % rAttr(index_g2x_Fogg_rofi,n) = avgCalvingFlux(i) g2x_g % rAttr(index_g2x_Fogg_rofi,n) = g2x_g % rAttr(index_g2x_Fogg_rofi,n) + avgFaceMeltFlux(i) if (trim(config_basal_mass_bal_float) == 'ismip6') then ! if MALI is calculating ISMF, add that to rofi @@ -1573,10 +1582,9 @@ subroutine glc_export_mct(g2x_g, errorCode) g2x_g % rAttr(index_g2x_Fogg_rofi,n) = g2x_g % rAttr(index_g2x_Fogg_rofi,n) - avgFloatingBMBFlux(i) endif + g2x_g % rAttr(index_g2x_Flgg_hflx, n) = basalConductiveFlux(i) g2x_g % rAttr(index_g2x_Sg_topo, n) = max(0.0, upperSurface(i)) !updated to avoid warning for values below sea level - g2x_g % rAttr(index_g2x_Sg_tbot, n) = temperature(nVertlevels,i) - SHR_CONST_TKTRIP ! layerThickness is not populated on init, so calculating like this instead: - g2x_g % rAttr(index_g2x_Sg_dztbot, n) = layerThicknessFractions(nVertLevels) * thickness(i) / 2.0 g2x_g % rAttr(index_g2x_Sg_lithop, n) = Thickness(i)*SHR_CONST_RHOICE*SHR_CONST_G !!!Mask configurations @@ -1648,10 +1656,10 @@ subroutine convert_seconds_to_timestamp(seconds, timeStamp)!{{{ character (len=StrKIND), intent(out) :: timeStamp real (kind=RKIND) :: secondsPerHour, secondsPerMinute, remaining integer :: minutes, hours, secondsLeft - + secondsPerHour = 3600 secondsPerMinute = 60 - + if(seconds < 0 .or. seconds > 86400) then secondsLeft = 00 minutes = 00 @@ -1659,16 +1667,16 @@ subroutine convert_seconds_to_timestamp(seconds, timeStamp)!{{{ else hours = int(seconds/secondsPerHour) remaining = seconds - real(hours) * secondsPerHour - + minutes = int(remaining/secondsPerMinute) remaining = remaining - real(minutes) * secondsPerMinute - + secondsLeft = int(remaining) end if - + write(timeStamp,"(a,i2.2,a,i2.2,a,i2.2)") "00_",hours,":",minutes,":",secondsLeft timeStamp = trim(timeStamp) - + end subroutine convert_seconds_to_timestamp!}}} subroutine add_stream_attributes(stream_manager, domain)!{{{ diff --git a/components/mpas-albany-landice/driver/glc_cpl_indices.F b/components/mpas-albany-landice/driver/glc_cpl_indices.F index 0bd17918ef00..30acb1efadc1 100644 --- a/components/mpas-albany-landice/driver/glc_cpl_indices.F +++ b/components/mpas-albany-landice/driver/glc_cpl_indices.F @@ -1,5 +1,5 @@ module glc_cpl_indices - + use seq_flds_mod use mct_mod ! use glc_constants, only : glc_nec, glc_smb ! TODO Will these be needed? If so, need to add MPAS version @@ -17,28 +17,20 @@ module glc_cpl_indices integer, public :: index_x2g_Flgl_qice = 0 !Ice sheet surface mass balance integer, public :: index_x2g_Sl_tsrf = 0 !Ice sheet upper surface boundary temperature - integer, public :: index_x2g_So_blt = 0 !Ice shelf boundary layer ocean temperature - integer, public :: index_x2g_So_bls = 0 !Ice shelf boundary layer ocean salinity - integer, public :: index_x2g_So_htv = 0 !Ice shelf ocean heat transfer velocity - integer, public :: index_x2g_So_stv = 0 !Ice shelf ocean salinity transfer velocity - integer, public :: index_x2g_So_rhoeff = 0 !Ocean effective pressure + integer, public :: index_x2g_So_tfrz_isf = 0 !Ice shelf-ocean interface temperature + integer, public :: index_x2g_So_liflfrac = 0 !Ice shelf fractional coverage from ocean integer, public :: index_x2g_So_tf3d(glc_nzoc_max) = 0 !Ocean thermal forcing at specified z-levels integer, public :: index_x2g_So_tf3d_mask(glc_nzoc_max) = 0 !mask of ocean thermal forcing at specified z-levels - integer, public :: index_x2g_Fogx_qiceli = 0 !Subshelf mass flux integer, public :: index_x2g_Fogx_qicehi = 0 !Subshelf heat flux for the ice sheet + integer, public :: index_x2g_Foxo_ismw = 0 !Ice shelf meltwater flux from ocean + integer, public :: index_x2g_Foxo_frazil_li = 0 !Frazil flux under ice shelves from ocean integer, public :: index_g2x_Fogg_rofi = 0 ! frozen runoff -> ocn integer, public :: index_g2x_Figg_rofi = 0 ! frozen runoff -> ice integer, public :: index_g2x_Fogg_rofl = 0 ! liquid runoff -> ocn integer, public :: index_g2x_Sg_topo = 0 ! surface elevation integer, public :: index_g2x_Flgg_hflx = 0 ! heat flux - integer, public :: index_g2x_Sg_tbot = 0 !bottom-layer temperature - integer, public :: index_g2x_Sg_dztbot = 0 !distance from bottom-layer temperature to ice-ocean interface (to compute T gradient) integer, public :: index_g2x_Sg_lithop = 0 !ice sheet lithostatic pressure - integer, public :: index_g2x_Fogx_qicelo = 0 !Subshelf liquid flux for ocean - integer, public :: index_g2x_Fogx_qiceho = 0 !Subshelf heat flux for the ocean - integer, public :: index_g2x_Sg_blis = 0 !Boundary layer interface salinity for ocean - integer, public :: index_g2x_Sg_blit = 0 !Boundary layer interface temperature for ocean integer, public :: index_g2x_Sg_icemask = 0 !complete grounded/floating ice mask integer, public :: index_g2x_Sg_ice_covered = 0 !ice covered ice mask integer, public :: index_g2x_Sg_icemask_grounded = 0 !grounded mask @@ -75,9 +67,12 @@ subroutine glc_cpl_indices_set( ) index_x2g_Sl_tsrf = mct_avect_indexra(x2g,'Sl_tsrf',perrwith='quiet') ! Surface temperature (deg C) - index_x2g_Fogx_qiceli = mct_avect_indexra(x2g,'Fogx_qiceli',perrwith='quiet') index_x2g_Fogx_qicehi = mct_avect_indexra(x2g,'Fogx_qicehi',perrwith='quiet') - index_x2g_So_rhoeff = mct_avect_indexra(x2g,'So_rhoeff',perrwith='quiet') + + ! Sub-ice-shelf fluxes from ocean + index_x2g_Foxo_ismw = mct_avect_indexra(x2g,'Foxo_ismw',perrwith='quiet') + index_x2g_Foxo_frazil_li = mct_avect_indexra(x2g,'Foxo_frazil_li',perrwith='quiet') + if (glc_nzoc > 0) then do iLev = 1, glc_nzoc @@ -92,12 +87,6 @@ subroutine glc_cpl_indices_set( ) !Following block of x2g/g2x vectors are used internally within coupler for subshelf melt flux !calculations (and so do not have directly-related export-side arrays) - index_x2g_So_blt = mct_avect_indexra(x2g,'So_blt',perrwith='quiet') - index_x2g_So_bls = mct_avect_indexra(x2g,'So_bls',perrwith='quiet') - index_x2g_So_htv = mct_avect_indexra(x2g,'So_htv',perrwith='quiet') - index_x2g_So_stv = mct_avect_indexra(x2g,'So_stv',perrwith='quiet') - index_g2x_Sg_tbot = mct_avect_indexra(g2x,'Sg_tbot',perrwith='quiet') - index_g2x_Sg_dztbot = mct_avect_indexra(g2x,'Sg_dztbot',perrwith='quiet') index_g2x_Sg_lithop = mct_avect_indexra(g2x,'Sg_lithop',perrwith='quiet') !Following block are GLC outputs for other components. @@ -151,11 +140,9 @@ subroutine glc_cpl_indices_set( ) ! but not evolving. The GLC_TWO_WAY_COUPLING xml variable is what controls ! how this should be handled in CLM and the ice sheet model. - !Following block are subshelf melt routine outputs for ocean import - index_g2x_Fogx_qicelo = mct_avect_indexra(g2x,'Fogx_qicelo',perrwith='quiet') - index_g2x_Fogx_qiceho = mct_avect_indexra(g2x,'Fogx_qiceho',perrwith='quiet') - index_g2x_Sg_blis = mct_avect_indexra(g2x,'Sg_blis',perrwith='quiet') - index_g2x_Sg_blit = mct_avect_indexra(g2x,'Sg_blit',perrwith='quiet') + ! ice shelf-ocean interface temperature from the ocean + index_x2g_So_tfrz_isf = mct_avect_indexra(x2g,'So_tfrz_isf',perrwith='quiet') + index_x2g_So_liflfrac = mct_avect_indexra(x2g,'So_liflfrac',perrwith='quiet') call mct_aVect_clean(x2g) call mct_aVect_clean(g2x) diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 3c4b562c4e3d..abee26859dee 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -490,11 +490,11 @@ /> + /> + + @@ -1564,10 +1572,10 @@ is the value of that variable from the *previous* time level! - - @@ -1772,12 +1780,18 @@ is the value of that variable from the *previous* time level! + + diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F index acc5048e3a4e..822f7b9a71ad 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_thermal.F @@ -133,7 +133,8 @@ subroutine li_thermal_init(domain, err) real (kind=RKIND), pointer :: & config_surface_air_temperature_value, & ! constant value of surface air temperature config_surface_air_temperature_lapse_rate, & ! optional lapse rate to apply to constant value of surface air temperature - config_basal_heat_flux_value ! constant value of basal heat flux + config_basal_heat_flux_value, & ! constant value of basal heat flux + config_invalid_value_ocn_intr_temp ! invalid ice-ocean interface temperature sentinel integer, pointer :: & config_stats_cell_ID ! cell ID for diagnostic output @@ -162,7 +163,9 @@ subroutine li_thermal_init(domain, err) real (kind=RKIND), dimension(:), pointer :: & surfaceAirTemperature, & ! surface air temperature (K) surfaceTemperature, & ! surface ice temperature (K) - basalTemperature ! basal ice temperature (K) + basalTemperature, & ! basal ice temperature (K) + iceOceanInterfaceTemperature, & ! ice-ocean interface temperature (K) + basalOceanHeatFlux ! ocean-driven basal heat flux for floating ice (W m^{-2}, positive down) real (kind=RKIND), dimension(:), pointer :: & basalHeatFlux ! basal heat flux into the ice (W m^{-2}, positive upward) @@ -248,6 +251,7 @@ subroutine li_thermal_init(domain, err) call mpas_pool_get_array(thermalPool, 'surfaceAirTemperature', surfaceAirTemperature) call mpas_pool_get_array(thermalPool, 'surfaceTemperature', surfaceTemperature) call mpas_pool_get_array(thermalPool, 'basalTemperature', basalTemperature) + call mpas_pool_get_array(thermalPool, 'basalOceanHeatFlux', basalOceanHeatFlux) call mpas_pool_get_array(thermalPool, 'basalHeatFlux', basalHeatFlux) if (config_print_thermal_info) then @@ -276,6 +280,8 @@ subroutine li_thermal_init(domain, err) call mpas_log_write('Initialize basal heat flux: $r', realArgs=(/config_basal_heat_flux_value/)) endif + basalOceanHeatFlux(:) = 0.0_RKIND + ! Precompute some grid quantities used in the vertical temperature solve allocate(dsigmaTerm(nVertLevels,2)) @@ -709,7 +715,9 @@ subroutine li_thermal_solver(domain, err) real (kind=RKIND), dimension(:), pointer :: & surfaceTemperature, & ! surface ice temperature (K) basalTemperature, & ! basal ice temperature (K) + iceOceanInterfaceTemperature, & ! ice-ocean interface temperature (K) surfaceAirTemperature, & ! surface air temperature (K) + basalOceanHeatFlux, & ! ocean-driven basal heat flux for floating ice (W m^{-2}, positive down) basalHeatFlux, & ! basal heat flux into the ice (W m^{-2}, positive upward) basalFrictionFlux, & ! basal frictional flux into the ice (W m^{-2}) surfaceConductiveFlux, & ! conductive heat flux at the upper surface (W m^{-2}, positive down) @@ -745,6 +753,7 @@ subroutine li_thermal_solver(domain, err) dTtop, dTbot, & ! temperature differences denth_top, denth_bot, & ! enthalpy differences columnHeatDissipation, & ! integrated heat dissipation in column + appliedBasalOceanHeatFlux, & ! imported ocean heat flux applied to floating ice maxtemp, mintemp, & ! max and min temperatures in column initialEnergy, & ! initial energy in ice column (J m^{-2}) finalEnergy, & ! final energy in ice column (J m^{-2}) @@ -834,9 +843,11 @@ subroutine li_thermal_solver(domain, err) call mpas_pool_get_array(thermalPool, 'drainedInternalMeltRate', drainedInternalMeltRate) call mpas_pool_get_array(thermalPool, 'surfaceTemperature', surfaceTemperature) call mpas_pool_get_array(thermalPool, 'basalTemperature', basalTemperature) + call mpas_pool_get_array(thermalPool, 'iceOceanInterfaceTemperature', iceOceanInterfaceTemperature) call mpas_pool_get_array(thermalPool, 'surfaceAirTemperature', surfaceAirTemperature) call mpas_pool_get_array(thermalPool, 'surfaceConductiveFlux', surfaceConductiveFlux) call mpas_pool_get_array(thermalPool, 'basalConductiveFlux', basalConductiveFlux) + call mpas_pool_get_array(thermalPool, 'basalOceanHeatFlux', basalOceanHeatFlux) call mpas_pool_get_array(thermalPool, 'basalHeatFlux', basalHeatFlux) call mpas_pool_get_array(thermalPool, 'basalFrictionFlux', basalFrictionFlux) call mpas_pool_get_array(thermalPool, 'heatDissipation', heatDissipation) @@ -851,6 +862,7 @@ subroutine li_thermal_solver(domain, err) call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) call mpas_pool_get_config(liConfigs, 'config_thermal_thickness', config_thermal_thickness) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) + call mpas_pool_get_config(liConfigs, 'config_invalid_value_ocn_intr_temp', config_invalid_value_ocn_intr_temp) call mpas_pool_get_config(liConfigs, 'config_stats_cell_ID', config_stats_cell_ID) if (config_print_thermal_info) then @@ -944,11 +956,16 @@ subroutine li_thermal_solver(domain, err) surfaceTemperature(iCell) = min(0.0_RKIND, surfaceAirTemperature(iCell)) - ! For floating ice, set the basal temperature to the freezing temperature of seawater. - ! Values based on Ocean Water Freezing Point Calculator with S = 35 PSU + ! For floating ice, keep the thermal solve on an imported ocean-temperature boundary. + ! The merged ocean heat flux remains separate from geothermal forcing and is applied + ! directly in floating melt/freeze bookkeeping. if (li_mask_is_floating_ice(cellMask(iCell))) then - depth = thickness(iCell) * rhoi/rhoo - basalTemperature(iCell) = oceanFreezingTempSurface + oceanFreezingTempDepthDependence * depth ! Celsius + if (iceOceanInterfaceTemperature(iCell) < 0.5_RKIND * config_invalid_value_ocn_intr_temp) then + basalTemperature(iCell) = iceOceanInterfaceTemperature(iCell) - kelvin_to_celsius + else + depth = thickness(iCell) * rhoi/rhoo + basalTemperature(iCell) = oceanFreezingTempSurface + oceanFreezingTempDepthDependence * depth ! Celsius + endif endif if (trim(config_thermal_solver) == 'enthalpy') then @@ -1198,6 +1215,10 @@ subroutine li_thermal_solver(domain, err) + heatDissipation(k,iCell) * (layerInterfaceSigma(k+1) - layerInterfaceSigma(k)) enddo columnHeatDissipation = columnHeatDissipation * thickness(iCell)*rhoi*cp_ice + appliedBasalOceanHeatFlux = 0.0_RKIND + if (li_mask_is_floating_ice(cellMask(iCell))) then + appliedBasalOceanHeatFlux = basalOceanHeatFlux(iCell) + endif if (verboseColumn) then call mpas_log_write(' ') @@ -1234,6 +1255,12 @@ subroutine li_thermal_solver(domain, err) call mpas_log_write('Basal fluxes:') call mpas_log_write('frictional =$r', realArgs=(/basalFrictionFlux(iCell)/)) call mpas_log_write('geothermal =$r', realArgs=(/basalHeatFlux(iCell)/)) + if (li_mask_is_floating_ice(cellMask(iCell))) then + call mpas_log_write('imported ocean basal heat flux (positive down) =$r', & + realArgs=(/basalOceanHeatFlux(iCell)/)) + call mpas_log_write('applied ocean basal heat flux =$r', & + realArgs=(/appliedBasalOceanHeatFlux/)) + endif call mpas_log_write('flux for bottom melting =$r', & realArgs=(/basalFrictionFlux(iCell) + basalHeatFlux(iCell) + basalConductiveFlux(iCell)/)) endif ! verboseColumn @@ -3156,4 +3183,3 @@ end subroutine tridiag_solver end module li_thermal !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| - diff --git a/components/mpas-ocean/bld/build-namelist b/components/mpas-ocean/bld/build-namelist index 11173ca120c2..dc8f5ce54b72 100755 --- a/components/mpas-ocean/bld/build-namelist +++ b/components/mpas-ocean/bld/build-namelist @@ -854,6 +854,7 @@ if ($OCN_GLC_ISMF_COUPLING eq 'coupler') { # Namelist group: land_ice # ############################ +add_default($nl, 'config_use_land_ice_pressure'); add_default($nl, 'config_land_ice_draft_mode'); add_default($nl, 'config_land_ice_rho_ocean'); @@ -862,10 +863,10 @@ add_default($nl, 'config_land_ice_rho_ocean'); ################################### if ($OCN_GLC_ISMF_COUPLING eq 'coupler') { - add_default($nl, 'config_land_ice_flux_mode', 'val'=>"coupled"); + add_default($nl, 'config_land_ice_flux_mode', 'val'=>"active"); add_default($nl, 'config_frazil_under_land_ice', 'val'=>".true."); } elsif ($OCN_GLC_ISMF_COUPLING eq 'internal_mpaso') { - add_default($nl, 'config_land_ice_flux_mode', 'val'=>"standalone"); + add_default($nl, 'config_land_ice_flux_mode', 'val'=>"active"); add_default($nl, 'config_frazil_under_land_ice', 'val'=>".true."); } elsif ($OCN_GLC_ISMF_COUPLING eq 'data_mpaso') { add_default($nl, 'config_land_ice_flux_mode', 'val'=>"data"); diff --git a/components/mpas-ocean/bld/build-namelist-section b/components/mpas-ocean/bld/build-namelist-section index b36e1c3f8ddb..2290a80d9f2a 100755 --- a/components/mpas-ocean/bld/build-namelist-section +++ b/components/mpas-ocean/bld/build-namelist-section @@ -330,6 +330,7 @@ add_default($nl, 'config_frazil_use_surface_pressure'); # Namelist group: land_ice # ############################ +add_default($nl, 'config_use_land_ice_pressure'); add_default($nl, 'config_land_ice_draft_mode'); add_default($nl, 'config_land_ice_rho_ocean'); 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 a20265916254..84295e25947f 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_defaults_mpaso.xml @@ -448,6 +448,22 @@ .false. +.false. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. +.true. 'pressure-dependent' 'data' 'data' @@ -463,21 +479,6 @@ 'off' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' -'pressure_only' 'Jenkins' .false. 10.0 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 42d7e4e85f4d..2dd2f1879393 100644 --- a/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml +++ b/components/mpas-ocean/bld/namelist_files/namelist_definition_mpaso.xml @@ -1740,6 +1740,14 @@ Default: Defined in namelist_defaults.xml + +If .true. then applies landIcePressure field + +Valid values: .true. or .false. +Default: Defined in namelist_defaults.xml + + Selects the mode in which land-ice draft is computed. @@ -1750,7 +1758,7 @@ Default: Defined in namelist_defaults.xml -ocean density used to calculate landIceDraft at floatation (assumed constant and uniform). Should be consistent with MALI's config_ocean_density when used to determine grounding line location. This is an alternative to the coupler variable effectiveDensityInLandIce which is not currently used. +ocean density used to calculate landIceDraft at floatation (assumed constant and uniform). Should be consistent with MALI's config_ocean_density when used to determine grounding line location. Valid values: Any positive real number Default: Defined in namelist_defaults.xml @@ -1763,7 +1771,7 @@ Default: Defined in namelist_defaults.xml category="land_ice_fluxes" group="land_ice_fluxes"> Selects the mode in which land-ice fluxes are computed. -Valid values: 'off','pressure_only', 'data', 'standalone','coupled' +Valid values: 'off', 'data', 'active' Default: Defined in namelist_defaults.xml diff --git a/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glcshelf/user_nl_mpaso b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glcshelf/user_nl_mpaso index 349facb8caaa..ae020e99a4b2 100644 --- a/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glcshelf/user_nl_mpaso +++ b/components/mpas-ocean/cime_config/testdefs/testmods_dirs/mpaso/ocn_glcshelf/user_nl_mpaso @@ -1,2 +1,2 @@ - config_land_ice_flux_mode = 'coupled' + config_land_ice_flux_mode = 'active' config_check_ssh_consistency = .false. diff --git a/components/mpas-ocean/driver/mpaso_cpl_indices.F b/components/mpas-ocean/driver/mpaso_cpl_indices.F index de86b439c2ea..5bf1937024f1 100644 --- a/components/mpas-ocean/driver/mpaso_cpl_indices.F +++ b/components/mpas-ocean/driver/mpaso_cpl_indices.F @@ -1,5 +1,5 @@ module mpaso_cpl_indices - + use seq_flds_mod use mct_mod @@ -12,7 +12,7 @@ module mpaso_cpl_indices ! ocn -> drv - integer :: index_o2x_So_t + integer :: index_o2x_So_t integer :: index_o2x_So_u integer :: index_o2x_So_v integer :: index_o2x_So_s @@ -35,11 +35,8 @@ module mpaso_cpl_indices ! ocn -> drv for calculation of ocean-ice sheet interactions - integer :: index_o2x_So_blt !boundary layer temperature - integer :: index_o2x_So_bls !boundary layer salinity - integer :: index_o2x_So_htv !ocean heat-transfer velocity - integer :: index_o2x_So_stv !ocean salt-transfer velocity - integer :: index_o2x_So_rhoeff !ocean effective density + integer :: index_o2x_So_tfrz_isf !boundary layer temperature + integer :: index_o2x_So_liflfrac !floating ice shelf fraction integer :: index_o2x_So_tf3d(glc_nzoc_max) !ocean thermal forcing at predefined z-levels integer :: index_o2x_So_tf3d_mask(glc_nzoc_max) !mask ofocean thermal forcing at predefined z-levels @@ -84,7 +81,7 @@ module mpaso_cpl_indices integer :: index_x2o_Foxx_tauy ! meridonal wind stress (tauy) (W/m2 ) integer :: index_x2o_Foxx_swnet ! net short-wave heat flux (W/m2 ) integer :: index_x2o_Foxx_sen ! sensible heat flux (W/m2 ) - integer :: index_x2o_Foxx_lat + integer :: index_x2o_Foxx_lat integer :: index_x2o_Foxx_lwup ! longwave radiation (up) (W/m2 ) integer :: index_x2o_Faxa_lwdn ! longwave radiation (down) (W/m2 ) integer :: index_x2o_Fioi_melth ! heat flux from snow & ice melt (W/m2 ) @@ -93,7 +90,7 @@ module mpaso_cpl_indices integer :: index_x2o_Fioi_bergw ! iceberg melt flux (kg/m2/s) integer :: index_x2o_Fioi_salt ! salt (kg(salt)/m2/s) integer :: index_x2o_Foxx_evap ! evaporation flux (kg/m2/s) - integer :: index_x2o_Faxa_prec + integer :: index_x2o_Faxa_prec integer :: index_x2o_Faxa_snow ! water flux due to snow (kg/m2/s) integer :: index_x2o_Faxa_rain ! water flux due to rain (kg/m2/s) integer :: index_x2o_Faxa_bcphidry ! flux: Black Carbon hydrophilic dry deposition @@ -123,17 +120,17 @@ module mpaso_cpl_indices integer :: index_x2o_Foxx_rofPN ! PN runoff flux (kg/m2/s) integer :: index_x2o_Foxx_rofDIC ! DIC runoff flux (kg/m2/s) integer :: index_x2o_Foxx_rofFe ! Fe runoff flux (kg/m2/s) - integer :: index_x2o_Sw_ustokes_wavenumber_1 ! partitioned Stokes drift wavenumber 1 + integer :: index_x2o_Sw_ustokes_wavenumber_1 ! partitioned Stokes drift wavenumber 1 integer :: index_x2o_Sw_vstokes_wavenumber_1 ! partitioned Stokes drift wavenumber 1 - integer :: index_x2o_Sw_ustokes_wavenumber_2 ! partitioned Stokes drift wavenumber 2 + integer :: index_x2o_Sw_ustokes_wavenumber_2 ! partitioned Stokes drift wavenumber 2 integer :: index_x2o_Sw_vstokes_wavenumber_2 ! partitioned Stokes drift wavenumber 2 - integer :: index_x2o_Sw_ustokes_wavenumber_3 ! partitioned Stokes drift wavenumber 3 + integer :: index_x2o_Sw_ustokes_wavenumber_3 ! partitioned Stokes drift wavenumber 3 integer :: index_x2o_Sw_vstokes_wavenumber_3 ! partitioned Stokes drift wavenumber 3 - integer :: index_x2o_Sw_ustokes_wavenumber_4 ! partitioned Stokes drift wavenumber 4 + integer :: index_x2o_Sw_ustokes_wavenumber_4 ! partitioned Stokes drift wavenumber 4 integer :: index_x2o_Sw_vstokes_wavenumber_4 ! partitioned Stokes drift wavenumber 4 - integer :: index_x2o_Sw_ustokes_wavenumber_5 ! partitioned Stokes drift wavenumber 5 + integer :: index_x2o_Sw_ustokes_wavenumber_5 ! partitioned Stokes drift wavenumber 5 integer :: index_x2o_Sw_vstokes_wavenumber_5 ! partitioned Stokes drift wavenumber 5 - integer :: index_x2o_Sw_ustokes_wavenumber_6 ! partitioned Stokes drift wavenumber 6 + integer :: index_x2o_Sw_ustokes_wavenumber_6 ! partitioned Stokes drift wavenumber 6 integer :: index_x2o_Sw_vstokes_wavenumber_6 ! partitioned Stokes drift wavenumber 6 integer :: index_x2o_Sw_Hs ! Significant wave height integer :: index_x2o_Sw_Fp ! Peak wave frequency @@ -150,10 +147,7 @@ module mpaso_cpl_indices ! drv -> glc and internal drv fields - integer :: index_x2o_Fogx_qicelo - integer :: index_x2o_Fogx_qiceho - integer :: index_x2o_Sg_blit - integer :: index_x2o_Sg_blis + integer :: index_x2o_Flgg_hflx integer :: index_x2o_Sg_lithop integer :: index_x2o_Sg_icemask integer :: index_x2o_Sg_icemask_grounded @@ -187,7 +181,7 @@ module mpaso_cpl_indices subroutine mpaso_cpl_indices_set( ) - use seq_flds_mod, only : wav_ocn_coup + use seq_flds_mod, only : wav_ocn_coup use glc_zocnclass_mod type(mct_aVect) :: o2x ! temporary @@ -231,11 +225,8 @@ subroutine mpaso_cpl_indices_set( ) index_o2x_Foxo_ismh = mct_avect_indexra(o2x,'Foxo_ismh',perrWith='quiet') index_o2x_Foxo_rrofih = mct_avect_indexra(o2x,'Foxo_rrofih',perrWith='quiet') - index_o2x_So_blt = mct_avect_indexra(o2x,'So_blt') - index_o2x_So_bls = mct_avect_indexra(o2x,'So_bls') - index_o2x_So_htv = mct_avect_indexra(o2x,'So_htv') - index_o2x_So_stv = mct_avect_indexra(o2x,'So_stv') - index_o2x_So_rhoeff = mct_avect_indexra(o2x,'So_rhoeff') + index_o2x_So_tfrz_isf = mct_avect_indexra(o2x,'So_tfrz_isf') + index_o2x_So_liflfrac = mct_avect_indexra(o2x,'So_liflfrac') if (glc_nzoc > 0) then do iLev = 1, glc_nzoc cnum = glc_zocnclass_as_string(iLev) @@ -283,14 +274,14 @@ subroutine mpaso_cpl_indices_set( ) index_x2o_Foxx_sen = mct_avect_indexra(x2o,'Foxx_sen') index_x2o_Foxx_lwup = mct_avect_indexra(x2o,'Foxx_lwup') index_x2o_Faxa_lwdn = mct_avect_indexra(x2o,'Faxa_lwdn') - index_x2o_Fioi_melth = mct_avect_indexra(x2o,'Fioi_melth') + index_x2o_Fioi_melth = mct_avect_indexra(x2o,'Fioi_melth') index_x2o_Fioi_meltw = mct_avect_indexra(x2o,'Fioi_meltw') index_x2o_Fioi_bergh = mct_avect_indexra(x2o,'PFioi_bergh') index_x2o_Fioi_bergw = mct_avect_indexra(x2o,'PFioi_bergw') - index_x2o_Fioi_salt = mct_avect_indexra(x2o,'Fioi_salt') - index_x2o_Faxa_prec = mct_avect_indexra(x2o,'Faxa_prec') - index_x2o_Faxa_snow = mct_avect_indexra(x2o,'Faxa_snow') - index_x2o_Faxa_rain = mct_avect_indexra(x2o,'Faxa_rain') + index_x2o_Fioi_salt = mct_avect_indexra(x2o,'Fioi_salt') + index_x2o_Faxa_prec = mct_avect_indexra(x2o,'Faxa_prec') + index_x2o_Faxa_snow = mct_avect_indexra(x2o,'Faxa_snow') + index_x2o_Faxa_rain = mct_avect_indexra(x2o,'Faxa_rain') index_x2o_Foxx_evap = mct_avect_indexra(x2o,'Foxx_evap') index_x2o_Foxx_rofl = mct_avect_indexra(x2o,'Foxx_rofl') index_x2o_Foxx_rofi = mct_avect_indexra(x2o,'Foxx_rofi') @@ -345,14 +336,11 @@ subroutine mpaso_cpl_indices_set( ) index_x2o_Fioi_fed2 = mct_avect_indexra(x2o,'Fioi_fed2',perrWith='quiet') index_x2o_Fioi_dust1 = mct_avect_indexra(x2o,'Fioi_dust1',perrWith='quiet') - index_x2o_Fogx_qicelo = mct_avect_indexra(x2o,'Fogx_qicelo') - index_x2o_Fogx_qiceho = mct_avect_indexra(x2o,'Fogx_qiceho') - index_x2o_Sg_blit = mct_avect_indexra(x2o,'Sg_blit') - index_x2o_Sg_blis = mct_avect_indexra(x2o,'Sg_blis') + index_x2o_Flgg_hflx = mct_avect_indexra(x2o,'Flgg_hflx') index_x2o_Sg_lithop = mct_avect_indexra(x2o,'Sg_lithop') index_x2o_Sg_icemask = mct_avect_indexra(x2o,'Sg_icemask') - index_x2o_Sg_icemask = mct_avect_indexra(x2o,'Sg_icemask_grounded') - index_x2o_Sg_icemask = mct_avect_indexra(x2o,'Sg_icemask_floating') + index_x2o_Sg_icemask_grounded = mct_avect_indexra(x2o,'Sg_icemask_grounded') + index_x2o_Sg_icemask_floating = mct_avect_indexra(x2o,'Sg_icemask_floating') if (wav_ocn_coup == 'twoway') then index_x2o_Sw_ustokes_wavenumber_1 = mct_avect_indexra(x2o,'Sw_ustokes_wavenumber_1') diff --git a/components/mpas-ocean/driver/ocn_comp_mct.F b/components/mpas-ocean/driver/ocn_comp_mct.F index c07ad9714542..6579b621dfc3 100644 --- a/components/mpas-ocean/driver/ocn_comp_mct.F +++ b/components/mpas-ocean/driver/ocn_comp_mct.F @@ -125,7 +125,7 @@ module ocn_comp_mct integer :: itimestep, & ! time step number for MPAS ocn_cpl_dt ! length of coupling interval in seconds - set by coupler/ESMF - real (kind=RKIND) :: precFactor + real (kind=RKIND) :: precFactor !======================================================================= @@ -241,10 +241,10 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename )!{{{ real (kind=RKIND), pointer :: & runningMeanRemovedIceRunoff ! the area integrated, running mean of removed ice runoff from the ocean - ! freshwater balance adjustment + ! freshwater balance adjustment logical, pointer :: config_use_precip_scaling character (len=StrKIND), pointer :: config_precip_scaling_mode - real (kind=RKIND), pointer :: SFWFScalingFactor ! scaling factor on freshwater fluxes + real (kind=RKIND), pointer :: SFWFScalingFactor ! scaling factor on freshwater fluxes real(kind=RKIND), pointer :: scalingFactor #ifdef HAVE_MOAB @@ -848,9 +848,9 @@ end subroutine xml_stream_get_attributes ! initialize scaled data ice-shelf melt fluxes based on remove ice runoff call ocn_init_scaled_dismf(domain) - ! Initialize the precipitation scaling scheme + ! Initialize the precipitation scaling scheme call ocn_scaled_sfwf_init(domain) - ! Ensure a valid first guess at the initial time + ! Ensure a valid first guess at the initial time call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling) if (config_use_precip_scaling) then call mpas_pool_get_config(domain % configs, 'config_precip_scaling_mode', config_precip_scaling_mode) @@ -858,10 +858,10 @@ end subroutine xml_stream_get_attributes block_ptr => domain % blocklist call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) call mpas_pool_get_array(forcingPool, "SFWFScalingFactor", SFWFScalingFactor) - if (trim(config_precip_scaling_mode) == 'time-dependent') then + if (trim(config_precip_scaling_mode) == 'time-dependent') then call mpas_pool_get_config(domain % configs, 'config_precip_scaling_initial_factor', scalingFactor) else if (trim(config_precip_scaling_mode) == 'constant') then - call mpas_pool_get_config(domain % configs, 'config_precip_scaling_constant_factor', scalingFactor) + call mpas_pool_get_config(domain % configs, 'config_precip_scaling_constant_factor', scalingFactor) end if ! Broadcast the scalar to the field SFWFScalingFactor = scalingFactor @@ -925,24 +925,14 @@ end subroutine xml_stream_get_attributes ! initialize scaled data ice-shelf melt fluxes based on remove ice runoff call ocn_init_scaled_dismf(domain) - ! initialize scaled data freshwater fluxes based on water balance relationship + ! initialize scaled data freshwater fluxes based on water balance relationship call ocn_scaled_sfwf_init(domain) end if call t_stopf ('mpaso_mct_init') - call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode) - if ( trim(config_land_ice_flux_mode) == 'off' .or. & - trim(config_land_ice_flux_mode) == 'pressure_only' .or. & - trim(config_land_ice_flux_mode) == 'data' .or. & - trim(config_land_ice_flux_mode) == 'standalone' ) then - ocn_c2_glcshelf = .false. - else if ( trim(config_land_ice_flux_mode) .eq. 'coupled' ) then - ocn_c2_glcshelf = .true. - else - call mpas_log_write('ERROR: unknown land_ice_flux_mode: ' // trim(config_land_ice_flux_mode), MPAS_LOG_CRIT) - end if + ocn_c2_glcshelf = .false. call seq_infodata_PutData(infodata, ocn_prognostic=.true., ocnrof_prognostic=.true., & ocn_c2_glcshelf=ocn_c2_glcshelf) @@ -964,13 +954,13 @@ end subroutine xml_stream_get_attributes !configuration for freshwater conservation call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling) - if (config_use_precip_scaling) then + if (config_use_precip_scaling) then ! independent of space so should be no need to loop over blocks block_ptr => domain % blocklist call mpas_pool_get_subpool(block_ptr % structs, 'forcing', forcingPool) call mpas_pool_get_array(forcingPool, "SFWFScalingFactor", SFWFScalingFactor) call seq_infodata_PutData(infodata, precip_fact=SFWFScalingFactor) - end if + end if !----------------------------------------------------------------------- ! @@ -1084,9 +1074,9 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ real (kind=RKIND), pointer :: & runningMeanRemovedIceRunoff ! the area integrated, running mean of removed ice runoff from the ocean - ! freshwater balance adjustment + ! freshwater balance adjustment logical, pointer :: config_use_precip_scaling - real (kind=RKIND), pointer :: SFWFScalingFactor ! scaling factor on freshwater fluxes + real (kind=RKIND), pointer :: SFWFScalingFactor ! scaling factor on freshwater fluxes #ifdef HAVE_MOAB integer :: ent_type @@ -1328,7 +1318,7 @@ subroutine ocn_run_mct( EClock, cdata_o, x2o_o, o2x_o)!{{{ call mpas_log_write(' Validating ocean state') endif - ! update scalling factors for freshwater fluxes + ! update scalling factors for freshwater fluxes call ocn_update_scaling_factor(domain, timeLevel=1) call ocn_validate_state(domain, timeLevel=1) @@ -2105,8 +2095,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ riverFluxDICField, & riverFluxALKField, & riverFluxFeField, & + landIceConductiveHeatFluxField, & landIceFreshwaterFluxField, & - landIceHeatFluxField, & landIceFractionField, & windSpeed10mField, & significantWaveHeightField, & @@ -2165,8 +2155,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ riverFluxDIC, & riverFluxALK, & riverFluxFe, & + landIceConductiveHeatFlux, & landIceFreshwaterFlux, & - landIceHeatFlux, & landIceFraction, & areaCell, & windSpeed10m, & @@ -2198,7 +2188,7 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ real (kind=RKIND) :: riverFactor - ! freshwater balance adjustment + ! freshwater balance adjustment logical, pointer :: config_use_precip_scaling !----------------------------------------------------------------------- @@ -2237,7 +2227,7 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ !configuration for freshwater balance constrains call mpas_pool_get_config(domain % configs, 'config_use_precip_scaling', config_use_precip_scaling) if (config_use_precip_scaling) then - !extract scaling factors + !extract scaling factors call seq_infodata_GetData(infodata, precip_fact=precFactor) end if @@ -2290,8 +2280,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ call mpas_pool_get_field(forcingPool, 'waveNetOceanZonalWindStress', waveNetOceanZonalWindStressField) call mpas_pool_get_field(forcingPool, 'waveNetOceanMeridionalWindStress', waveNetOceanMeridionalWindStressField) + call mpas_pool_get_field(forcingPool, 'landIceConductiveHeatFlux', landIceConductiveHeatFluxField) call mpas_pool_get_field(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFluxField) - call mpas_pool_get_field(forcingPool, 'landIceHeatFlux', landIceHeatFluxField) call mpas_pool_get_field(forcingPool, 'landIceFraction', landIceFractionField) call mpas_pool_get_field(forcingPool, 'landIceInterfaceTracers', landIceInterfaceTracersField) @@ -2325,8 +2315,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ iceRunoffFlux => iceRunoffFluxField % array removedRiverRunoffFlux => removedRiverRunoffFluxField % array removedIceRunoffFlux => removedIceRunoffFluxField % array + landIceConductiveHeatFlux => landIceConductiveHeatFluxField % array landIceFreshwaterFlux => landIceFreshwaterFluxField % array - landIceHeatFlux => landIceHeatFluxField % array landIceInterfaceTracers => landIceInterfaceTracersField % array landIceFraction => landIceFractionField % array windSpeed10m => windSpeed10mField % array @@ -2451,14 +2441,14 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ do i = 1, nCellsSolve n = n + 1 if ( windStressZonalField % isActive ) then - if (config_momentum_use_active_wave) then + if (config_momentum_use_active_wave) then windStressZonal(i) = x2o_o(index_x2o_Foxx_taux, n) - (x2o_o(index_x2o_Faww_Tawx, n) - x2o_o(index_x2o_Fwow_Twox, n)) else windStressZonal(i) = x2o_o(index_x2o_Foxx_taux, n) endif end if if ( windStressMeridionalField % isActive ) then - if (config_momentum_use_active_wave) then + if (config_momentum_use_active_wave) then windStressMeridional(i) = x2o_o(index_x2o_Foxx_tauy, n) - (x2o_o(index_x2o_Faww_Tawy, n) - x2o_o(index_x2o_Fwow_Twoy, n)) else windStressMeridional(i) = x2o_o(index_x2o_Foxx_tauy, n) @@ -2597,15 +2587,8 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ end if endif - if ( landIceFreshwaterFluxField % isActive ) then - !landIceFreshwaterFlux(i) = x2o_o(index_x2o_Fogx_qicelo, n) - end if - if ( landIceHeatFluxField % isActive ) then - !landIceHeatFlux(i) = x2o_o(index_x2o_Fogx_qiceho, n) - end if - if ( landIceInterfaceTracersField % isActive ) then - !landIceInterfaceTracers(indexIT, i) = x2o_o(index_x2o_Sg_blit, n) - !landIceInterfaceTracers(indexIS, i) = x2o_o(index_x2o_Sg_blis, n) + if ( landIceConductiveHeatFluxField % isActive ) then + landIceConductiveHeatFlux(i) = x2o_o(index_x2o_Flgg_hflx, n) end if if ( landIceFractionField % isActive ) then !landIceFraction(i) = x2o_o(index_x2o_Sg_icemask, n) @@ -2778,7 +2761,6 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ call mpas_pool_get_field(forcingPool, 'waveNetOceanMeridionalWindStress', waveNetOceanMeridionalWindStressField) call mpas_pool_get_field(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFluxField) - call mpas_pool_get_field(forcingPool, 'landIceHeatFlux', landIceHeatFluxField) call mpas_pool_get_field(forcingPool, 'landIceFraction', landIceFractionField) call mpas_pool_get_field(forcingPool, 'landIceInterfaceTracers', landIceInterfaceTracersField) @@ -2940,12 +2922,12 @@ subroutine ocn_import_mct(x2o_o, errorCode)!{{{ call mpas_dmpar_exch_halo_field(waveNetOceanMeridionalWindStressField) end if + if ( landIceConductiveHeatFluxField % isActive ) then + call mpas_dmpar_exch_halo_field(landIceConductiveHeatFluxField) + end if if ( landIceFreshwaterFluxField % isActive ) then call mpas_dmpar_exch_halo_field(landIceFreshwaterFluxField) end if - if ( landIceHeatFluxField % isActive ) then - call mpas_dmpar_exch_halo_field(landIceHeatFluxField) - end if if ( landIceInterfaceTracersField % isActive ) then call mpas_dmpar_exch_halo_field(landIceInterfaceTracersField) end if @@ -3110,6 +3092,7 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ integer, dimension(:), pointer :: landIceMask real (kind=RKIND), dimension(:), pointer :: seaIceEnergy, accumulatedFrazilIceMass, frazilSurfacePressure, & + landIceFloatingFraction, & avgTotalFreshWaterTemperatureFlux, & avgCO2_gas_flux, DMSFlux, surfaceUpwardCO2Flux, & avgOceanSurfaceDIC, & @@ -3127,16 +3110,17 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ avgLandIceFreshwaterFlux, & avgRemovedRiverRunoffFlux, & avgRemovedIceRunoffFlux, & - avgLandIceHeatFlux, & + avgLandIceCouplingHeatFlux, & avgRemovedIceRunoffHeatFlux real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, & avgSSHGradient, avgOceanSurfacePhytoC, & avgOceanSurfaceDOC, layerThickness, & avgThermalForcingAtZLevels, & - avgThermalForcingAtZLevelsMask + avgThermalForcingAtZLevelsMask, & + avgLandIceInterfaceTracers - real (kind=RKIND) :: surfaceFreezingTemp + real (kind=RKIND) :: surfaceFreezingTemp, landIceInterfaceTemperatureC logical, pointer :: frazilIceActive, & config_remove_ais_river_runoff, & @@ -3221,6 +3205,7 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ call mpas_pool_get_array(statePool, 'layerThickness', layerThickness, 1) call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask) + call mpas_pool_get_array(forcingPool, 'landIceFloatingFraction', landIceFloatingFraction) call mpas_pool_get_array(forcingPool, 'avgTracersSurfaceValue', avgTracersSurfaceValue) call mpas_pool_get_array(forcingPool, 'avgSurfaceVelocity', avgSurfaceVelocity) call mpas_pool_get_array(forcingPool, 'avgSSHGradient', avgSSHGradient) @@ -3233,9 +3218,9 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ end if ! Cryo fields - if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then + if (trim(config_land_ice_flux_mode) == 'active' .or. trim(config_land_ice_flux_mode) == 'data') then call mpas_pool_get_array(forcingPool, 'avgLandIceFreshwaterFlux', avgLandIceFreshwaterFlux) - call mpas_pool_get_array(forcingPool, 'avgLandIceHeatFlux', avgLandIceHeatFlux) + call mpas_pool_get_array(forcingPool, 'avgLandIceCouplingHeatFlux', avgLandIceCouplingHeatFlux) endif if (config_remove_ais_river_runoff) then call mpas_pool_get_array(forcingPool, 'avgRemovedRiverRunoffFlux', avgRemovedRiverRunoffFlux) @@ -3283,6 +3268,10 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ call mpas_pool_get_array(MacroMoleculesSeaIceCoupling, 'avgOceanSurfaceDOC', avgOceanSurfaceDOC) call mpas_pool_get_array(MacroMoleculesSeaIceCoupling, 'avgOceanSurfaceDON', avgOceanSurfaceDON) endif + if(trim(config_land_ice_flux_mode) == 'active') then + call mpas_pool_get_array(forcingPool, 'avgLandIceInterfaceTracers', avgLandIceInterfaceTracers) + end if + ! call mpas_pool_get_array(forcingPool, 'CO2Flux', CO2Flux) ! call mpas_pool_get_array(forcingPool, 'DMSFlux', DMSFlux) ! call mpas_pool_get_array(forcingPool, 'surfaceUpwardCO2Flux', surfaceUpwardCO2Flux) @@ -3302,9 +3291,9 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ o2x_o(index_o2x_Faoo_h2otemp, n) = avgTotalFreshWaterTemperatureFlux(i) * rho_sw * cp_sw ! Cryo fields - if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then - o2x_o(index_o2x_Foxo_ismw, n) = avgLandIceFreshwaterFlux(i) - o2x_o(index_o2x_Foxo_ismh, n) = avgLandIceHeatFlux(i) + if (trim(config_land_ice_flux_mode) == 'active' .or. trim(config_land_ice_flux_mode) == 'data') then + o2x_o % rAttr(index_o2x_Foxo_ismw, n) = avgLandIceFreshwaterFlux(i) + o2x_o % rAttr(index_o2x_Foxo_ismh, n) = avgLandIceCouplingHeatFlux(i) endif if (config_remove_ais_river_runoff) then o2x_o(index_o2x_Foxo_rrofl, n) = avgRemovedRiverRunoffFlux(i) @@ -3352,11 +3341,19 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ else - o2x_o(index_o2x_Fioo_q, n) = 0.0_RKIND - o2x_o(index_o2x_Fioo_frazil, n) = 0.0_RKIND - if (trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then - o2x_o(index_o2x_Foxo_q_li, n) = accumulatedFrazilIceMass(i) * config_frazil_heat_of_fusion / ocn_cpl_dt - o2x_o(index_o2x_Foxo_frazil_li, n) = accumulatedFrazilIceMass(i) / ocn_cpl_dt + o2x_o % rAttr(index_o2x_Fioo_q, n) = 0.0_RKIND + o2x_o % rAttr(index_o2x_Fioo_frazil, n) = 0.0_RKIND + if (trim(config_land_ice_flux_mode) == 'active' .or. trim(config_land_ice_flux_mode) == 'data') then + if (trim(config_land_ice_flux_mode) == 'active') then + ! Export the frazil enthalpy carried with the frozen freshwater transferred to MALI. + ! The latent heat associated with frazil formation stays in MPAS-Ocean. + landIceInterfaceTemperatureC = avgLandIceInterfaceTracers(indexIT, i) - T0_Kelvin + o2x_o % rAttr(index_o2x_Foxo_q_li, n) = -accumulatedFrazilIceMass(i) * cp_sw & + * landIceInterfaceTemperatureC / ocn_cpl_dt + else + o2x_o % rAttr(index_o2x_Foxo_q_li, n) = 0.0_RKIND + endif + o2x_o % rAttr(index_o2x_Foxo_frazil_li, n) = accumulatedFrazilIceMass(i) / ocn_cpl_dt endif end if @@ -3402,19 +3399,12 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ ! o2x_o(index_o2x_Faoo_fdms_ocn, n) = DMSFlux(i) ! o2x_o(index_o2x_Faoo_fco2_ocn, n) = surfaceUpwardCO2Flux(i) -!JW o2x_o(index_o2x_So_blt, n) = landIceBoundaryLayerTemperature(i) -!JW o2x_o(index_o2x_So_bls, n) = landIceBoundaryLayerSalinity(i) -!JW o2x_o(index_o2x_So_htv, n) = landIceHeatTransferVelocity(i) -!JW o2x_o(index_o2x_So_stv, n) = landIceSaltTransferVelocity(i) - - if ( trim(config_land_ice_flux_mode) .eq. 'standalone' .or. & - trim(config_land_ice_flux_mode) .eq. 'coupled' ) then - o2x_o(index_o2x_So_blt, n) = landIceBoundaryLayerTracers(indexBLT,i) - o2x_o(index_o2x_So_bls, n) = landIceBoundaryLayerTracers(indexBLS,i) - o2x_o(index_o2x_So_htv, n) = landIceTracerTransferVelocities(indexHeatTrans,i) - o2x_o(index_o2x_So_stv, n) = landIceTracerTransferVelocities(indexSaltTrans,i) - o2x_o(index_o2x_So_rhoeff, n) = 0.0_RKIND + if ( trim(config_land_ice_flux_mode) .eq. 'active' ) then + o2x_o % rAttr(index_o2x_So_tfrz_isf, n) = landIceFloatingFraction(i) & + * avgLandIceInterfaceTracers(indexIT,i) + o2x_o % rAttr(index_o2x_So_liflfrac, n) = landIceFloatingFraction(i) endif + if (trim(config_glc_thermal_forcing_coupling_mode) == '3d') then do iLevel = 1, nGlcZLevels if (avgThermalForcingAtZLevelsMask(ilevel, i) > 0.9_RKIND) then @@ -3432,20 +3422,6 @@ subroutine ocn_export_mct(o2x_o, errorCode) !{{{ enddo endif - - !Fyke: test - !write(stderrUnit,*) 'n=',n - !write(stderrUnit,*) 'o2x_o(index_o2x_So_blt, n)=',o2x_o(index_o2x_So_blt, n) - !write(stderrUnit,*) 'o2x_o(index_o2x_So_bls, n)=',o2x_o(index_o2x_So_bls, n) - !write(stderrUnit,*) 'o2x_o(index_o2x_So_htv, n)=',o2x_o(index_o2x_So_htv, n) - !write(stderrUnit,*) 'o2x_o(index_o2x_So_stv, n)=',o2x_o(index_o2x_So_stv, n) - !write(stderrUnit,*) 'o2x_o(index_o2x_So_rhoeff, n)=',o2x_o(index_o2x_So_rhoeff, n) - !o2x_o(index_o2x_So_blt, n) = 0._r8 - !o2x_o(index_o2x_So_bls, n) = 34.5_r8 - !o2x_o(index_o2x_So_htv, n) = 1.e-4_r8 - !o2x_o(index_o2x_So_stv, n) = 3.e-6_r8 - !o2x_o(index_o2x_So_rhoeff, n) = 1000._r8*9.81_r8*918._r8 !lithostatic pressure of 1km of ice - end do block_ptr => block_ptr % next diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml index 1adc60ceb6f2..72c5c0c5023c 100644 --- a/components/mpas-ocean/src/Registry.xml +++ b/components/mpas-ocean/src/Registry.xml @@ -1099,19 +1099,23 @@ /> + - @@ -1904,7 +1907,6 @@ - @@ -2017,7 +2019,6 @@ - @@ -2136,7 +2137,6 @@ - @@ -2168,7 +2168,7 @@ - + - - - + @@ -4214,36 +4212,32 @@ description="Total flux of mass through the ocean surface that includes interface melt freeze and column averaged frazil freeze rate. Positive into ocean." packages="landIceFluxesPKG;dataLandIceFluxesPKG" /> + + - - - + - - - - - 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 95a9bca8f70c..bafc437a0a10 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 @@ -597,6 +597,7 @@ subroutine energy_conservation(domain, err) if (config_use_frazil_ice_formation) then do iCell = 1, nCellsSolve + ! Frazil latent heat remains internal to MPAS-Ocean. ! Frazil ice mass is negative. Negative coefficient makes heat ! flux positive, because freezing ice releases heat. sumArray(10) = sumArray(10) - areaCell(iCell) * config_frazil_heat_of_fusion & @@ -606,6 +607,7 @@ subroutine energy_conservation(domain, err) if (landIceFreshwaterFluxesOn) then do iCell = 1, nCellsSolve + ! Export only the non-latent interfacial shelf heat exchange here. sumArray(17) = sumArray(17) + areaCell(iCell) * landIceHeatFlux(iCell) enddo end if @@ -614,6 +616,7 @@ subroutine energy_conservation(domain, err) call mpas_pool_get_array(statePool, 'accumulatedLandIceFrazilMass', accumulatedLandIceFrazilMassOld, 1) call mpas_pool_get_array(statePool, 'accumulatedLandIceFrazilMass', accumulatedLandIceFrazilMassNew, 2) do iCell = 1, nCellsSolve + ! Under-land-ice frazil latent heat also remains internal to MPAS-Ocean. ! Frazil ice mass is negative. Negative coefficient makes heat ! flux positive, because freezing ice releases heat. sumArray(18) = sumArray(18) - areaCell(iCell) * config_frazil_heat_of_fusion & @@ -765,10 +768,10 @@ subroutine energy_conservation(domain, err) v=accumulatedSensibleHeatFlux ; write(m,"('sensibleHeatFlux ',es16.8,' x2o_Foxx_sen hsen ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v v=accumulatedIcebergHeatFlux ; write(m,"('icebergHeatFlux ',es16.8,' x2o_Fioi_bergh hberg ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v if (landIceFreshwaterFluxesOn) then -v=accumulatedLandIceHeatFlux ; write(m,"('landIceHeatFlux ',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v +v=accumulatedLandIceHeatFlux ; write(m,"('landIceHeatFlux (non-lat)',es16.8,' ',f16.8)") v,v/A; call mpas_log_write(m); s=s+v end if if (config_use_frazil_ice_formation .and. config_frazil_under_land_ice) then -v=accumulatedLandIceFrazilHeatFlux ; write(m,"(' landIceFrazilHeatFlux ',es16.8,' (already in hfreeze, do not sum )',f16.8)") v,v/A; call mpas_log_write(m); ! no sum: s=s+v +v=accumulatedLandIceFrazilHeatFlux ; write(m,"(' landIceFrazilHeatFlux ',es16.8,' (latent, ocean-internal; do not sum )',f16.8)") v,v/A; call mpas_log_write(m); ! no sum: s=s+v end if write(m,"('SUM EXPLICIT HEAT FLUXES ',es16.8,' ',f16.8)") s, s/A; call mpas_log_write(m) explicitHeatFluxSum = s @@ -2503,8 +2506,7 @@ subroutine ocn_restart_conservation_check(domain, err)!{{{ err = 0 - if ( trim(config_land_ice_flux_mode) == 'standalone' .or. & - trim(config_land_ice_flux_mode) == 'coupled' .or. & + if ( trim(config_land_ice_flux_mode) == 'active' .or. & trim(config_land_ice_flux_mode) == 'data' ) then landIceFreshwaterFluxesOn = .true. end if @@ -2563,4 +2565,3 @@ end subroutine ocn_finalize_conservation_check!}}} end module ocn_conservation_check ! vim: foldmethod=marker - 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 4b249153853b..750596f71557 100644 --- a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F +++ b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F @@ -122,7 +122,6 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: landIcePressurePKGActive logical, pointer :: dataLandIceFluxesPKGActive logical, pointer :: landIceFluxesPKGActive - logical, pointer :: landIceCouplingPKGActive logical, pointer :: landIceThermalForcing3DPKGActive logical, pointer :: dataSubglacialRunoffFluxPKGActive logical, pointer :: thicknessBulkPKGActive @@ -193,6 +192,7 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ logical, pointer :: config_use_bulk_wind_stress logical, pointer :: config_use_bulk_thickness_flux logical, pointer :: config_compute_active_tracer_budgets + logical, pointer :: config_use_land_ice_pressure character (len=StrKIND), pointer :: config_land_ice_flux_mode character (len=StrKIND), pointer :: config_glc_thermal_forcing_coupling_mode character (len=StrKIND), pointer :: config_subglacial_runoff_mode @@ -307,26 +307,20 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{ ! ! test for land ice pressure, landIcePressurePKG ! test for land ice fluxes, landIceFluxesPKG - ! test for land ice coupling, landIceCouplingPKG ! call mpas_pool_get_package(packagePool, 'landIcePressurePKGActive', landIcePressurePKGActive) call mpas_pool_get_package(packagePool, 'dataLandIceFluxesPKGActive', & dataLandIceFluxesPKGActive) call mpas_pool_get_package(packagePool, 'landIceFluxesPKGActive', landIceFluxesPKGActive) - call mpas_pool_get_package(packagePool, 'landIceCouplingPKGActive', landIceCouplingPKGActive) + call mpas_pool_get_config(configPool, 'config_use_land_ice_pressure', config_use_land_ice_pressure) call mpas_pool_get_config(configPool, 'config_land_ice_flux_mode', config_land_ice_flux_mode) - if ( trim(config_land_ice_flux_mode) == 'pressure_only' ) then - landIcePressurePKGActive = .true. - else if ( trim(config_land_ice_flux_mode) == 'data' ) then + if ( config_use_land_ice_pressure ) then landIcePressurePKGActive = .true. + end if + if ( trim(config_land_ice_flux_mode) == 'data' ) then dataLandIceFluxesPKGActive = .true. - else if ( trim(config_land_ice_flux_mode) == 'standalone' ) then - landIcePressurePKGActive = .true. - landIceFluxesPKGActive = .true. - else if ( trim(config_land_ice_flux_mode) == 'coupled' ) then - landIcePressurePKGActive = .true. + else if ( trim(config_land_ice_flux_mode) == 'active' ) then landIceFluxesPKGActive = .true. - landIceCouplingPKGActive = .true. end if diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F index 9c4cda353dff..263bdafd9960 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F @@ -40,7 +40,6 @@ module ocn_time_integration_rk4 use ocn_time_average_coupled use ocn_wetting_drying - use ocn_effective_density_in_land_ice use ocn_surface_land_ice_fluxes use ocn_transport_tests use ocn_time_varying_forcing @@ -143,7 +142,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ logical, pointer :: config_verify_not_dry logical, pointer :: config_prevent_drying logical, pointer :: config_zero_drying_velocity - character (len=StrKIND), pointer :: config_land_ice_flux_mode + logical, pointer :: config_use_land_ice_pressure + character (len=StrKIND), pointer :: config_land_ice_draft_mode ! Diagnostics Indices @@ -172,7 +172,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ real (kind=RKIND), dimension(:,:,:), pointer :: tracerGroup, tracersCur, tracersNew ! Diagnostics Field Pointers - type (field1DReal), pointer :: boundaryLayerDepthField, effectiveDensityField + type (field1DReal), pointer :: boundaryLayerDepthField type (field2DReal), pointer :: normalizedRelativeVorticityEdgeField, divergenceField, relativeVorticityField ! State/Tend Field Pointers @@ -211,7 +211,8 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ call mpas_pool_get_config(domain % configs, 'config_use_freq_filtered_thickness', config_use_freq_filtered_thickness) call mpas_pool_get_config(domain % configs, 'config_use_GM', config_use_GM) call mpas_pool_get_config(domain % configs, 'config_use_cvmix_kpp', config_use_cvmix_kpp) - call mpas_pool_get_config(domain % configs, 'config_land_ice_flux_mode', config_land_ice_flux_mode) + call mpas_pool_get_config(domain % configs, 'config_use_land_ice_pressure', config_use_land_ice_pressure) + call mpas_pool_get_config(domain % configs, 'config_land_ice_draft_mode', config_land_ice_draft_mode) call mpas_pool_get_config(domain % configs, 'config_disable_vel_all_tend', config_disable_vel_all_tend) call mpas_pool_get_config(domain % configs, 'config_disable_thick_all_tend', config_disable_thick_all_tend) call mpas_pool_get_config(domain % configs, 'config_disable_tr_all_tend', config_disable_tr_all_tend) @@ -757,10 +758,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ #endif call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, scratchPool, tracersPool, 2) - ! Update the effective density in land ice if we're coupling to land ice - call ocn_effective_density_in_land_ice_update(forcingPool, & - statePool, err) - #ifdef MPAS_OPENACC !$acc update host(layerThickEdgeFlux, layerThickEdgeMean) !$acc update host(relativeVorticity, circulation) @@ -856,14 +853,6 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{ block => block % next end do - if (trim(config_land_ice_flux_mode) == 'coupled') then - call mpas_timer_start("RK4-effective density halo") - call mpas_pool_get_subpool(domain % blocklist % structs, 'state', statePool) - call mpas_pool_get_field(statePool, 'effectiveDensityInLandIce', effectiveDensityField, 2) - call mpas_dmpar_exch_halo_field(effectiveDensityField) - call mpas_timer_stop("RK4-effective density halo") - end if - call mpas_timer_stop("RK4-cleanup phase") end subroutine ocn_time_integrator_rk4!}}} diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F index 9e3923d577a5..2db81d584a5d 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F @@ -51,7 +51,6 @@ module ocn_time_integration_si use ocn_vertical_remap use ocn_time_average_coupled - use ocn_effective_density_in_land_ice use ocn_surface_land_ice_fluxes use ocn_transport_tests @@ -330,9 +329,6 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ lowFreqDivergenceNew, &! low frq div at new time tracerGroupIdealAgeMask ! mask for resetting surface ideal age - type (field1DReal), pointer :: & - effectiveDensityField ! field pointer for halo update - type (dm_info) :: dminfo ! State/Tracer arrays, info @@ -3030,11 +3026,6 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, & scratchPool, tracersPool, 2) - ! Update the effective desnity in land ice if we're coupling to - ! land ice - call ocn_effective_density_in_land_ice_update(forcingPool, & - statePool, err) - call mpas_timer_start('si final mpas reconstruct', .false.) #ifdef MPAS_OPENACC @@ -3141,14 +3132,6 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{ call ocn_reconstruct_eddy_vectors(meshPool) end if - if (trim(config_land_ice_flux_mode) == 'coupled') then - call mpas_timer_start("si effective density halo") - call mpas_pool_get_field(statePool, & - 'effectiveDensityInLandIce', & - effectiveDensityField, 2) - call mpas_dmpar_exch_halo_field(effectiveDensityField) - call mpas_timer_stop("si effective density halo") - end if deallocate(bottomDepthEdge) call mpas_timer_stop('si fini') diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F index 74ec6fa83c8c..7ffb674df356 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F @@ -45,7 +45,6 @@ module ocn_time_integration_split use ocn_vertical_remap use ocn_time_average_coupled - use ocn_effective_density_in_land_ice use ocn_surface_land_ice_fluxes use ocn_transport_tests @@ -229,9 +228,6 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ lowFreqDivergenceNew, &! low frq div at new time tracerGroupIdealAgeMask ! mask for resetting surface ideal age - type (field1DReal), pointer :: & - effectiveDensityField ! field pointer for halo update - ! State/Tracer arrays, info real (kind=RKIND), dimension(:,:,:), pointer :: &! tracersGroupCur, &! old, new tracer arrays @@ -2680,11 +2676,6 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, & scratchPool, tracersPool, 2) - ! Update the effective desnity in land ice if we're coupling to - ! land ice - call ocn_effective_density_in_land_ice_update(forcingPool, & - statePool, err) - call mpas_timer_start('se final mpas reconstruct', .false.) #ifdef MPAS_OPENACC @@ -2791,14 +2782,6 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{ call ocn_reconstruct_eddy_vectors(meshPool) end if - if (trim(config_land_ice_flux_mode) == 'coupled') then - call mpas_timer_start("se effective density halo") - call mpas_pool_get_field(statePool, & - 'effectiveDensityInLandIce', & - effectiveDensityField, 2) - call mpas_dmpar_exch_halo_field(effectiveDensityField) - call mpas_timer_stop("se effective density halo") - end if deallocate(baroclinicThickness) call mpas_timer_stop('se fini') diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F index acd717e75bd3..5657b0bc3cf7 100644 --- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F +++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split_ab2.F @@ -49,7 +49,6 @@ module ocn_time_integration_split_ab2 use ocn_vertical_remap use ocn_time_average_coupled - use ocn_effective_density_in_land_ice use ocn_surface_land_ice_fluxes use ocn_transport_tests @@ -237,9 +236,6 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ lowFreqDivergenceNew, &! low frq div at new time tracerGroupIdealAgeMask ! mask for resetting surface ideal age - type (field1DReal), pointer :: & - effectiveDensityField ! field pointer for halo update - ! State/Tracer arrays, info real (kind=RKIND), dimension(:,:,:), pointer :: &! tracersGroupCur, &! old, new tracer arrays @@ -2828,11 +2824,6 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ call ocn_diagnostic_solve(dt, statePool, forcingPool, meshPool, verticalMeshPool, & scratchPool, tracersPool, 2) - ! Update the effective desnity in land ice if we're coupling to - ! land ice - call ocn_effective_density_in_land_ice_update(forcingPool, & - statePool, err) - call mpas_timer_start('se final mpas reconstruct', .false.) #ifdef MPAS_OPENACC @@ -2939,14 +2930,6 @@ subroutine ocn_time_integrator_split_ab2(domain, dt)!{{{ call ocn_reconstruct_eddy_vectors(meshPool) end if - if (trim(config_land_ice_flux_mode) == 'coupled') then - call mpas_timer_start("se effective density halo") - call mpas_pool_get_field(statePool, & - 'effectiveDensityInLandIce', & - effectiveDensityField, 2) - call mpas_dmpar_exch_halo_field(effectiveDensityField) - call mpas_timer_stop("se effective density halo") - end if deallocate(bottomDepthEdge) ! get current time & setup num TS iteration diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_isomip_plus.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_isomip_plus.F index 78436a13af4b..34b5a940f5fe 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_isomip_plus.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_isomip_plus.F @@ -107,7 +107,7 @@ subroutine ocn_init_setup_isomip_plus(domain, iErr)!{{{ integer, dimension(:), pointer :: minLevelCell, maxLevelCell real(kind=RKIND), dimension(:), pointer :: refBottomDepth, & - fCell, fEdge, fVertex, effectiveDensityInLandIce, xCell, & + fCell, fEdge, fVertex, xCell, & refZMid, refLayerThickness real(kind=RKIND), dimension(:,:), pointer :: layerThickness @@ -301,11 +301,6 @@ subroutine ocn_init_setup_isomip_plus(domain, iErr)!{{{ fEdge(:) = config_isomip_plus_coriolis_parameter fVertex(:) = config_isomip_plus_coriolis_parameter - if(config_land_ice_flux_mode == 'coupled') then - call mpas_pool_get_array(statePool, 'effectiveDensityInLandIce', effectiveDensityInLandIce, 1) - effectiveDensityInLandIce(:) = config_isomip_plus_effective_density - end if - block_ptr => block_ptr % next end do @@ -684,9 +679,9 @@ subroutine ocn_init_validate_isomip_plus(configPool, packagePool, iocontext, iEr call mpas_pool_get_config(configPool, 'config_land_ice_flux_mode', config_land_ice_flux_mode) - if(config_land_ice_flux_mode .ne. 'standalone' .and. config_land_ice_flux_mode .ne. 'coupled') then - call mpas_log_write( 'Validation failed for isomip_plus. config_land_ice_flux_mode must either'// & - ' be standalone or coupled.', MPAS_LOG_CRIT) + if(config_land_ice_flux_mode .ne. 'active') then + call mpas_log_write( 'Validation failed for isomip_plus. config_land_ice_flux_mode must be'// & + ' active.', MPAS_LOG_CRIT) iErr = 1 return end if diff --git a/components/mpas-ocean/src/mode_init/mpas_ocn_init_ssh_and_landIcePressure.F b/components/mpas-ocean/src/mode_init/mpas_ocn_init_ssh_and_landIcePressure.F index 09e1db7333f8..2f31781516ab 100644 --- a/components/mpas-ocean/src/mode_init/mpas_ocn_init_ssh_and_landIcePressure.F +++ b/components/mpas-ocean/src/mode_init/mpas_ocn_init_ssh_and_landIcePressure.F @@ -96,8 +96,7 @@ subroutine ocn_init_ssh_and_landIcePressure_balance(domain, iErr)!{{{ real (kind=RKIND), dimension(:), pointer :: ssh - real (kind=RKIND), dimension(:), pointer :: landIcePressure, landIceDraft, & - effectiveDensityInLandIce + real (kind=RKIND), dimension(:), pointer :: landIcePressure, landIceDraft integer, dimension(:), pointer :: minLevelCell @@ -145,7 +144,6 @@ subroutine ocn_init_ssh_and_landIcePressure_balance(domain, iErr)!{{{ call mpas_pool_get_array(forcingPool, 'landIceDraft', landIceDraft) call mpas_pool_get_array(statePool, 'ssh', ssh, 1) - call mpas_pool_get_array(statePool, 'effectiveDensityInLandIce', effectiveDensityInLandIce, 1) call ocn_equation_of_state_density(statePool, meshPool, tracersSurfaceValue, & nCells, 0, 'relative', density, iErr, & @@ -160,17 +158,10 @@ subroutine ocn_init_ssh_and_landIcePressure_balance(domain, iErr)!{{{ if(sshAdjustmentMask(iCell) == 0) then ssh(iCell) = 0.0_RKIND landIcePressure(iCell) = 0.0_RKIND - - if (associated(effectiveDensityInLandIce)) & - ! effective density cannot be determined - effectiveDensityInLandIce(iCell) = 0.0_RKIND - cycle end if landIcePressure(iCell) = & max(0.0_RKIND, -density(minLevelCell(iCell), iCell)*gravity*ssh(iCell)) - if (associated(effectiveDensityInLandIce)) & - effectiveDensityInLandIce(iCell) = density(minLevelCell(iCell), iCell) end do ! copy the SSH into the landIceDraft so we can use it later to remove it when diff --git a/components/mpas-ocean/src/ocean.cmake b/components/mpas-ocean/src/ocean.cmake index 60a9c0915d51..552ac8b57cc7 100644 --- a/components/mpas-ocean/src/ocean.cmake +++ b/components/mpas-ocean/src/ocean.cmake @@ -103,7 +103,6 @@ list(APPEND RAW_SOURCES core_ocean/shared/mpas_ocn_forcing.F core_ocean/shared/mpas_ocn_surface_bulk_forcing.F core_ocean/shared/mpas_ocn_surface_land_ice_fluxes.F - core_ocean/shared/mpas_ocn_effective_density_in_land_ice.F core_ocean/shared/mpas_ocn_frazil_forcing.F core_ocean/shared/mpas_ocn_tidal_forcing.F core_ocean/shared/mpas_ocn_time_average_coupled.F diff --git a/components/mpas-ocean/src/shared/Makefile b/components/mpas-ocean/src/shared/Makefile index 92019f57efe8..48f68a1e02b7 100644 --- a/components/mpas-ocean/src/shared/Makefile +++ b/components/mpas-ocean/src/shared/Makefile @@ -62,7 +62,6 @@ OBJS = mpas_ocn_init_routines.o \ mpas_ocn_forcing.o \ mpas_ocn_surface_bulk_forcing.o \ mpas_ocn_surface_land_ice_fluxes.o \ - mpas_ocn_effective_density_in_land_ice.o \ mpas_ocn_frazil_forcing.o \ mpas_ocn_tidal_forcing.o \ mpas_ocn_forcing_restoring.o \ @@ -189,8 +188,6 @@ mpas_ocn_tidal_forcing.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_equati mpas_ocn_transport_tests.o: mpas_ocn_config.o -mpas_ocn_effective_density_in_land_ice.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o - mpas_ocn_forcing_restoring.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_mesh.o mpas_ocn_tracer_surface_restoring.o: mpas_ocn_constants.o mpas_ocn_config.o mpas_ocn_framework_forcing.o mpas_ocn_diagnostics_variables.o diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F index a86c24386c8f..24984705d1a4 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics.F @@ -3733,8 +3733,6 @@ subroutine ocn_compute_land_ice_flux_input_fields(layerThickness, normalVelocity !----------------------------------------------------------------- - if ( trim(config_land_ice_flux_mode) == "off" ) return - if ( trim(config_land_ice_draft_mode) == "data" ) then landIceDraftForSsh = landIceDraft @@ -3758,7 +3756,7 @@ subroutine ocn_compute_land_ice_flux_input_fields(layerThickness, normalVelocity #endif endif - if ( trim(config_land_ice_flux_mode) == "pressure_only") return + if ( trim(config_land_ice_flux_mode) == "off") return call mpas_timer_start("land_ice_diagnostic_fields", .false.) 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 abdc7c1ff9b4..b6bac9e911a6 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F @@ -75,9 +75,9 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:), pointer :: pressureAdjustedSSH real (kind=RKIND), dimension(:), pointer :: pgf_sal real (kind=RKIND), dimension(:), pointer :: temp_twd - real (kind=RKIND), dimension(:), pointer :: normalCoeffTWD - real (kind=RKIND), dimension(:), pointer :: tangentialCoeffTWD - real (kind=RKIND), dimension(:), pointer :: topographicWaveDrag + real (kind=RKIND), dimension(:), pointer :: normalCoeffTWD + real (kind=RKIND), dimension(:), pointer :: tangentialCoeffTWD + real (kind=RKIND), dimension(:), pointer :: topographicWaveDrag real (kind=RKIND), dimension(:,:), pointer :: tracersSurfaceValue real (kind=RKIND), dimension(:,:), pointer :: tracersSurfaceLayerValue @@ -101,10 +101,12 @@ module ocn_diagnostics_variables real (kind=RKIND), dimension(:), pointer :: bed_slope_edges real (kind=RKIND), dimension(:), pointer :: lonGradEdge real (kind=RKIND), dimension(:), pointer :: latGradEdge + real (kind=RKIND), dimension(:,:), pointer :: landIceInterfaceTracers real (kind=RKIND), dimension(:,:), pointer :: landIceBoundaryLayerTracers real (kind=RKIND), dimension(:,:), pointer :: landIceTracerTransferVelocities - integer :: indexBLT, indexBLS, & + integer :: indexIT, indexIS, & + indexBLT, indexBLS, & indexHeatTrans, indexSaltTrans real (kind=RKIND), dimension(:), pointer :: penetrativeTemperatureFluxOBL @@ -260,6 +262,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e ! Local pointers for pool retrievals integer, pointer :: & + indexITPtr, indexISPtr, & indexBLTPtr, indexBLSPtr, & indexHeatTransPtr, indexSaltTransPtr, & index_vertNonLocalFluxTempPtr, & @@ -363,6 +366,10 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e call mpas_pool_get_array(diagnosticsPool, 'topDrag', topDrag) call mpas_pool_get_array(diagnosticsPool, 'landIceDraftForSsh', landIceDraftForSsh) call mpas_pool_get_array(diagnosticsPool, 'topDragMagnitude', topDragMag) + call mpas_pool_get_dimension(diagnosticsPool, & + 'index_landIceInterfaceTemperature', indexITPtr) + call mpas_pool_get_dimension(diagnosticsPool, & + 'index_landIceInterfaceSalinity', indexISPtr) call mpas_pool_get_dimension(diagnosticsPool, & 'index_landIceBoundaryLayerTemperature', indexBLTPtr) call mpas_pool_get_dimension(diagnosticsPool, & @@ -371,11 +378,15 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e 'index_landIceHeatTransferVelocity', indexHeatTransPtr) call mpas_pool_get_dimension(diagnosticsPool, & 'index_landIceSaltTransferVelocity', indexSaltTransPtr) + indexIT = indexITPtr + indexIS = indexISPtr indexBLT = indexBLTPtr indexBLS = indexBLSPtr indexHeatTrans = indexHeatTransPtr indexSaltTrans = indexSaltTransPtr + call mpas_pool_get_array(diagnosticsPool, & + 'landIceInterfaceTracers', landIceInterfaceTracers) call mpas_pool_get_array(diagnosticsPool, & 'landIceBoundaryLayerTracers', landIceBoundaryLayerTracers) call mpas_pool_get_array(diagnosticsPool, & @@ -755,6 +766,7 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e !$acc enter data create( & !$acc topDrag, & !$acc topDragMag, & + !$acc landIceInterfaceTracers, & !$acc landIceBoundaryLayerTracers, & !$acc landIceTracerTransferVelocities, & !$acc landIceFrictionVelocity, & @@ -1017,6 +1029,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ !$acc exit data delete( & !$acc topDrag, & !$acc topDragMag, & + !$acc landIceInterfaceTracers, & !$acc landIceBoundaryLayerTracers, & !$acc landIceTracerTransferVelocities, & !$acc landIceDraftForSsh, & @@ -1247,6 +1260,7 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{ if ( trim(config_land_ice_flux_mode) /= 'off' ) then nullify(topDrag, & topDragMag, & + landIceInterfaceTracers, & landIceBoundaryLayerTracers, & landIceTracerTransferVelocities, & landIceFrictionVelocity, & diff --git a/components/mpas-ocean/src/shared/mpas_ocn_effective_density_in_land_ice.F b/components/mpas-ocean/src/shared/mpas_ocn_effective_density_in_land_ice.F deleted file mode 100644 index a8b5ca9f053c..000000000000 --- a/components/mpas-ocean/src/shared/mpas_ocn_effective_density_in_land_ice.F +++ /dev/null @@ -1,216 +0,0 @@ -! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) -! and the University Corporation for Atmospheric Research (UCAR). -! -! Unless noted otherwise source code is licensed under the BSD license. -! Additional copyright and license information can be found in the LICENSE file -! distributed with this code, or at http://mpas-dev.github.io/license.html -! -!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| -! -! ocn_effective_density_in_land_ice -! -!> \brief MPAS ocean effective density in land ice -!> \author Xylar Asay-Davis -!> \date 10/03/2015 -!> \details -!> This module contains routines for computing the effective seawater -!> density in land ice using Arhimedes' principle. -! -!----------------------------------------------------------------------- - -module ocn_effective_density_in_land_ice - - use mpas_constants - use mpas_kind_types - use mpas_derived_types - use mpas_pool_routines - - use ocn_constants - use ocn_config - use ocn_mesh - - implicit none - private - save - - !-------------------------------------------------------------------- - ! - ! Public parameters - ! - !-------------------------------------------------------------------- - - !-------------------------------------------------------------------- - ! - ! Public member functions - ! - !-------------------------------------------------------------------- - - public :: ocn_effective_density_in_land_ice_update - - !-------------------------------------------------------------------- - ! - ! Private module variables - ! - !-------------------------------------------------------------------- - -!*********************************************************************** - -contains - -!*********************************************************************** -! -! routine ocn_effective_density_in_land_ice_update -! -!> \brief updates effective density in land ice -!> \author Xylar Asay-Davis -!> \date 10/03/2015 -!> \details -!> This routine updates the value of the effective seawater density -!> displaced by land ice, based on Archimedes' principle. The -!> effective density is smoothed and extrapolated by averaging with -!> nearest neighbors (cellsOnCell). -! -!----------------------------------------------------------------------- - - subroutine ocn_effective_density_in_land_ice_update(forcingPool, & - statePool, ierr)!{{{ - - !----------------------------------------------------------------- - ! input variables - !----------------------------------------------------------------- - - type (mpas_pool_type), intent(in) :: & - forcingPool !< [in] Forcing information - - !----------------------------------------------------------------- - ! input/output variables - !----------------------------------------------------------------- - - type (mpas_pool_type), intent(inout) :: & - statePool !< [inout] state information - - !----------------------------------------------------------------- - ! output variables - !----------------------------------------------------------------- - - integer, intent(out) :: ierr !< [out] Error flag - - !----------------------------------------------------------------- - ! local variables - !----------------------------------------------------------------- - - integer :: & - i, iCell, &! loop indices for neighbor loop and cell loop - cell2 ! neighbor cell address - - real (kind=RKIND) :: & - weightSum ! local sum of weights for masked points - - real (kind=RKIND), dimension(:), pointer :: & - ssh, &! sea surface height - landIcePressure, &! land ice pressure - effectiveDensityCur, &! effective density at current time - effectiveDensityNew ! effective density at new time level - - integer, dimension(:), pointer :: & - landIceFloatingMask ! mask for land ice presence on cell - - ! Scratch Arrays - ! effectiveDensityScratch: effective seawater density in land - ! ice before horizontal averaging units: kg m^{-3} - real (kind=RKIND), dimension(:), allocatable :: & - effectiveDensityScratch - - ! End preamble - !----------------------------------------------------------------- - ! Begin code - - ! Set default return code and exit if not needed - ierr = 0 - if ( (trim(config_land_ice_flux_mode) /= 'coupled') ) return - - ! Retrieve forcing and state variables - call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', & - landIceFloatingMask) - call mpas_pool_get_array(forcingPool, 'landIcePressure', & - landIcePressure) - call mpas_pool_get_array(statePool, 'ssh', ssh, 2) - call mpas_pool_get_array(statePool, 'effectiveDensityInLandIce', & - effectiveDensityCur, 1) - call mpas_pool_get_array(statePool, 'effectiveDensityInLandIce', & - effectiveDensityNew, 2) - !$acc enter data copyin(landIceFloatingMask, landIcePressure, ssh, & - !$acc effectiveDensityCur, effectiveDensityNew) - - allocate(effectiveDensityScratch(nCellsAll)) - !$acc enter data create(effectiveDensityScratch) - - ! Compute effective density in each cell - ! TODO: should only apply to floating land ice, - ! once wetting/drying is supported -#ifdef MPAS_OPENACC - !$acc parallel loop & - !$acc present(effectiveDensityScratch, effectiveDensityCur, & - !$acc landIcePressure, landIceFloatingMask, ssh) -#else - !$omp parallel - !$omp do schedule(runtime) -#endif - do iCell = 1, nCellsAll - if (landIceFloatingMask(iCell) == 1) then - ! land ice is present to update the effective density - effectiveDensityScratch(iCell) = -landIcePressure(iCell)/ & - (ssh(iCell)*gravity) - else - ! we copy the previous effective density - effectiveDensityScratch(iCell) = effectiveDensityCur(iCell) - end if - end do -#ifndef MPAS_OPENACC - !$omp end do -#endif - - ! smooth/extrapolate density by averaging with nearest neighbors -#ifdef MPAS_OPENACC - !$acc parallel loop & - !$acc present(effectiveDensityNew, effectiveDensityScratch, & - !$acc nEdgesOnCell, cellsOnCell, cellMask) & - !$acc private(weightSum) -#else - !$omp do schedule(runtime) private(weightSum, i, cell2) -#endif - do iCell = 1, nCellsOwned - weightSum = 1.0_RKIND - effectiveDensityNew(iCell) = effectiveDensityScratch(iCell) - do i = 1, nEdgesOnCell(iCell) - cell2 = cellsOnCell(i,iCell) - if (cell2 <= nCellsAll) then - effectiveDensityNew(iCell) = effectiveDensityNew(iCell) & - + cellMask(1,cell2)* & - effectiveDensityScratch(cell2) - weightSum = weightSum + cellMask(1,cell2) - end if - end do - effectiveDensityNew(iCell) = effectiveDensityNew(iCell)/ & - weightSum - end do -#ifndef MPAS_OPENACC - !$omp end do - !$omp end parallel -#endif - - !$acc exit data delete(effectiveDensityScratch) - !$acc exit data delete(landIceFloatingMask, landIcePressure, ssh, & - !$acc effectiveDensityCur, effectiveDensityNew) - deallocate(effectiveDensityScratch) - - !-------------------------------------------------------------------- - - end subroutine ocn_effective_density_in_land_ice_update !}}} - -!*********************************************************************** - -end module ocn_effective_density_in_land_ice - -!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| -! vim: foldmethod=marker diff --git a/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F b/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F index db5df3dfbe00..b3b72037e1ca 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_surface_land_ice_fluxes.F @@ -69,7 +69,6 @@ module ocn_surface_land_ice_fluxes logical :: & landIceDataOn , &! flag to determine if using data (prescribed) fluxes - landIceStandaloneOn, &! flag to determine if using standalone fluxes landIceCoupledOn, &! flag to determine if using fluxes from the coupler useHollandJenkinsAdvDiff ! use Holland-Jenkins advection-diffusion @@ -392,6 +391,7 @@ subroutine ocn_surface_land_ice_fluxes_active_tracers(meshPool, forcingPool, tra integer, pointer :: nCells real (kind=RKIND), dimension(:), pointer :: landIceHeatFlux, & + landIceCouplingHeatFlux, & dataLandIceHeatFlux real (kind=RKIND), pointer :: runningMeanRemovedIceRunoff, & @@ -407,6 +407,7 @@ subroutine ocn_surface_land_ice_fluxes_active_tracers(meshPool, forcingPool, tra call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux) + call mpas_pool_get_array(forcingPool, 'landIceCouplingHeatFlux', landIceCouplingHeatFlux) if ( landIceDataOn ) then @@ -424,6 +425,7 @@ subroutine ocn_surface_land_ice_fluxes_active_tracers(meshPool, forcingPool, tra !$omp do schedule(runtime) do iCell = 1, nCells landIceHeatFlux(iCell) = scaling * dataLandIceHeatFlux(iCell) + landIceCouplingHeatFlux(iCell) = landIceHeatFlux(iCell) end do !$omp end do !$omp end parallel @@ -498,18 +500,14 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & real (kind=RKIND), dimension(:), pointer :: landIcePressure, landIceFloatingFraction, & landIceSurfaceTemperature, & landIceFreshwaterFlux, & - landIceHeatFlux, heatFluxToLandIce + landIceConductiveHeatFlux, & + landIceHeatFlux, landIceCouplingHeatFlux, & + heatFluxToLandIce real (kind=RKIND), dimension(:), pointer :: frazilIceFreshwaterFlux, landIceFreshwaterFluxTotal integer, dimension(:), pointer :: landIceFloatingMask - real (kind=RKIND), dimension(:,:), pointer :: & - landIceInterfaceTracers - - integer, pointer :: indexITPtr, indexISPtr - integer :: indexIT, indexIS - ! Scratch Arrays !--- ! freezeInterfaceSalinity: salinity at land ice-ocean interface where freezing is occurring @@ -524,7 +522,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & err = 0 - if (.not. (landIceStandaloneOn .or. landIceDataOn)) return + if (.not. (landIceFluxesOn .or. landIceDataOn)) return call mpas_timer_start("land_ice_build_arrays") @@ -534,26 +532,17 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & call mpas_pool_get_array(forcingPool, 'frazilIceFreshwaterFlux', frazilIceFreshwaterFlux) call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFluxTotal', landIceFreshwaterFluxTotal) call mpas_pool_get_array(forcingPool, 'landIceFloatingMask', landIceFloatingMask) + call mpas_pool_get_array(forcingPool, 'landIceCouplingHeatFlux', landIceCouplingHeatFlux) + call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce) nCells = nCellsArray( size(nCellsArray) ) + landIceCouplingHeatFlux(:) = 0.0_RKIND - if (landIceStandaloneOn) then + if (landIceFluxesOn .and. .not. landIceDataOn) then call mpas_pool_get_array(forcingPool, 'landIcePressure', landIcePressure) call mpas_pool_get_array(forcingPool, 'landIceFloatingFraction', landIceFloatingFraction) call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux) - call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce) - - - call mpas_pool_get_array(forcingPool, 'landIceInterfaceTracers', landIceInterfaceTracers) - call mpas_pool_get_dimension(forcingPool, & - 'index_landIceInterfaceTemperature', & - indexITPtr) - call mpas_pool_get_dimension(forcingPool, & - 'index_landIceInterfaceSalinity', & - indexISPtr) - indexIT = indexITPtr - indexIS = indexISPtr if (useHollandJenkinsAdvDiff) then call mpas_pool_get_array(forcingPool, 'landIceSurfaceTemperature', landIceSurfaceTemperature) @@ -594,6 +583,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & * (landIceInterfaceTracers(indexIT,iCell)-landIceBoundaryLayerTracers(indexBLT,iCell))) landIceHeatFlux(iCell) = landIceFloatingFraction(iCell)*heatFlux + landIceCouplingHeatFlux(iCell) = 0.0_RKIND heatFluxToLandIce(iCell) = 0.0_RKIND end do @@ -616,6 +606,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & landIceInterfaceTracers(indexIS,:), & landIceFreshwaterFlux, & landIceHeatFlux, & + landIceCouplingHeatFlux, & heatFluxToLandIce, & nCells, & err) @@ -653,11 +644,14 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & landIceInterfaceTracers(indexIT,iCell) = freezeInterfaceTemperature(iCell) landIceFreshwaterFlux(iCell) = freezeFreshwaterFlux(iCell) landIceHeatFlux(iCell) = freezeHeatFlux(iCell) + landIceCouplingHeatFlux(iCell) = freezeIceHeatFlux(iCell) heatFluxToLandIce(iCell) = freezeIceHeatFlux(iCell) end do else ! not using Holland and Jenkins advection/diffusion - call compute_melt_fluxes( & + call mpas_pool_get_array(forcingPool, 'landIceConductiveHeatFlux', landIceConductiveHeatFlux) + call compute_melt_fluxes_with_conduction( & landIceFloatingMask, & + landIceConductiveHeatFlux, & landIceBoundaryLayerTracers(indexBLT,:), & landIceBoundaryLayerTracers(indexBLS,:), & landIceTracerTransferVelocities(indexHeatTrans,:), & @@ -667,6 +661,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & landIceInterfaceTracers(indexIS,:), & landIceFreshwaterFlux, & landIceHeatFlux, & + landIceCouplingHeatFlux, & heatFluxToLandIce, & nCells, & err) @@ -683,6 +678,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & landIceFreshwaterFlux(iCell) = landIceFloatingFraction(iCell)*landIceFreshwaterFlux(iCell) landIceHeatFlux(iCell) = landIceFloatingFraction(iCell)*landIceHeatFlux(iCell) + landIceCouplingHeatFlux(iCell) = landIceFloatingFraction(iCell)*landIceCouplingHeatFlux(iCell) heatFluxToLandIce(iCell) = landIceFloatingFraction(iCell)*heatFluxToLandIce(iCell) end do @@ -696,7 +692,7 @@ subroutine ocn_surface_land_ice_fluxes_build_arrays(meshPool, & freezeIceHeatFlux) end if - end if ! landIceStandaloneOn + end if ! landIceFluxesOn ! Add frazil and interface melt/freeze to get total fresh water flux if ( associated(frazilIceFreshwaterFlux) ) then @@ -779,7 +775,7 @@ subroutine ocn_surface_land_ice_fluxes_accumulate_fluxes(meshPool, & call mpas_pool_get_dimension(meshPool, 'nCellsArray', nCellsArray) nCells = nCellsArray( size(nCellsArray) ) - if (landIceStandaloneOn .or. landIceDataOn) then + if (landIceFluxesOn .or. landIceDataOn) then call mpas_pool_get_array(statePool, 'accumulatedLandIceMass', accumulatedLandIceMassNew, 2) call mpas_pool_get_array(statePool, 'accumulatedLandIceMass', accumulatedLandIceMassOld, 1) call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFlux) @@ -789,13 +785,13 @@ subroutine ocn_surface_land_ice_fluxes_accumulate_fluxes(meshPool, & end do end if - if (landIceStandaloneOn) then + if (landIceFluxesOn .and. .not. landIceDataOn) then call mpas_pool_get_array(forcingPool, 'heatFluxToLandIce', heatFluxToLandIce) call mpas_pool_get_array(statePool, 'accumulatedLandIceHeat', accumulatedLandIceHeatNew, 2) call mpas_pool_get_array(statePool, 'accumulatedLandIceHeat', accumulatedLandIceHeatOld, 1) ! accumulate land-ice heat do iCell = 1, nCells - accumulatedLandIceHeatNew(iCell) = accumulatedLandIceHeatOld(iCell) + dt*heatFluxToLandIce(iCell) + accumulatedLandIceHeatNew(iCell) = accumulatedLandIceHeatOld(iCell) - dt*heatFluxToLandIce(iCell) end do end if @@ -835,48 +831,32 @@ subroutine ocn_surface_land_ice_fluxes_init(err)!{{{ isomipOn = .false. jenkinsOn = .false. hollandJenkinsOn = .false. + landIcePressureOn = .false. cpLandIce = 0.0_RKIND rhoLandIce = 0.0_RKIND ISOMIPgammaT = 0.0_RKIND !*** Determine whether land ice fluxes are on + if ( config_use_land_ice_pressure ) & + landIcePressureOn = .true. + select case (trim(config_land_ice_flux_mode)) case ('data','Data','DATA') - landIcePressureOn = .true. landIceFluxesOn = .true. landIceDataOn = .true. - landIceStandaloneOn = .false. - landIceCoupledOn = .false. - case ('standalone','Standalone','STANDALONE') - landIcePressureOn = .true. - landIceFluxesOn = .true. - landIceDataOn = .false. - landIceStandaloneOn = .true. - landIceCoupledOn = .false. - case ('coupled','Coupled','COUPLED') - landIcePressureOn = .true. + case ('active','Active','ACTIVE') landIceFluxesOn = .true. landIceDataOn = .false. - landIceStandaloneOn = .false. - landIceCoupledOn = .true. - case ('pressure_only','Pressure_Only','PRESSURE_ONLY') - landIcePressureOn = .true. - landIceFluxesOn = .false. - landIceDataOn = .false. - landIceStandaloneOn = .false. - landIceCoupledOn = .false. case ('off','Off','OFF') - landIcePressureOn = .false. landIceFluxesOn = .false. landIceDataOn = .false. - landIceStandaloneOn = .false. - landIceCoupledOn = .false. case default call mpas_log_write( ' Error: Incorrect choice of config_land_ice_flux_mode: '// config_land_ice_flux_mode) err = 1 return end select + call mpas_log_write( ' Choice of config_land_ice_flux_mode: '// config_land_ice_flux_mode) !*** Determine which flux formulation to use @@ -1041,7 +1021,12 @@ subroutine compute_melt_fluxes( & ! 2. the diffusion (if any) of heat into the ice, based on temperature difference between ! the reference point in the ice (either the surface or the middle of the bottom layer) ! and the interface - outIceHeatFlux(iCell) = -cpLandIce*outFreshwaterFlux(iCell)*outInterfaceTemperature(iCell) + ! NB: sign convention is the same as outOceanHeatFlux, positive into the ocean + if (outFreshwaterFlux(iCell) > 0) then ! melting + outIceHeatFlux(iCell) = cpLandIce*outInterfaceTemperature(iCell)*outFreshwaterFlux(iCell) + else + outIceHeatFlux(iCell) = cp_sw*outInterfaceTemperature(iCell)*outFreshwaterFlux(iCell) + endif end do !$omp end do !$omp end parallel @@ -1050,6 +1035,139 @@ subroutine compute_melt_fluxes( & end subroutine compute_melt_fluxes !}}} + subroutine compute_melt_fluxes_with_conduction( & + mask, & + landIceConductiveHeatFlux, & + oceanTemperature, & + oceanSalinity, & + oceanHeatTransferVelocity, & + oceanSaltTransferVelocity, & + interfacePressure, & + outInterfaceTemperature, & + outInterfaceSalinity, & + outFreshwaterFlux, & + outOceanHeatFlux, & + outCouplingHeatFlux, & + outIceHeatFlux, & + nCells, & + err) !{{{ + + !----------------------------------------------------------------- + ! + ! input variables + ! + !----------------------------------------------------------------- + + integer, dimension(:), intent(in) :: & + mask !< Input: mask for land-ice fluxes + + real (kind=RKIND), dimension(:), intent(in) :: & + landIceConductiveHeatFlux, & !< Input: conductive heat flux into the base of ice + oceanTemperature, & !< Input: ocean temperature in top layer + oceanSalinity, & !< Input: ocean salinity in top layer + oceanHeatTransferVelocity, & !< Input: ocean heat transfer velocity + oceanSaltTransferVelocity, & !< Input: ocean salt transfer velocity + interfacePressure !< Input: pressure at the ice-ocean interface + + integer, intent(in) :: nCells !< Input: number of cells in each array + + !----------------------------------------------------------------- + ! + ! output variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:), intent(out) :: & + outInterfaceTemperature, & !< Output: ice/ocean temperature at the interface + outInterfaceSalinity, & !< Output: ocean salinity at the interface + outFreshwaterFlux, & !< Output: ocean thickness flux (melt rate) + outOceanHeatFlux, & !< Output: the temperature flux into the ocean + outCouplingHeatFlux, & !< Output: the glacier-facing coupled heat flux + outIceHeatFlux !< Output: the temperature flux into the ice + + integer, intent(out) :: err !< Output: Error flag + + !----------------------------------------------------------------- + ! + ! local variables + ! + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:), allocatable :: outAdvectiveHeatFlux + + real (kind=RKIND) :: T0, transferVelocityRatio, Tlatent, a, b, c, & + dTf_dS + integer :: iCell + + real (kind=RKIND), parameter :: minInterfaceSalinity = 0.001_RKIND + + err = 0 + Tlatent = latent_heat_fusion_mks/cp_sw + + allocate(outAdvectiveHeatFlux(nCells)) + + !$omp parallel + !$omp do schedule(runtime) private(T0, dTf_dS, transferVelocityRatio, a, b, c) + do iCell = 1, nCells + if (mask(iCell) == 0) cycle + + T0 = ocn_freezing_temperature(salinity=0.0_RKIND, pressure=interfacePressure(iCell), & + inLandIceCavity=.true.) + dTf_dS = ocn_freezing_temperature_salinity_deriv(salinity=0.0_RKIND, pressure=interfacePressure(iCell), & + inLandIceCavity=.true.) + + transferVelocityRatio = oceanSaltTransferVelocity(iCell)/oceanHeatTransferVelocity(iCell) + + a = -dTf_dS + b = transferVelocityRatio*Tlatent + oceanTemperature(iCell) - T0 + & + landIceConductiveHeatFlux(iCell) / (rho_sw * cp_sw * oceanHeatTransferVelocity(iCell)) + c = -transferVelocityRatio*Tlatent*max(oceanSalinity(iCell), 0.0_RKIND) + + ! a is non-negative; c is strictly non-positive so we never get imaginary roots. + ! Since a can be zero, we need a solution of the quadratic equation for 1/Si instead of Si. + ! Following: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf + ! Since a and -c are are non-negative, the term in the square root is also always >= |b|. + ! In all reasonable cases, b will be strictly positive, since transferVelocityRatio*Tlatent ~ 2 C, + ! T0 ~ -1.8 C and oceanTemperature should never be able to get below about -3 C + ! As long as either b or both a and c are greater than zero, the strictly non-negative root is + outInterfaceSalinity(iCell) = max(-(2.0_RKIND*c)/(b + sqrt(b**2 - 4.0_RKIND*a*c)), minInterfaceSalinity) + + outInterfaceTemperature(iCell) = dTf_dS*outInterfaceSalinity(iCell)+T0 + + outFreshwaterFlux(iCell) = rho_sw*oceanSaltTransferVelocity(iCell) & + * (oceanSalinity(iCell)/outInterfaceSalinity(iCell) - 1.0_RKIND) + + if (outFreshwaterFlux(iCell) > 0) then ! melting + outAdvectiveHeatFlux(iCell) = cpLandIce*outInterfaceTemperature(iCell)*outFreshwaterFlux(iCell) + else + outAdvectiveHeatFlux(iCell) = cp_sw*outInterfaceTemperature(iCell)*outFreshwaterFlux(iCell) + endif + ! The heat fluxes that determine the temperature tendency: + ! 1. the advection of meltwater into the top layer (or removal for freezing) + ! 2. the turbulent transfer of heat across the boundary layer, based on the thermal driving + ! Applied as: tracersSurfaceFlux(1, iCell) += landIceHeatFlux(iCell)/(rho_sw*cp_sw) + outOceanHeatFlux(iCell) = outAdvectiveHeatFlux(iCell) + & + -rho_sw*cp_sw*oceanHeatTransferVelocity(iCell)*(oceanTemperature(iCell)-outInterfaceTemperature(iCell)) + + ! the temperature fluxes across the ice-ocean interface are: + ! 1. the advection of ice at the interface temperature out of the domain due to melting + ! (or in due to freezing) + ! 2. the diffusion (if any) of heat into the ice, based on temperature difference between + ! the reference point in the ice (either the surface or the middle of the bottom layer) + ! and the interface + ! NB: sign convention is the same as outOceanHeatFlux, positive into the ocean + ! Export only the advective part to land ice through the coupler; keep the + ! full ice-side term for MPAS-Ocean internal energy accounting. + outCouplingHeatFlux(iCell) = outAdvectiveHeatFlux(iCell) + outIceHeatFlux(iCell) = outCouplingHeatFlux(iCell) + landIceConductiveHeatFlux(iCell) + end do + !$omp end do + !$omp end parallel + + !-------------------------------------------------------------------- + + end subroutine compute_melt_fluxes_with_conduction !}}} + !*********************************************************************** ! @@ -1187,7 +1305,8 @@ subroutine compute_HJ99_melt_fluxes( & ! Since we're considering only melting and ignoring diffusion, ! the ice loses heat simply by the loss of ice mass at the data ! (surface?) ice temperature - outIceHeatFlux(iCell) = -cpLandIce*outFreshwaterFlux(iCell)*iceTemperature(iCell) + ! NB: sign convention is the same as outOceanHeatFlux, positive into the ocean + outIceHeatFlux(iCell) = cpLandIce*outFreshwaterFlux(iCell)*iceTemperature(iCell) end do !$omp end do !$omp end parallel diff --git a/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F b/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F index 91ceb0de8a4b..db2b4f29e5d2 100644 --- a/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F +++ b/components/mpas-ocean/src/shared/mpas_ocn_time_average_coupled.F @@ -53,12 +53,13 @@ subroutine ocn_time_average_coupled_init(forcingPool)!{{{ type (mpas_pool_type), intent(inout) :: forcingPool real (kind=RKIND), dimension(:,:), pointer :: avgTracersSurfaceValue, avgSurfaceVelocity, avgSSHGradient, & - avgLandIceBoundaryLayerTracers, avgLandIceTracerTransferVelocities + avgLandIceInterfaceTracers - real (kind=RKIND), dimension(:), pointer :: avgEffectiveDensityInLandIce, avgTotalFreshWaterTemperatureFlux, & + real (kind=RKIND), dimension(:), pointer :: avgTotalFreshWaterTemperatureFlux, & avgLandIceFreshwaterFlux, & avgRemovedRiverRunoffFlux, avgRemovedIceRunoffFlux, & - avgLandIceHeatFlux, avgRemovedIceRunoffHeatFlux + avgLandIceHeatFlux, avgLandIceCouplingHeatFlux, & + avgRemovedIceRunoffHeatFlux real (kind=RKIND), dimension(:,:), pointer :: avgThermalForcingAtZLevels real (kind=RKIND), dimension(:,:), pointer :: avgThermalForcingAtZLevelsMask @@ -108,32 +109,30 @@ subroutine ocn_time_average_coupled_init(forcingPool)!{{{ !$omp end do !$omp end parallel - if(trim(config_land_ice_flux_mode) == 'coupled') then - call mpas_pool_get_array(forcingPool, 'avgLandIceBoundaryLayerTracers', avgLandIceBoundaryLayerTracers) - call mpas_pool_get_array(forcingPool, 'avgLandIceTracerTransferVelocities', avgLandIceTracerTransferVelocities) - call mpas_pool_get_array(forcingPool, 'avgEffectiveDensityInLandIce', avgEffectiveDensityInLandIce) + if(trim(config_land_ice_flux_mode) == 'active') then + call mpas_pool_get_array(forcingPool, 'avgLandIceInterfaceTracers', avgLandIceInterfaceTracers) !$omp parallel !$omp do schedule(runtime) do iCell = 1, nCells - avgLandIceBoundaryLayerTracers(:, iCell) = 0.0_RKIND - avgLandIceTracerTransferVelocities(:, iCell) = 0.0_RKIND - avgEffectiveDensityInLandIce(iCell) = 0.0_RKIND + avgLandIceInterfaceTracers(:, iCell) = 0.0_RKIND end do !$omp end do !$omp end parallel end if ! Set up polar fields if necessary - if(trim(config_land_ice_flux_mode)=='standalone' .or. trim(config_land_ice_flux_mode) == 'data') then + if(trim(config_land_ice_flux_mode) == 'active' .or. trim(config_land_ice_flux_mode) == 'data') then call mpas_pool_get_array(forcingPool, 'avgLandIceFreshwaterFlux', avgLandIceFreshwaterFlux) call mpas_pool_get_array(forcingPool, 'avgLandIceHeatFlux', avgLandIceHeatFlux) + call mpas_pool_get_array(forcingPool, 'avgLandIceCouplingHeatFlux', avgLandIceCouplingHeatFlux) !$omp parallel !$omp do schedule(runtime) do iCell = 1, nCells avgLandIceFreshwaterFlux(iCell) = 0.0_RKIND avgLandIceHeatFlux(iCell) = 0.0_RKIND + avgLandIceCouplingHeatFlux(iCell) = 0.0_RKIND end do !$omp end do !$omp end parallel @@ -278,12 +277,11 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel integer, pointer :: index_temperaturePtr, index_SSHzonalPtr, & index_SSHmeridionalPtr, nAccumulatedCoupled, nCells integer :: index_temperature, index_SSHzonal, index_SSHmeridional - real (kind=RKIND), dimension(:,:), pointer :: & - avgLandIceBoundaryLayerTracers, avgLandIceTracerTransferVelocities - real (kind=RKIND), dimension(:), pointer :: effectiveDensityInLandIce, avgEffectiveDensityInLandIce, & - totalFreshWaterTemperatureFlux, avgTotalFreshWaterTemperatureFlux, & + real (kind=RKIND), dimension(:,:), pointer :: avgLandIceInterfaceTracers + real (kind=RKIND), dimension(:), pointer :: totalFreshWaterTemperatureFlux, avgTotalFreshWaterTemperatureFlux, & landIceFreshwaterFlux, avgLandIceFreshwaterFlux, & landIceHeatFlux, avgLandIceHeatFlux, & + landIceCouplingHeatFlux, avgLandIceCouplingHeatFlux, & removedRiverRunoffFlux, avgRemovedRiverRunoffFlux, & removedIceRunoffFlux, avgRemovedIceRunoffFlux, & avgRemovedIceRunoffHeatFlux @@ -361,34 +359,31 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel !$omp end do !$omp end parallel - if(trim(config_land_ice_flux_mode) == 'coupled') then - call mpas_pool_get_array(statePool, 'effectiveDensityInLandIce', effectiveDensityInLandIce, timeLevel) - - call mpas_pool_get_array(forcingPool, 'avgLandIceBoundaryLayerTracers', avgLandIceBoundaryLayerTracers) - call mpas_pool_get_array(forcingPool, 'avgLandIceTracerTransferVelocities', avgLandIceTracerTransferVelocities) - call mpas_pool_get_array(forcingPool, 'avgEffectiveDensityInLandIce', avgEffectiveDensityInLandIce) - + if(trim(config_land_ice_flux_mode) == 'active') then + call mpas_pool_get_array(forcingPool, 'avgLandIceInterfaceTracers', avgLandIceInterfaceTracers) !$omp parallel !$omp do schedule(runtime) do iCell = 1, nCells - avgLandIceBoundaryLayerTracers(:, iCell) = ( avgLandIceBoundaryLayerTracers(:, iCell) * nAccumulatedCoupled & - + landIceBoundaryLayerTracers(:, iCell) ) / ( nAccumulatedCoupled + 1 ) - avgLandIceTracerTransferVelocities(:, iCell) = ( avgLandIceTracerTransferVelocities(:, iCell) * nAccumulatedCoupled & - + landIceTracerTransferVelocities(:, iCell) ) / ( nAccumulatedCoupled + 1) - avgEffectiveDensityInLandIce(iCell) = ( avgEffectiveDensityInLandIce(iCell) * nAccumulatedCoupled & - + effectiveDensityInLandIce(iCell) ) / ( nAccumulatedCoupled + 1) + avgLandIceInterfaceTracers(:, iCell) = avgLandIceInterfaceTracers(:, iCell) * nAccumulatedCoupled & + + landIceInterfaceTracers(:, iCell) + ! avgLandIceInterfaceTemperature needs to be in Kelvin, like avgTemperatureSurfaceValue + avgLandIceInterfaceTracers(index_temperature, iCell) = avgLandIceInterfaceTracers(index_temperature, iCell) & + + T0_Kelvin + avgLandIceInterfaceTracers(:, iCell) = avgLandIceInterfaceTracers(:, iCell) / ( nAccumulatedCoupled + 1 ) end do !$omp end do !$omp end parallel end if ! Accumulate polar fields if necessary - if(trim(config_land_ice_flux_mode) == 'standalone' .or. trim(config_land_ice_flux_mode) == 'data') then + if(trim(config_land_ice_flux_mode) == 'active' .or. trim(config_land_ice_flux_mode) == 'data') then call mpas_pool_get_array(forcingPool, 'avgLandIceFreshwaterFlux', avgLandIceFreshwaterFlux) call mpas_pool_get_array(forcingPool, 'landIceFreshwaterFlux', landIceFreshwaterFlux) call mpas_pool_get_array(forcingPool, 'avgLandIceHeatFlux', avgLandIceHeatFlux) call mpas_pool_get_array(forcingPool, 'landIceHeatFlux', landIceHeatFlux) + call mpas_pool_get_array(forcingPool, 'avgLandIceCouplingHeatFlux', avgLandIceCouplingHeatFlux) + call mpas_pool_get_array(forcingPool, 'landIceCouplingHeatFlux', landIceCouplingHeatFlux) !$omp parallel !$omp do schedule(runtime) @@ -397,6 +392,8 @@ subroutine ocn_time_average_coupled_accumulate(statePool, forcingPool, timeLevel + landIceFreshwaterFlux(iCell) ) / ( nAccumulatedCoupled + 1) avgLandIceHeatFlux(iCell) = ( avgLandIceHeatFlux(iCell) * nAccumulatedCoupled & + landIceHeatFlux(iCell) ) / ( nAccumulatedCoupled + 1) + avgLandIceCouplingHeatFlux(iCell) = ( avgLandIceCouplingHeatFlux(iCell) * nAccumulatedCoupled & + + landIceCouplingHeatFlux(iCell) ) / ( nAccumulatedCoupled + 1) end do !$omp end do !$omp end parallel diff --git a/driver-mct/main/cime_comp_mod.F90 b/driver-mct/main/cime_comp_mod.F90 index b0b4822141b3..c09c2efff12c 100644 --- a/driver-mct/main/cime_comp_mod.F90 +++ b/driver-mct/main/cime_comp_mod.F90 @@ -146,7 +146,7 @@ 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_glc_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 @@ -754,9 +754,9 @@ subroutine cime_pre_init1(esmf_log_option) integer(i8) :: irtc_rate ! factor to convert time to seconds integer :: tmode ! Thread mode provided by the MPI library - + beg_count = shr_sys_irtc(irtc_rate) - + #if defined(MPINIT_WORKAROUND) && (MPINIT_WORKAROUND == 1) call atm_init_hip_mct() #endif @@ -770,7 +770,7 @@ subroutine cime_pre_init1(esmf_log_option) end_count = shr_sys_irtc(irtc_rate) mpi_init_time = real( (end_count-beg_count), r8)/real(irtc_rate, r8) - + call mpi_comm_dup(MPI_COMM_WORLD, global_comm, ierr) call shr_mpi_chkerr(ierr,subname//' mpi_comm_dup') @@ -1452,9 +1452,9 @@ subroutine cime_pre_init2() !Print BFBFLAG value in the log file if (iamroot_CPLID) then call seq_infodata_GetData(infodata, bfbflag=bfbflag) - write(logunit,'(2A,L4)') subname,'BFBFLAG is:',bfbflag + write(logunit,'(2A,L4)') subname,'BFBFLAG is:',bfbflag endif - + call t_stopf('CPL:cime_pre_init2') ! CPL:cime_pre_init2 timer elapsed time will be double counted @@ -2361,7 +2361,7 @@ subroutine cime_init() call seq_nlmap_init_a2oi_cons(mapper_lcl, fractions_ax) end if end if - + !---------------------------------------------------------- !| ATM PREP for recalculation of initial solar @@ -2490,7 +2490,7 @@ subroutine cime_init() write(logunit,103) subname,' Reading restart file ',trim(rest_file) call shr_sys_flush(logunit) end if - + call t_startf('CPL:seq_rest_read-init') call seq_rest_read(rest_file, infodata, & atm, lnd, ice, ocn, rof, glc, wav, esp, iac, & @@ -4237,7 +4237,7 @@ subroutine cime_run_iac_recv_post() info_debug=info_debug, timer_diag='CPL:iacpost_diagav') ! Need to reset our max vector to zero for the following year - if (lnd_present) then + if (lnd_present) then call prep_iac_zero_max() endif @@ -4305,28 +4305,11 @@ subroutine cime_run_ocnglc_coupling() if (glc_present) then - ! create o2x_gx for either ocn-glc coupling or ocn-glc shelf coupling - if (ocn_c2_glctf .or. (ocn_c2_glcshelf .and. glcshelf_c2_ocn)) then + ! create o2x_gx for ocn-glc shelf coupling + if (ocn_c2_glctf) then call prep_glc_calc_o2x_gx(ocn_c2_glctf, ocn_c2_glcshelf, timer='CPL:glcprep_ocn2glc') !remap ocean fields to o2x_g at ocean couping interval endif - ! if ice-shelf coupling is on, now proceed to handle those calculations here in the coupler - if (ocn_c2_glcshelf .and. glcshelf_c2_ocn) then - ! the boundary flux calculations done in the coupler require inputs from both GLC and OCN, - ! so they will only be valid if both OCN->GLC and GLC->OCN - - call prep_glc_calculate_subshelf_boundary_fluxes ! this is actual boundary layer flux calculation - !this outputs - !x2g_g/g2x_g, where latter is going - !to ocean, so should get remapped to - !ocean grid in prep_ocn_shelf_calc_g2x_ox - call prep_ocn_shelf_calc_g2x_ox(timer='CPL:glcpost_glcshelf2ocn') - !Map g2x_gx shelf fields that were updated above, to g2x_ox. - !Do this at intrinsic coupling - !frequency - call prep_glc_accum_ocn(timer='CPL:glcprep_accum_ocn') !accum x2g_g fields here into x2g_gacc - endif - if (glcshelf_c2_ice) then call prep_ice_shelf_calc_g2x_ix(timer='CPL:glcpost_glcshelf2ice') !Map g2x_gx shelf fields to g2x_ix. @@ -4362,7 +4345,7 @@ subroutine cime_run_lnd_setup_send() call prep_lnd_calc_z2x_lx(timer='CPL:lndprep_iac2lnd') endif - if (ocn_c2_lnd) then + if (ocn_c2_lnd) then call prep_lnd_accum_avg(timer='CPL:lndprep_o2xavg') call prep_lnd_calc_o2x_lx(timer='CPL:lndprep_ocn2lnd') endif @@ -4835,7 +4818,7 @@ subroutine cime_run_calc_budgets1(in_cplrun) if (present(in_cplrun)) then lcplrun = .not. in_cplrun endif - + if (iamin_CPLID) then call cime_comp_barriers(mpicom=mpicom_CPLID, timer='CPL:BUDGET1_BARRIER') call t_drvstartf ('CPL:BUDGET1',cplrun=lcplrun,budget=.true.,barrier=mpicom_CPLID) @@ -4879,7 +4862,7 @@ subroutine cime_run_calc_budgets2(in_cplrun) if (present(in_cplrun)) then lcplrun = .not. in_cplrun endif - + if (iamin_CPLID) then call cime_comp_barriers(mpicom=mpicom_CPLID, timer='CPL:BUDGET2_BARRIER') @@ -4951,7 +4934,7 @@ subroutine cime_run_calc_budgets3(in_cplrun) if (present(in_cplrun)) then lcplrun = .not. in_cplrun endif - + if (iamin_CPLID) then call cime_comp_barriers(mpicom=mpicom_CPLID, timer='CPL:BUDGET0_BARRIER') call t_drvstartf ('CPL:BUDGET0',cplrun=lcplrun,budget=.true.,barrier=mpicom_CPLID) @@ -5556,7 +5539,7 @@ subroutine rpointer_manage(force_remove) ! the steps required for rpointer file consistency based on state. logical, intent(in) :: force_remove ! force removal of .prev files in this call - + integer :: i, n, rcode, unit logical :: previous_rings, file_exists character(32) :: buf diff --git a/driver-mct/main/prep_glc_mod.F90 b/driver-mct/main/prep_glc_mod.F90 index 614df63c6479..ab672ec796e6 100644 --- a/driver-mct/main/prep_glc_mod.F90 +++ b/driver-mct/main/prep_glc_mod.F90 @@ -58,8 +58,6 @@ module prep_glc_mod public :: prep_glc_get_mapper_So2g_shelf public :: prep_glc_get_mapper_Fo2g_shelf - public :: prep_glc_calculate_subshelf_boundary_fluxes - !-------------------------------------------------------------------------- ! Private interfaces !-------------------------------------------------------------------------- @@ -88,9 +86,9 @@ module prep_glc_mod type(mct_aVect), pointer :: o2x_gx(:) ! Ocn export, glc grid, cpl pes - allocated in driver ! accumulation variables - + type(mct_aVect), pointer :: x2gacc_gx(:) ! Glc export, glc grid, cpl pes - allocated in driver - integer , target :: x2gacc_gx_cnt ! x2gacc_gx: number of time samples accumulated + integer , target :: x2gacc_gx_cnt ! x2gacc_gx: number of time samples accumulated 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 @@ -119,20 +117,6 @@ module prep_glc_mod type(mct_aVect), pointer :: o2gacc_ox(:) ! Ocn export, lnd grid, cpl pes - allocated in driver integer , target :: o2gacc_ox_cnt ! number of time samples accumulated - real(r8), allocatable :: oceanTemperature(:) - real(r8), allocatable :: oceanSalinity(:) - real(r8), allocatable :: oceanHeatTransferVelocity(:) - real(r8), allocatable :: oceanSaltTransferVelocity(:) - real(r8), allocatable :: interfacePressure(:) - real(r8), allocatable :: iceTemperature(:) - real(r8), allocatable :: iceTemperatureDistance(:) - integer, allocatable :: iceFloatingMask(:) - real(r8), allocatable :: outInterfaceSalinity(:) - real(r8), allocatable :: outInterfaceTemperature(:) - real(r8), allocatable :: outFreshwaterFlux(:) - real(r8), allocatable :: outOceanHeatFlux(:) - real(r8), allocatable :: outIceHeatFlux(:) - !================================================================================================ contains @@ -311,22 +295,6 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glctf, ocn_c2_glcshelf) call seq_map_init_rcfile(mapper_Fo2g_shelf, ocn(1), glc(1), & 'seq_maps.rc','ocn2glc_shelf_fmapname:','ocn2glc_shelf_fmaptype:',samegrid_go, & 'mapper_Fo2g_shelf initialization',esmf_map_flag) - !Initialize module-level arrays associated with compute_melt_fluxes - allocate(oceanTemperature(lsize_g)) - allocate(oceanSalinity(lsize_g)) - allocate(oceanHeatTransferVelocity(lsize_g)) - allocate(oceanSaltTransferVelocity(lsize_g)) - allocate(interfacePressure(lsize_g)) - allocate(iceTemperature(lsize_g)) - allocate(iceTemperatureDistance(lsize_g)) - allocate(iceFloatingMask(lsize_g)) - allocate(outInterfaceSalinity(lsize_g)) - allocate(outInterfaceTemperature(lsize_g)) - allocate(outFreshwaterFlux(lsize_g)) - allocate(outOceanHeatFlux(lsize_g)) - allocate(outIceHeatFlux(lsize_g)) - ! TODO: Can we allocate these only while used or are we worried about performance hit? - ! TODO: add deallocates! end if @@ -595,12 +563,16 @@ subroutine prep_glc_merge_ocn_forcing( o2x_g, fractions_g, x2g_g ) !----------------------------------------------------------------------- integer :: num_flux_fields - integer :: num_state_fields + integer :: num_shelf_state_fields + integer :: num_tf_state_fields integer :: nflds integer :: i,n integer :: mrgstr_index integer :: index_o2x integer :: index_x2g + integer :: index_x2g_Fogx_qicehi + integer :: index_x2g_Foxo_ismh + integer :: index_x2g_Foxo_q_li integer :: index_ofrac integer :: lsize logical :: iamroot @@ -614,19 +586,22 @@ subroutine prep_glc_merge_ocn_forcing( o2x_g, fractions_g, x2g_g ) call seq_comm_getdata(CPLID, iamroot=iamroot) lsize = mct_aVect_lsize(x2g_g) - !num_flux_fields = shr_string_listGetNum(trim(seq_flds_x2g_fluxes_from_ocn)) - num_flux_fields = 0 - num_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_tf_states_from_ocn)) + num_flux_fields = shr_string_listGetNum(trim(seq_flds_x2g_shelf_fluxes_from_ocn)) + num_shelf_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_shelf_states_from_ocn)) + num_tf_state_fields = shr_string_listGetNum(trim(seq_flds_x2g_tf_states_from_ocn)) + index_x2g_Fogx_qicehi = mct_aVect_indexRA(x2g_g, 'Fogx_qicehi') + index_x2g_Foxo_ismh = mct_aVect_indexRA(x2g_g, 'Foxo_ismh', perrWith='quiet') + index_x2g_Foxo_q_li = mct_aVect_indexRA(x2g_g, 'Foxo_q_li', perrWith='quiet') if (first_time) then - nflds = num_flux_fields + num_state_fields + nflds = num_flux_fields + num_shelf_state_fields + num_tf_state_fields + 1 allocate(mrgstr(nflds)) end if mrgstr_index = 1 - do i = 1, num_state_fields - call seq_flds_getField(field, i, seq_flds_x2g_tf_states_from_ocn) + do i = 1, num_shelf_state_fields + call seq_flds_getField(field, i, seq_flds_x2g_shelf_states_from_ocn) index_o2x = mct_aVect_indexRA(o2x_g, trim(field)) index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) @@ -642,38 +617,52 @@ subroutine prep_glc_merge_ocn_forcing( o2x_g, fractions_g, x2g_g ) mrgstr_index = mrgstr_index + 1 enddo - !index_lfrac = mct_aVect_indexRA(fractions_g,"lfrac") - !do i = 1, num_flux_fields + do i = 1, num_tf_state_fields + call seq_flds_getField(field, i, seq_flds_x2g_tf_states_from_ocn) + index_o2x = mct_aVect_indexRA(o2x_g, trim(field)) + index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) + + if (first_time) then + mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & + ' = o2x%'//trim(field) + end if + + do n = 1, lsize + x2g_g%rAttr(index_x2g,n) = o2x_g%rAttr(index_o2x,n) + end do + + mrgstr_index = mrgstr_index + 1 + enddo - ! call seq_flds_getField(field, i, seq_flds_x2g_fluxes_from_lnd) - ! index_l2x = mct_aVect_indexRA(l2x_g, trim(field)) - ! index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) + do i = 1, num_flux_fields + call seq_flds_getField(field, i, seq_flds_x2g_shelf_fluxes_from_ocn) + index_o2x = mct_aVect_indexRA(o2x_g, trim(field)) + index_x2g = mct_aVect_indexRA(x2g_g, trim(field)) - ! if (trim(field) == qice_fieldname) then + if (first_time) then + mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & + ' = o2x%'//trim(field) + end if - ! if (first_time) then - ! mrgstr(mrgstr_index) = subname//'x2g%'//trim(field)//' =' // & - ! ' = l2x%'//trim(field) - ! end if + do n = 1, lsize + x2g_g%rAttr(index_x2g,n) = o2x_g%rAttr(index_o2x,n) + end do - ! ! treat qice as if it were a state variable, with a simple copy. - ! do n = 1, lsize - ! x2g_g%rAttr(index_x2g,n) = l2x_g%rAttr(index_l2x,n) - ! end do + mrgstr_index = mrgstr_index + 1 + end do - ! else - ! write(logunit,*) subname,' ERROR: Flux fields other than ', & - ! qice_fieldname, ' currently are not handled in lnd2glc remapping.' - ! write(logunit,*) '(Attempt to handle flux field <', trim(field), '>.)' - ! write(logunit,*) 'Substantial thought is needed to determine how to remap other fluxes' - ! write(logunit,*) 'in a smooth, conservative manner.' - ! call shr_sys_abort(subname//& - ! ' ERROR: Flux fields other than qice currently are not handled in lnd2glc remapping.') - ! endif ! qice_fieldname + if (index_x2g_Fogx_qicehi > 0 .and. index_x2g_Foxo_ismh > 0 .and. index_x2g_Foxo_q_li > 0) then + if (first_time) then + mrgstr(mrgstr_index) = subname//'x2g%Fogx_qicehi = x2g%Foxo_ismh + x2g%Foxo_q_li' + end if - ! mrgstr_index = mrgstr_index + 1 + do n = 1, lsize + x2g_g%rAttr(index_x2g_Fogx_qicehi, n) = x2g_g%rAttr(index_x2g_Foxo_ismh, n) + & + x2g_g%rAttr(index_x2g_Foxo_q_li, n) + end do - !end do + mrgstr_index = mrgstr_index + 1 + endif if (first_time) then if (iamroot) then @@ -862,6 +851,8 @@ subroutine prep_glc_calc_o2x_gx(ocn_c2_glctf, ocn_c2_glcshelf, timer) if (ocn_c2_glcshelf) then call seq_map_map(mapper_So2g_shelf, o2x_ox, o2x_gx(eoi), & fldlist=seq_flds_x2g_shelf_states_from_ocn,norm=.true.) + call seq_map_map(mapper_Fo2g_shelf, o2x_ox, o2x_gx(eoi), & + fldlist=seq_flds_x2g_shelf_fluxes_from_ocn,norm=.true.) end if enddo @@ -980,136 +971,6 @@ end subroutine prep_glc_map_one_state_field_lnd2glc !================================================================================================ - subroutine prep_glc_calculate_subshelf_boundary_fluxes - - !--------------------------------------------------------------- - ! Description - ! On the ice sheet grid, calculate shelf boundary fluxes - - use shr_const_mod , only: SHR_CONST_KAPPA_LAND_ICE - - ! Local Variables - - integer :: gsize, n, egi - type(mct_aVect), pointer :: o2x_ox ! Ocn export, ocn grid, cpl pes - type(mct_aVect), pointer :: x2g_gx ! Glc import, glc grid, cpl pes - type(mct_aVect), pointer :: g2x_gx ! Glc import, glc grid, cpl pes - - integer :: index_x2g_So_blt - integer :: index_x2g_So_bls - integer :: index_x2g_So_htv - integer :: index_x2g_So_stv - integer :: index_x2g_So_rhoeff - integer :: index_g2x_Sg_tbot - integer :: index_g2x_Sg_dztbot - integer :: index_g2x_Sg_lithop - integer :: index_g2x_Sg_icemask_floating - - integer :: index_g2x_Sg_blis - integer :: index_g2x_Sg_blit - integer :: index_g2x_Fogx_qiceho - integer :: index_g2x_Fogx_qicelo - integer :: index_x2g_Fogx_qiceli - integer :: index_x2g_Fogx_qicehi - - character(*), parameter :: subname = '(prep_glc_calculate_subshelf_boundary_fluxes)' - !--------------------------------------------------------------- - - if (.not.(glc_present)) return - - do egi = 1,num_inst_glc - - o2x_ox => component_get_c2x_cx(ocn(egi)) - g2x_gx => component_get_c2x_cx(glc(egi)) - x2g_gx => component_get_x2c_cx(glc(egi)) - - !Remap relevant ocean variables to ice sheet grid. - !Done here instead of in glc-frequency mapping so it happens within ocean coupling interval. - ! Also could map o2x_ox->o2x_gx(1) but using x2g_gx as destination allows us to see - ! these fields on the GLC grid of the coupler history file, which helps with debugging. - call seq_map_map(mapper_So2g_shelf, o2x_ox, x2g_gx, & - fldlist=seq_flds_x2g_shelf_states_from_ocn,norm=.true.) - - ! inputs to melt flux calculation - index_x2g_So_blt = mct_avect_indexra(x2g_gx,'So_blt',perrwith='quiet') - index_x2g_So_bls = mct_avect_indexra(x2g_gx,'So_bls',perrwith='quiet') - index_x2g_So_htv = mct_avect_indexra(x2g_gx,'So_htv',perrwith='quiet') - index_x2g_So_stv = mct_avect_indexra(x2g_gx,'So_stv',perrwith='quiet') - index_x2g_So_rhoeff = mct_avect_indexra(x2g_gx,'So_rhoeff',perrwith='quiet') - - index_g2x_Sg_tbot = mct_avect_indexra(g2x_gx,'Sg_tbot',perrwith='quiet') - index_g2x_Sg_dztbot = mct_avect_indexra(g2x_gx,'Sg_dztbot',perrwith='quiet') - index_g2x_Sg_lithop = mct_avect_indexra(g2x_gx,'Sg_lithop',perrwith='quiet') - index_g2x_Sg_icemask_floating = mct_avect_indexra(g2x_gx,'Sg_icemask_floating',perrwith='quiet') - - ! outputs to melt flux calculation - index_g2x_Sg_blis = mct_avect_indexra(g2x_gx,'Sg_blis',perrwith='quiet') - index_g2x_Sg_blit = mct_avect_indexra(g2x_gx,'Sg_blit',perrwith='quiet') - index_g2x_Fogx_qiceho = mct_avect_indexra(g2x_gx,'Fogx_qiceho',perrwith='quiet') - index_g2x_Fogx_qicelo = mct_avect_indexra(g2x_gx,'Fogx_qicelo',perrwith='quiet') - index_x2g_Fogx_qiceli = mct_avect_indexra(x2g_gx,'Fogx_qiceli',perrwith='quiet') - index_x2g_Fogx_qicehi = mct_avect_indexra(x2g_gx,'Fogx_qicehi',perrwith='quiet') - - gsize = mct_aVect_lsize(g2x_gx) - - do n=1,gsize - !Extract glc and ocn-sourced coupler fields used as input to compute_melt_fluxes to local arrays... - - ! Fields from the ocean, now on the GLC grid - oceanTemperature(n) = x2g_gx%rAttr(index_x2g_So_blt,n) - oceanSalinity(n) = x2g_gx%rAttr(index_x2g_So_bls,n) - oceanHeatTransferVelocity(n) = x2g_gx%rAttr(index_x2g_So_htv,n) - oceanSaltTransferVelocity(n) = x2g_gx%rAttr(index_x2g_So_stv,n) - - ! Fields from the ice sheet model (still on the GLC grid) - iceTemperature(n) = g2x_gx%rAttr(index_g2x_Sg_tbot,n) - iceTemperatureDistance(n) = g2x_gx%rAttr(index_g2x_Sg_dztbot,n) - interfacePressure(n) = g2x_gx%rAttr(index_g2x_Sg_lithop,n) - iceFloatingMask(n) = g2x_gx%rAttr(index_g2x_Sg_icemask_floating,n) - - !... initialize local compute_melt_fluxes output arrays... - outInterfaceSalinity(n) = 0.0_r8 - outInterfaceTemperature(n) = 0.0_r8 - outFreshwaterFlux(n) = 0.0_r8 - outOceanHeatFlux(n) = 0.0_r8 - outIceHeatFlux(n) = 0.0_r8 - end do - - !...calculate fluxes... - call compute_melt_fluxes(oceanTemperature=oceanTemperature,& - oceanSalinity=oceanSalinity,& - oceanHeatTransferVelocity=oceanHeatTransferVelocity,& - oceanSaltTransferVelocity=oceanSaltTransferVelocity,& - interfacePressure=interfacePressure,& - iceTemperature=iceTemperature,& - iceTemperatureDistance=iceTemperatureDistance, & - iceFloatingMask=iceFloatingMask, & - outInterfaceSalinity=outInterfaceSalinity,& - outInterfaceTemperature=outInterfaceTemperature,& - outFreshwaterFlux=outFreshwaterFlux,& - outOceanHeatFlux=outOceanHeatFlux,& - outIceHeatFlux=outIceHeatFlux,& - gsize=gsize) - - !...and assign fluxes to glc and ocn-directed coupler fields - do n=1,gsize - !Assign outputs from compute_melt_fluxes back into coupler attributes - g2x_gx%rAttr(index_g2x_Sg_blis,n) = outInterfaceSalinity(n) !to ocean - g2x_gx%rAttr(index_g2x_Sg_blit,n) = outInterfaceTemperature(n) !to ocean - g2x_gx%rAttr(index_g2x_Fogx_qiceho,n) = outOceanHeatFlux(n) !to ocean - g2x_gx%rAttr(index_g2x_Fogx_qicelo,n)= outFreshwaterFlux(n) !to ocean - x2g_gx%rAttr(index_x2g_Fogx_qicehi,n) = outIceHeatFlux(n) !to ice sheet - x2g_gx%rAttr(index_x2g_Fogx_qiceli,n) = -1.0_r8 * outFreshwaterFlux(n) !to ice sheet - end do - - !Note: remap ocean-side outputs back onto ocean grid done in call to prep_ocn_shelf_calc_g2x_ox - - end do ! loop over GLC instances - - end subroutine prep_glc_calculate_subshelf_boundary_fluxes - - !================================================================================================ - subroutine prep_glc_zero_fields() !--------------------------------------------------------------- @@ -1650,165 +1511,4 @@ function prep_glc_get_mapper_Fo2g_shelf() prep_glc_get_mapper_Fo2g_shelf=> mapper_Fo2g_shelf end function prep_glc_get_mapper_Fo2g_shelf -!*********************************************************************** -! -! routine compute_melt_fluxes -! -!> \brief Computes ocean and ice melt fluxes, etc. -!> \author Xylar Asay-Davis -!> \date 3/27/2015 -!> This routine computes melt fluxes (melt rate, temperature fluxes -!> into the ice and the ocean, and salt flux) as well as the interface -!> temperature and salinity. This routine expects an ice temperature -!> in the bottom layer of ice and ocean temperature and salinity in -!> the top ocean layer as well as the pressure at the ice/ocean interface. -!> -!> The ocean heat and salt transfer velocities are determined based on -!> observations of turbulent mixing rates in the under-ice boundary layer. -!> They should be the product of the friction velocity and a (possibly -!> spatially variable) non-dimenional transfer coefficient. -!> -!> The iceTemperatureDistance is the distance between the location -!> where the iceTemperature is supplied and the ice-ocean interface, -!> used to compute a temperature gradient. The ice thermal conductivity, -!> SHR_CONST_KAPPA_LAND_ICE, is zero for the freezing solution from Holland and Jenkins -!> (1999) in which the ice is purely insulating. -! -!----------------------------------------------------------------------- - - subroutine compute_melt_fluxes( & - oceanTemperature, & - oceanSalinity, & - oceanHeatTransferVelocity, & - oceanSaltTransferVelocity, & - interfacePressure, & - iceTemperature, & - iceTemperatureDistance, & - iceFloatingMask, & - outInterfaceSalinity, & - outInterfaceTemperature, & - outFreshwaterFlux, & - outOceanHeatFlux, & - outIceHeatFlux, & - gsize) - - use shr_const_mod, only: SHR_CONST_CPICE, & - SHR_CONST_CPSW, & - SHR_CONST_LATICE, & - SHR_CONST_RHOICE, & - SHR_CONST_RHOSW, & - SHR_CONST_DTF_DP, & - SHR_CONST_DTF_DS, & - SHR_CONST_DTF_DPDS, & - SHR_CONST_TF0, & - SHR_CONST_KAPPA_LAND_ICE - - !----------------------------------------------------------------- - ! - ! input variables - ! - !----------------------------------------------------------------- - - real (kind=r8), dimension(:), intent(in) :: & - oceanTemperature, & !< Input: ocean temperature in top layer - oceanSalinity, & !< Input: ocean salinity in top layer - oceanHeatTransferVelocity, & !< Input: ocean heat transfer velocity - oceanSaltTransferVelocity, & !< Input: ocean salt transfer velocity - interfacePressure, & !< Input: pressure at the ice-ocean interface - iceTemperature, & !< Input: ice temperature in bottom layer - iceTemperatureDistance !< Input: distance to ice temperature from ice-ocean interface - integer, dimension(:), intent(in) :: & - iceFloatingMask !< Input: mask of cells that contain floating ice - - integer, intent(in) :: gsize !< Input: number of values in each array - - !----------------------------------------------------------------- - ! - ! output variables - ! - !----------------------------------------------------------------- - - real (kind=r8), dimension(:), intent(out) :: & - outInterfaceSalinity, & !< Output: ocean salinity at the interface - outInterfaceTemperature, & !< Output: ice/ocean temperature at the interface - outFreshwaterFlux, & !< Output: ocean thickness flux (melt rate) - outOceanHeatFlux, & !< Output: the temperature flux into the ocean - outIceHeatFlux !< Output: the temperature flux into the ice - - !----------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------- - - real (kind=r8) :: T0, transferVelocityRatio, Tlatent, nu, a, b, c, eta, & - iceHeatFluxCoeff, iceDeltaT, dTf_dS - integer :: n - character(*), parameter :: subname = '(compute_melt_fluxes)' - - real (kind=r8), parameter :: minInterfaceSalinity = 0.001_r8 - - real (kind=r8), parameter :: referencePressure = 0.0_r8 ! Using reference pressure of 0 - - real (kind=r8) :: pressureOffset - - Tlatent = SHR_CONST_LATICE/SHR_CONST_CPSW - do n = 1, gsize - if (iceFloatingMask(n) == 0) cycle ! Only calculate on floating cells - - if (oceanHeatTransferVelocity(n) == 0.0_r8) then - write(logunit,*) 'compute_melt_fluxes ERROR: oceanHeatTransferVelocity value of 0 causes divide by 0 at index ', n - call shr_sys_abort('compute_melt_fluxes ERROR: oceanHeatTransferVelocity value of 0 causes divide by 0') - end if - - iceHeatFluxCoeff = SHR_CONST_RHOICE*SHR_CONST_CPICE*SHR_CONST_KAPPA_LAND_ICE/iceTemperatureDistance(n) - nu = iceHeatFluxCoeff/(SHR_CONST_RHOSW*SHR_CONST_CPSW*oceanHeatTransferVelocity(n)) - pressureOffset = max(interfacePressure(n) - referencePressure, 0.0_r8) - T0 = SHR_CONST_TF0 + SHR_CONST_DTF_DP * pressureOffset - !Note: These two terms for T0 are not needed because we are evaluating at salinity=0: - !+ SHR_CONST_DTF_DS * oceanSalinity(n) + SHR_CONST_DTF_DPDS * pressureOffset * oceanSalinity(n) - iceDeltaT = T0 - iceTemperature(n) - dTf_dS = SHR_CONST_DTF_DS + SHR_CONST_DTF_DPDS * pressureOffset - - transferVelocityRatio = oceanSaltTransferVelocity(n)/oceanHeatTransferVelocity(n) - - a = -1.0_r8 * dTf_dS * (1.0_r8 + nu) - b = transferVelocityRatio*Tlatent - nu*iceDeltaT + oceanTemperature(n) - T0 - c = -transferVelocityRatio*Tlatent*max(oceanSalinity(n), 0.0_r8) - ! a is non-negative; c is strictly non-positive so we never get imaginary roots. - ! Since a can be zero, we need a solution of the quadratic equation for 1/Si instead of Si. - ! Following: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf - ! Since a and -c are are non-negative, the term in the square root is also always >= |b|. - ! In all reasonable cases, b will be strictly positive, since transferVelocityRatio*Tlatent ~ 2 C, - ! T0 ~ -1.8 C and oceanTemperature should never be able to get below about -3 C - ! As long as either b or both a and c are greater than zero, the strictly non-negative root is - outInterfaceSalinity(n) = max(-(2.0_r8*c)/(b + sqrt(b**2 - 4.0_r8*a*c)), minInterfaceSalinity) - - outInterfaceTemperature(n) = dTf_dS*outInterfaceSalinity(n)+T0 - - outFreshwaterFlux(n) = SHR_CONST_RHOSW*oceanSaltTransferVelocity(n) & - * (oceanSalinity(n)/outInterfaceSalinity(n) - 1.0_r8) - - ! According to Jenkins et al. (2001), the temperature fluxes into the ocean are: - ! 1. the advection of meltwater into the top layer (or removal for freezing) - ! 2. the turbulent transfer of heat across the boundary layer, based on the termal driving - outOceanHeatFlux(n) = SHR_CONST_CPSW*(outFreshwaterFlux(n)*outInterfaceTemperature(n) & - - SHR_CONST_RHOSW*oceanHeatTransferVelocity(n)*(oceanTemperature(n)-outInterfaceTemperature(n))) - - ! the temperature fluxes into the ice are: - ! 1. the advection of ice at the interface temperature out of the domain due to melting - ! (or in due to freezing) - ! 2. the diffusion (if any) of heat into the ice, based on temperature difference between - ! the reference point in the ice (either the surface or the middle of the bottom layer) - ! and the interface - outIceHeatFlux(n) = -SHR_CONST_CPICE*outFreshwaterFlux(n)*outInterfaceTemperature(n) - - outIceHeatFlux(n) = outIceHeatFlux(n) & - - iceHeatFluxCoeff*(iceTemperature(n) - outInterfaceTemperature(n)) - - end do - - !-------------------------------------------------------------------- - end subroutine compute_melt_fluxes - end module prep_glc_mod diff --git a/driver-mct/main/prep_ocn_mod.F90 b/driver-mct/main/prep_ocn_mod.F90 index 32fc58992372..959eb9559b11 100644 --- a/driver-mct/main/prep_ocn_mod.F90 +++ b/driver-mct/main/prep_ocn_mod.F90 @@ -543,7 +543,6 @@ end subroutine prep_ocn_mrg subroutine prep_ocn_merge( flux_epbalfact, a2x_o, i2x_o, r2x_o, w2x_o, g2x_o, xao_o, & fractions_o, x2o_o ) - use prep_glc_mod, only: prep_glc_calculate_subshelf_boundary_fluxes use seq_flds_mod, only: wav_ocn_coup !----------------------------------------------------------------------- @@ -1312,7 +1311,7 @@ subroutine prep_ocn_calc_a2x_ox(timer) #ifdef COMPARE_TO_NUOPC call seq_map_mapvect(mapper_Va2o, vect_map, a2x_ax, a2x_ox(eai), 'Sa_u', 'Sa_v', norm=.true.) -#else +#else !--- tcx the norm should be true below, it's false for bfb backwards compatability call seq_map_mapvect(mapper_Va2o, vect_map, a2x_ax, a2x_ox(eai), 'Sa_u', 'Sa_v', norm=.false.) #endif diff --git a/driver-mct/main/seq_diag_mct.F90 b/driver-mct/main/seq_diag_mct.F90 index 3c3872d48e50..d03fbe6ff3e1 100644 --- a/driver-mct/main/seq_diag_mct.F90 +++ b/driver-mct/main/seq_diag_mct.F90 @@ -49,7 +49,7 @@ module seq_diag_mct use seq_diagBGC_mct, only : seq_diagBGC_preprint_mct, seq_diagBGC_print_mct use prep_glc_mod, only : prep_glc_get_l2gacc_lx_cnt_avg - use glc_elevclass_mod, only: glc_get_num_elevation_classes + use glc_elevclass_mod, only: glc_get_num_elevation_classes implicit none save @@ -82,17 +82,6 @@ module seq_diag_mct !---------------------------------------------------------------------------- !----- local constants ----- - real(r8),parameter :: HFLXtoWFLX = & ! water flux implied by latent heat of fusion - & - (shr_const_ocn_ref_sal-shr_const_ice_ref_sal) / & - & (shr_const_ocn_ref_sal*shr_const_latice) - - real(r8),parameter :: SFLXtoWFLX = & ! water flux implied by salt flux - ! WFLX (kg/m^2s) = -SFLX (kg/m^2s) - ! / ocn_ref_sal (psu) (34.7g/kg) - ! / 1.e-3 kg/g - -1._r8/(shr_const_ocn_ref_sal*1.e-3_r8) - - !--- C for component --- !--- "r" is recieve in the coupler, "s" is send from the coupler @@ -141,14 +130,14 @@ module seq_diag_mct integer(in),parameter :: f_hlatf = 8 ! heat : latent, fusion, snow integer(in),parameter :: f_hioff = 9 ! heat : latent, fusion, frozen runoff integer(in),parameter :: f_hsen =10 ! heat : sensible - integer(in),parameter :: f_hpolar =11 ! heat : AIS imbalance + integer(in),parameter :: f_hpolar =11 ! heat : non-latent shelf exchange integer(in),parameter :: f_hh2ot =12 ! heat : water temperature 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_wpolar =18 ! water: frozen-freshwater shelf exchange 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 @@ -269,7 +258,7 @@ module seq_diag_mct integer :: index_l2x_Flrl_wslake integer :: index_x2l_Sg_icemask - integer, allocatable :: index_l2x_Flgl_qice(:) + integer, allocatable :: index_l2x_Flgl_qice(:) integer, allocatable :: index_x2l_Sg_ice_covered(:) integer :: index_x2l_Faxa_lwdn @@ -347,7 +336,7 @@ module seq_diag_mct integer :: index_g2x_Fogg_rofi integer :: index_g2x_Figg_rofi - integer :: index_x2g_Flgl_qice + integer :: index_x2g_Flgl_qice integer :: index_g2x_Sg_icemask integer :: index_x2o_Foxx_rofl_16O @@ -889,12 +878,12 @@ 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 ! for summing fluxes over no. of elev. classes + 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 - integer(in) :: num + character(len= 2) :: cnum + integer(in) :: num !----- formats ----- character(*),parameter :: subName = '(seq_diag_lnd_mct) ' @@ -916,7 +905,7 @@ 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 + ! get number of elevation classes and allocate relevant sets of indices glc_nec = glc_get_num_elevation_classes() if (glc_nec.ge.1) then if (first_time) then @@ -924,7 +913,7 @@ subroutine seq_diag_lnd_mct( lnd, frac_l, infodata, do_l2x, do_x2l) allocate(index_x2l_Sg_ice_covered(0:glc_nec)) end if end if - + if (present(do_l2x)) then if (first_time) then index_l2x_Fall_swnet = mct_aVect_indexRA(l2x_l,'Fall_swnet') @@ -989,9 +978,9 @@ 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 + 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) * effective_area end do @@ -1031,7 +1020,7 @@ 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 ! contribution from land ice calving currently zero + 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 ! Nneeded? Not sure if / when these should be deallocated @@ -1321,7 +1310,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 (may not be used / needed 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 @@ -1391,7 +1380,7 @@ 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_) - l2gacc_lx_cnt_avg = prep_glc_get_l2gacc_lx_cnt_avg() ! counter for how many times SMB flux accumulation has occured + 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(x2g_g) @@ -1403,7 +1392,7 @@ subroutine seq_diag_glc_mct( glc, frac_g, infodata, do_x2g, do_g2x ) 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 + budg_dataL(f_hgsmb,ic,ip) = budg_dataL(f_wgsmb,ic,ip)*shr_const_latice end if ! end do fields from coupler to glc (x2g_) @@ -1503,6 +1492,8 @@ subroutine seq_diag_ocn_mct( ocn, xao_o, frac_o, infodata, do_o2x, do_x2o, do_xa nf = f_hfrz; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*max(0.0_r8,o2x_o%rAttr(index_o2x_Fioo_q,n)) nf = f_hh2ot; budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + (ca_o+ca_i)*o2x_o%rAttr(index_o2x_Faoo_h2otemp,n) if (flds_polar) then + ! Shelf coupling budgets track frozen-freshwater mass exchange in wpolar + ! and only non-latent shelf heat exchange in hpolar. nf = f_wpolar;budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - ca_c*o2x_o%rAttr(index_o2x_Foxo_frazil_li,n) nf = f_wpolar;budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) + ca_c*o2x_o%rAttr(index_o2x_Foxo_ismw,n) nf = f_wpolar;budg_dataL(nf,ic,ip) = budg_dataL(nf,ic,ip) - (ca_o+ca_i)*o2x_o%rAttr(index_o2x_Foxo_rrofl,n) @@ -1973,7 +1964,7 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, & if (plev > 0) then ! ---- doprint ---- doprint ---- doprint ---- - + if (.not.sumdone) then if (do_bgc_budg) then @@ -2340,7 +2331,7 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, & ! ---- doprint ---- doprint ---- doprint ---- if (do_bgc_budg) then - call seq_diagBGC_print_mct(EClock, ip, plev) + call seq_diagBGC_print_mct(EClock, ip, plev) endif endif ! plev > 0 diff --git a/driver-mct/shr/seq_flds_mod.F90 b/driver-mct/shr/seq_flds_mod.F90 index f9960c1761ec..e151056756a2 100644 --- a/driver-mct/shr/seq_flds_mod.F90 +++ b/driver-mct/shr/seq_flds_mod.F90 @@ -221,6 +221,7 @@ module seq_flds_mod character(CXX) :: seq_flds_x2g_states character(CXX) :: seq_flds_x2g_states_from_lnd character(CXX) :: seq_flds_x2g_shelf_states_from_ocn + character(CXX) :: seq_flds_x2g_shelf_fluxes_from_ocn character(CXX) :: seq_flds_x2g_tf_states_from_ocn character(CXX) :: seq_flds_x2g_fluxes character(CXX) :: seq_flds_x2g_fluxes_from_lnd @@ -364,6 +365,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) character(CXX) :: x2g_states = '' character(CXX) :: x2g_states_from_lnd = '' character(CXX) :: x2g_shelf_states_from_ocn = '' + character(CXX) :: x2g_shelf_fluxes_from_ocn = '' character(CXX) :: x2g_tf_states_from_ocn = '' character(CXX) :: x2g_fluxes = '' character(CXX) :: x2g_fluxes_from_lnd = '' @@ -821,7 +823,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'kg m-3' attname = 'Sa_dens' call metadata_set(attname, longname, stdname, units) - + ! UoverN for use by topounits call seq_flds_add(a2x_states,"Sa_uovern") call seq_flds_add(x2l_states,"Sa_uovern") @@ -1644,16 +1646,20 @@ subroutine seq_flds_set(nmlfile, ID, infodata) if (flds_polar) then - ! Ocean Land ice freeze potential + ! Non-latent frazil enthalpy flux exported from ocean to land ice call seq_flds_add(o2x_fluxes,"Foxo_q_li") - longname = 'Ocean land ice freeze potential' - stdname = 'ice_shelf_cavity_ice_heat_flux' + call seq_flds_add(x2g_fluxes,"Foxo_q_li") + call seq_flds_add(x2g_shelf_fluxes_from_ocn,"Foxo_q_li") + longname = 'Non-latent frazil enthalpy flux at the land-ice interface' + stdname = 'ocean_land_ice_frazil_enthalpy_flux' units = 'W m-2' attname = 'Foxo_q_li' call metadata_set(attname, longname, stdname, units) ! Ocean land ice frazil production call seq_flds_add(o2x_fluxes,"Foxo_frazil_li") + call seq_flds_add(x2g_fluxes,"Foxo_frazil_li") + call seq_flds_add(x2g_shelf_fluxes_from_ocn,"Foxo_frazil_li") longname = 'Ocean land ice frazil production' stdname = 'ocean_land_ice_frazil_ice_production' units = 'kg m-2 s-1' @@ -1662,16 +1668,20 @@ subroutine seq_flds_set(nmlfile, ID, infodata) ! Water flux from ice shelf melt call seq_flds_add(o2x_fluxes,"Foxo_ismw") + call seq_flds_add(x2g_fluxes,"Foxo_ismw") + call seq_flds_add(x2g_shelf_fluxes_from_ocn,"Foxo_ismw") longname = 'Water flux due to basal melting of ice shelves' stdname = 'basal_iceshelf_melt_flux' units = 'kg m-2 s-1' attname = 'Foxo_ismw' call metadata_set(attname, longname, stdname, units) - ! Heat flux from ice shelf melt + ! Non-latent heat flux from ice shelf melt/freeze call seq_flds_add(o2x_fluxes,"Foxo_ismh") - longname = 'Heat flux due to basal melting of ice shelves' - stdname = 'basal_iceshelf_heat_flux' + call seq_flds_add(x2g_fluxes,"Foxo_ismh") + call seq_flds_add(x2g_shelf_fluxes_from_ocn,"Foxo_ismh") + longname = 'Non-latent interfacial heat flux due to basal melting or freezing of ice shelves' + stdname = 'basal_iceshelf_non_latent_heat_flux' units = 'W m-2' attname = 'Foxo_ismh' call metadata_set(attname, longname, stdname, units) @@ -2236,7 +2246,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) endif !------------------------------ - ! ice<->wav exchange + ! ice<->wav exchange !------------------------------ ! Sea ice thickness @@ -2371,7 +2381,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = ' ' attname = 'coszen_str' call metadata_set(attname, longname, stdname, units) - + if (rof_sed) then call seq_flds_add(l2x_fluxes,'Flrl_rofmud') call seq_flds_add(l2x_fluxes_to_rof,'Flrl_rofmud') @@ -2451,7 +2461,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'kg m-2 s-1' attname = 'Flrr_supply' call metadata_set(attname, longname, stdname, units) - + call seq_flds_add(r2x_fluxes,'Flrr_deficit') call seq_flds_add(x2l_fluxes,'Flrr_deficit') longname = 'River model supply deficit' @@ -2726,7 +2736,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'deg' attname = 'Sw_Dp' call metadata_set(attname, longname, stdname, units) - + call seq_flds_add(w2x_fluxes,'Faww_Tawx') call seq_flds_add(x2o_fluxes,'Faww_Tawx') longname = 'Zonal wave supported stress' @@ -2792,7 +2802,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = '' attname = 'Sw_Ustar' call metadata_set(attname, longname, stdname, units) - + call seq_flds_add(w2x_states,'Sw_Z0') if (wav_ocn_coup == 'twoway') call seq_flds_add(x2o_states,'Sw_Z0') longname = 'Surface Roughness Length based on wave state' @@ -2805,12 +2815,12 @@ subroutine seq_flds_set(nmlfile, ID, infodata) !---------------------------- ! lnd->iac, iac->lnd, iac->atm !---------------------------- - + ! lnd/iac coupling needs one field in each class per pft ! Note that this ends up as 17*4=68 coupled fields... ! also send harvest fraction from iac to land ! use the same loop and index string - + ! these two variables are hardcoded above because there is not an appropriate ! namelist to put them in: iac_npft and iac_nharvest @@ -2829,7 +2839,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'Sl_hr_pft' // pftstr attname = trim(attname) call metadata_set(attname, longname, stdname, units) - + if(add_iac_to_cplstate)call seq_flds_add(l2x_states,'Sl_npp_pft' // pftstr) call seq_flds_add(x2z_states,'Sl_npp_pft' // pftstr) longname = 'Net primary production for pft ' // pftstr @@ -2837,7 +2847,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'gC/m^2/s' attname = 'Sl_npp_pft' // pftstr call metadata_set(attname, longname, stdname, units) - + ! Review if(add_iac_to_cplstate)call seq_flds_add(l2x_states,'Sl_pftwgt_pft' //pftstr) call seq_flds_add(x2z_states,'Sl_pftwgt_pft' //pftstr) @@ -2847,7 +2857,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'Sl_pftwgt_pft' //pftstr call metadata_set(attname, longname, stdname, units) - ! iac->lnd + ! iac->lnd ! This is pft for beginning of model year + 1 ! ts wonders if landfrac should go as well - just to @@ -2871,8 +2881,8 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'percent' attname = 'Sz_pct_pft_prev' //pftstr attname = trim(attname) - call metadata_set(attname, longname, stdname, units) - + call metadata_set(attname, longname, stdname, units) + ! send the harvest data also, these are for model year if (i <= iac_nharvest) then call seq_flds_add(z2x_states,trim('Sz_harvest_frac' //pftstr)) @@ -2886,7 +2896,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call metadata_set(attname, longname, stdname, units) end if - end do + end do ! iac->atm flux. ! Monthly values of surface, low alt, high alt co2 fluxes, so we ! loop over 36 total fields. @@ -3217,54 +3227,24 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call set_glc_elevclass_field(name, attname, longname, stdname, units, x2l_fluxes_from_glc, & additional_list = .true.) - name = 'So_blt' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf boundary layer ocean temperature' - stdname = 'Ice_shelf_boundary_layer_ocean_temperature' - units = 'C' - attname = 'So_blt' - call metadata_set(attname, longname, stdname, units) - - name = 'So_bls' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf boundary layer ocean salinity' - stdname = 'Ice_shelf_boundary_layer_ocean_salinity' - units = 'psu' - attname = 'So_bls' - call metadata_set(attname, longname, stdname, units) - - name = 'So_htv' + name = 'So_tfrz_isf' call seq_flds_add(o2x_states,trim(name)) call seq_flds_add(x2g_states,trim(name)) call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf ocean heat transfer velocity' - stdname = 'Ice_shelf_ocean_heat_transfer_velocity' - units = 'm/s' - attname = 'So_htv' - call metadata_set(attname, longname, stdname, units) - - name = 'So_stv' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf ocean salinity transfer velocity' - stdname = 'Ice_shelf_ocean_salinity_transfer_velocity' - units = 'm/s' - attname = 'So_stv' + longname = 'Ice shelf-ocean freezing temperature' + stdname = 'Ice_shelf_ocean_freezing_temperature' + units = 'K' + attname = name call metadata_set(attname, longname, stdname, units) - name = 'So_rhoeff' + name = 'So_liflfrac' call seq_flds_add(o2x_states,trim(name)) call seq_flds_add(x2g_states,trim(name)) call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ocean effective pressure' - stdname = 'Ocean_effective_pressure' - units = 'Pa' - attname = 'So_rhoeff' + longname = 'Floating ice shelf area fraction' + stdname = 'floating_ice_shelf_area_fraction' + units = '1' + attname = name call metadata_set(attname, longname, stdname, units) if ((flds_tf) .and. (glc_nzoc > 0)) then @@ -3296,42 +3276,6 @@ subroutine seq_flds_set(nmlfile, ID, infodata) additional_list = .true.) end if - name = 'Fogx_qicelo' - call seq_flds_add(g2x_fluxes,trim(name)) - call seq_flds_add(x2o_fluxes,trim(name)) - longname = 'Subshelf liquid flux for ocean' - stdname = 'Subshelf_liquid_flux_for_ocean' - units = 'kg m-2 s-1' - attname = 'Fogx_qicelo' - call metadata_set(attname, longname, stdname, units) - - name = 'Fogx_qiceho' - call seq_flds_add(g2x_fluxes,trim(name)) - call seq_flds_add(x2o_fluxes,trim(name)) - longname = 'Subshelf heat flux for the ocean' - stdname = 'Subshelf_heat_flux_for_the_ocean' - units = 'W m-2' - attname = 'Fogx_qiceho' - call metadata_set(attname, longname, stdname, units) - - name = 'Sg_blit' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Boundary layer interface temperature for ocean' - stdname = 'Boundary_layer_interface_temperature_for_ocean' - units = 'C' - attname = 'Sg_blit' - call metadata_set(attname, longname, stdname, units) - - name = 'Sg_blis' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Boundary layer interface salinity for ocean' - stdname = 'Boundary_layer_interface_salinity_for_ocean' - units = 'psu' - attname = 'Sg_blis' - call metadata_set(attname, longname, stdname, units) - name = 'Sg_lithop' call seq_flds_add(g2x_states,trim(name)) call seq_flds_add(x2o_states,trim(name)) @@ -3359,32 +3303,6 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'Sg_icemask_floating' call metadata_set(attname, longname, stdname, units) - name = 'Sg_tbot' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Bottom layer ice temperature' - stdname = 'Bottom_layer_ice_temperature' - units = 'C' - attname = 'Sg_tbot' - call metadata_set(attname, longname, stdname, units) - - name = 'Sg_dztbot' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Bottom layer ice layer half thickness' - stdname = 'Bottom_layer_ice_layer_half_thickness' - units = 'm' - attname = 'Sg_dztbot' - call metadata_set(attname, longname, stdname, units) - - name = 'Fogx_qiceli' - call seq_flds_add(x2g_fluxes,trim(name)) - longname = 'Subshelf mass flux for ice sheet' - stdname = 'Subshelf_mass_flux_for_ice_sheet' - units = 'kg m-2 s-1' - attname = 'Fogx_qiceli' - call metadata_set(attname, longname, stdname, units) - name = 'Fogx_qicehi' call seq_flds_add(x2g_fluxes,trim(name)) longname = 'Subshelf heat flux for ice sheet' @@ -4166,7 +4084,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) !----------------------------------------------------------------------------- ! Read namelist for FAN NH3 emissions - ! If specified, the NH3 surface emission is sent to CAM. + ! If specified, the NH3 surface emission is sent to CAM. !----------------------------------------------------------------------------- call shr_fan_readnl(nlfilename='drv_flds_in', ID=ID, fan_fields=fan_fields, have_fields=fan_have_fields) @@ -4247,6 +4165,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) seq_flds_x2g_states = trim(x2g_states) seq_flds_x2g_states_from_lnd = trim(x2g_states_from_lnd) seq_flds_x2g_shelf_states_from_ocn = trim(x2g_shelf_states_from_ocn) + seq_flds_x2g_shelf_fluxes_from_ocn = trim(x2g_shelf_fluxes_from_ocn) seq_flds_x2g_tf_states_from_ocn = trim(x2g_tf_states_from_ocn) seq_flds_xao_states = trim(xao_states) seq_flds_xao_albedo = trim(xao_albedo) @@ -4321,6 +4240,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) write(logunit,*) subname//': seq_flds_x2g_states_from_lnd= ',trim(seq_flds_x2g_states_from_lnd) write(logunit,*) subname//': seq_flds_l2x_states_to_glc= ',trim(seq_flds_l2x_states_to_glc) write(logunit,*) subname//': seq_flds_x2g_shelf_states_from_ocn= ',trim(seq_flds_x2g_shelf_states_from_ocn) + write(logunit,*) subname//': seq_flds_x2g_shelf_fluxes_from_ocn= ',trim(seq_flds_x2g_shelf_fluxes_from_ocn) write(logunit,*) subname//': seq_flds_x2g_tf_states_from_ocn= ',trim(seq_flds_x2g_tf_states_from_ocn) write(logunit,*) subname//': seq_flds_x2g_fluxes= ',trim(seq_flds_x2g_fluxes) write(logunit,*) subname//': seq_flds_x2g_fluxes_from_lnd= ',trim(seq_flds_x2g_fluxes_from_lnd) diff --git a/driver-moab/main/cime_comp_mod.F90 b/driver-moab/main/cime_comp_mod.F90 index 78939f6f5075..89dfab6c1fd8 100644 --- a/driver-moab/main/cime_comp_mod.F90 +++ b/driver-moab/main/cime_comp_mod.F90 @@ -4395,22 +4395,9 @@ subroutine cime_run_ocnglc_coupling() if (glc_present) then - if (ocn_c2_glcshelf .and. glcshelf_c2_ocn) then - ! the boundary flux calculations done in the coupler require inputs from both GLC and OCN, - ! so they will only be valid if both OCN->GLC and GLC->OCN - - call prep_glc_calc_o2x_gx(timer='CPL:glcprep_ocn2glc') !remap ocean fields to o2x_g at ocean couping interval - - call prep_glc_calculate_subshelf_boundary_fluxes ! this is actual boundary layer flux calculation - !this outputs - !x2g_g/g2x_g, where latter is going - !to ocean, so should get remapped to - !ocean grid in prep_ocn_shelf_calc_g2x_ox - call prep_ocn_shelf_calc_g2x_ox(timer='CPL:glcpost_glcshelf2ocn') - !Map g2x_gx shelf fields that were updated above, to g2x_ox. - !Do this at intrinsic coupling - !frequency - call prep_glc_accum_ocn(timer='CPL:glcprep_accum_ocn') !accum x2g_g fields here into x2g_gacc + ! create o2x_gx for ocn-glc shelf coupling + if (ocn_c2_glctf) then + call prep_glc_calc_o2x_gx(ocn_c2_glctf, ocn_c2_glcshelf, timer='CPL:glcprep_ocn2glc') !remap ocean fields to o2x_g at ocean couping interval endif if (glcshelf_c2_ice) then diff --git a/driver-moab/main/prep_glc_mod.F90 b/driver-moab/main/prep_glc_mod.F90 index 339acaaaac04..333420dfa6a4 100644 --- a/driver-moab/main/prep_glc_mod.F90 +++ b/driver-moab/main/prep_glc_mod.F90 @@ -56,8 +56,6 @@ module prep_glc_mod public :: prep_glc_get_mapper_So2g public :: prep_glc_get_mapper_Fo2g - public :: prep_glc_calculate_subshelf_boundary_fluxes - !-------------------------------------------------------------------------- ! Private interfaces !-------------------------------------------------------------------------- @@ -85,9 +83,9 @@ module prep_glc_mod type(mct_aVect), pointer :: o2x_gx(:) ! Ocn export, glc grid, cpl pes - allocated in driver ! accumulation variables - + type(mct_aVect), pointer :: x2gacc_gx(:) ! Glc export, glc grid, cpl pes - allocated in driver - integer , target :: x2gacc_gx_cnt ! x2gacc_gx: number of time samples accumulated + integer , target :: x2gacc_gx_cnt ! x2gacc_gx: number of time samples accumulated 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 @@ -115,20 +113,6 @@ module prep_glc_mod type(mct_aVect), pointer :: o2gacc_ox(:) ! Ocn export, lnd grid, cpl pes - allocated in driver integer , target :: o2gacc_ox_cnt ! number of time samples accumulated - real(r8), allocatable :: oceanTemperature(:) - real(r8), allocatable :: oceanSalinity(:) - real(r8), allocatable :: oceanHeatTransferVelocity(:) - real(r8), allocatable :: oceanSaltTransferVelocity(:) - real(r8), allocatable :: interfacePressure(:) - real(r8), allocatable :: iceTemperature(:) - real(r8), allocatable :: iceTemperatureDistance(:) - integer, allocatable :: iceFloatingMask(:) - real(r8), allocatable :: outInterfaceSalinity(:) - real(r8), allocatable :: outInterfaceTemperature(:) - real(r8), allocatable :: outFreshwaterFlux(:) - real(r8), allocatable :: outOceanHeatFlux(:) - real(r8), allocatable :: outIceHeatFlux(:) - !================================================================================================ contains @@ -290,23 +274,6 @@ subroutine prep_glc_init(infodata, lnd_c2_glc, ocn_c2_glcshelf) 'seq_maps.rc','ocn2glc_fmapname:','ocn2glc_fmaptype:',samegrid_go, & 'mapper_Fo2g initialization',esmf_map_flag) - !Initialize module-level arrays associated with compute_melt_fluxes - allocate(oceanTemperature(lsize_g)) - allocate(oceanSalinity(lsize_g)) - allocate(oceanHeatTransferVelocity(lsize_g)) - allocate(oceanSaltTransferVelocity(lsize_g)) - allocate(interfacePressure(lsize_g)) - allocate(iceTemperature(lsize_g)) - allocate(iceTemperatureDistance(lsize_g)) - allocate(iceFloatingMask(lsize_g)) - allocate(outInterfaceSalinity(lsize_g)) - allocate(outInterfaceTemperature(lsize_g)) - allocate(outFreshwaterFlux(lsize_g)) - allocate(outOceanHeatFlux(lsize_g)) - allocate(outIceHeatFlux(lsize_g)) - ! TODO: Can we allocate these only while used or are we worried about performance hit? - ! TODO: add deallocates! - call shr_sys_flush(logunit) end if @@ -810,136 +777,6 @@ end subroutine prep_glc_map_one_state_field_lnd2glc !================================================================================================ - subroutine prep_glc_calculate_subshelf_boundary_fluxes - - !--------------------------------------------------------------- - ! Description - ! On the ice sheet grid, calculate shelf boundary fluxes - - use shr_const_mod , only: SHR_CONST_KAPPA_LAND_ICE - - ! Local Variables - - integer :: gsize, n, egi - type(mct_aVect), pointer :: o2x_ox ! Ocn export, ocn grid, cpl pes - type(mct_aVect), pointer :: x2g_gx ! Glc import, glc grid, cpl pes - type(mct_aVect), pointer :: g2x_gx ! Glc import, glc grid, cpl pes - - integer :: index_x2g_So_blt - integer :: index_x2g_So_bls - integer :: index_x2g_So_htv - integer :: index_x2g_So_stv - integer :: index_x2g_So_rhoeff - integer :: index_g2x_Sg_tbot - integer :: index_g2x_Sg_dztbot - integer :: index_g2x_Sg_lithop - integer :: index_g2x_Sg_icemask_floating - - integer :: index_g2x_Sg_blis - integer :: index_g2x_Sg_blit - integer :: index_g2x_Fogx_qiceho - integer :: index_g2x_Fogx_qicelo - integer :: index_x2g_Fogx_qiceli - integer :: index_x2g_Fogx_qicehi - - character(*), parameter :: subname = '(prep_glc_calculate_subshelf_boundary_fluxes)' - !--------------------------------------------------------------- - - if (.not.(glc_present)) return - - do egi = 1,num_inst_glc - - o2x_ox => component_get_c2x_cx(ocn(egi)) - g2x_gx => component_get_c2x_cx(glc(egi)) - x2g_gx => component_get_x2c_cx(glc(egi)) - - !Remap relevant ocean variables to ice sheet grid. - !Done here instead of in glc-frequency mapping so it happens within ocean coupling interval. - ! Also could map o2x_ox->o2x_gx(1) but using x2g_gx as destination allows us to see - ! these fields on the GLC grid of the coupler history file, which helps with debugging. - call seq_map_map(mapper_So2g, o2x_ox, x2g_gx, & - fldlist=seq_flds_x2g_shelf_states_from_ocn,norm=.true.) - - ! inputs to melt flux calculation - index_x2g_So_blt = mct_avect_indexra(x2g_gx,'So_blt',perrwith='quiet') - index_x2g_So_bls = mct_avect_indexra(x2g_gx,'So_bls',perrwith='quiet') - index_x2g_So_htv = mct_avect_indexra(x2g_gx,'So_htv',perrwith='quiet') - index_x2g_So_stv = mct_avect_indexra(x2g_gx,'So_stv',perrwith='quiet') - index_x2g_So_rhoeff = mct_avect_indexra(x2g_gx,'So_rhoeff',perrwith='quiet') - - index_g2x_Sg_tbot = mct_avect_indexra(g2x_gx,'Sg_tbot',perrwith='quiet') - index_g2x_Sg_dztbot = mct_avect_indexra(g2x_gx,'Sg_dztbot',perrwith='quiet') - index_g2x_Sg_lithop = mct_avect_indexra(g2x_gx,'Sg_lithop',perrwith='quiet') - index_g2x_Sg_icemask_floating = mct_avect_indexra(g2x_gx,'Sg_icemask_floating',perrwith='quiet') - - ! outputs to melt flux calculation - index_g2x_Sg_blis = mct_avect_indexra(g2x_gx,'Sg_blis',perrwith='quiet') - index_g2x_Sg_blit = mct_avect_indexra(g2x_gx,'Sg_blit',perrwith='quiet') - index_g2x_Fogx_qiceho = mct_avect_indexra(g2x_gx,'Fogx_qiceho',perrwith='quiet') - index_g2x_Fogx_qicelo = mct_avect_indexra(g2x_gx,'Fogx_qicelo',perrwith='quiet') - index_x2g_Fogx_qiceli = mct_avect_indexra(x2g_gx,'Fogx_qiceli',perrwith='quiet') - index_x2g_Fogx_qicehi = mct_avect_indexra(x2g_gx,'Fogx_qicehi',perrwith='quiet') - - gsize = mct_aVect_lsize(g2x_gx) - - do n=1,gsize - !Extract glc and ocn-sourced coupler fields used as input to compute_melt_fluxes to local arrays... - - ! Fields from the ocean, now on the GLC grid - oceanTemperature(n) = x2g_gx%rAttr(index_x2g_So_blt,n) - oceanSalinity(n) = x2g_gx%rAttr(index_x2g_So_bls,n) - oceanHeatTransferVelocity(n) = x2g_gx%rAttr(index_x2g_So_htv,n) - oceanSaltTransferVelocity(n) = x2g_gx%rAttr(index_x2g_So_stv,n) - - ! Fields from the ice sheet model (still on the GLC grid) - iceTemperature(n) = g2x_gx%rAttr(index_g2x_Sg_tbot,n) - iceTemperatureDistance(n) = g2x_gx%rAttr(index_g2x_Sg_dztbot,n) - interfacePressure(n) = g2x_gx%rAttr(index_g2x_Sg_lithop,n) - iceFloatingMask(n) = g2x_gx%rAttr(index_g2x_Sg_icemask_floating,n) - - !... initialize local compute_melt_fluxes output arrays... - outInterfaceSalinity(n) = 0.0_r8 - outInterfaceTemperature(n) = 0.0_r8 - outFreshwaterFlux(n) = 0.0_r8 - outOceanHeatFlux(n) = 0.0_r8 - outIceHeatFlux(n) = 0.0_r8 - end do - - !...calculate fluxes... - call compute_melt_fluxes(oceanTemperature=oceanTemperature,& - oceanSalinity=oceanSalinity,& - oceanHeatTransferVelocity=oceanHeatTransferVelocity,& - oceanSaltTransferVelocity=oceanSaltTransferVelocity,& - interfacePressure=interfacePressure,& - iceTemperature=iceTemperature,& - iceTemperatureDistance=iceTemperatureDistance, & - iceFloatingMask=iceFloatingMask, & - outInterfaceSalinity=outInterfaceSalinity,& - outInterfaceTemperature=outInterfaceTemperature,& - outFreshwaterFlux=outFreshwaterFlux,& - outOceanHeatFlux=outOceanHeatFlux,& - outIceHeatFlux=outIceHeatFlux,& - gsize=gsize) - - !...and assign fluxes to glc and ocn-directed coupler fields - do n=1,gsize - !Assign outputs from compute_melt_fluxes back into coupler attributes - g2x_gx%rAttr(index_g2x_Sg_blis,n) = outInterfaceSalinity(n) !to ocean - g2x_gx%rAttr(index_g2x_Sg_blit,n) = outInterfaceTemperature(n) !to ocean - g2x_gx%rAttr(index_g2x_Fogx_qiceho,n) = outOceanHeatFlux(n) !to ocean - g2x_gx%rAttr(index_g2x_Fogx_qicelo,n)= outFreshwaterFlux(n) !to ocean - x2g_gx%rAttr(index_x2g_Fogx_qicehi,n) = outIceHeatFlux(n) !to ice sheet - x2g_gx%rAttr(index_x2g_Fogx_qiceli,n) = -1.0_r8 * outFreshwaterFlux(n) !to ice sheet - end do - - !Note: remap ocean-side outputs back onto ocean grid done in call to prep_ocn_shelf_calc_g2x_ox - - end do ! loop over GLC instances - - end subroutine prep_glc_calculate_subshelf_boundary_fluxes - - !================================================================================================ - subroutine prep_glc_zero_fields() !--------------------------------------------------------------- @@ -1466,165 +1303,4 @@ function prep_glc_get_mapper_Fo2g() prep_glc_get_mapper_Fo2g=> mapper_Fo2g end function prep_glc_get_mapper_Fo2g -!*********************************************************************** -! -! routine compute_melt_fluxes -! -!> \brief Computes ocean and ice melt fluxes, etc. -!> \author Xylar Asay-Davis -!> \date 3/27/2015 -!> This routine computes melt fluxes (melt rate, temperature fluxes -!> into the ice and the ocean, and salt flux) as well as the interface -!> temperature and salinity. This routine expects an ice temperature -!> in the bottom layer of ice and ocean temperature and salinity in -!> the top ocean layer as well as the pressure at the ice/ocean interface. -!> -!> The ocean heat and salt transfer velocities are determined based on -!> observations of turbulent mixing rates in the under-ice boundary layer. -!> They should be the product of the friction velocity and a (possibly -!> spatially variable) non-dimenional transfer coefficient. -!> -!> The iceTemperatureDistance is the distance between the location -!> where the iceTemperature is supplied and the ice-ocean interface, -!> used to compute a temperature gradient. The ice thermal conductivity, -!> SHR_CONST_KAPPA_LAND_ICE, is zero for the freezing solution from Holland and Jenkins -!> (1999) in which the ice is purely insulating. -! -!----------------------------------------------------------------------- - - subroutine compute_melt_fluxes( & - oceanTemperature, & - oceanSalinity, & - oceanHeatTransferVelocity, & - oceanSaltTransferVelocity, & - interfacePressure, & - iceTemperature, & - iceTemperatureDistance, & - iceFloatingMask, & - outInterfaceSalinity, & - outInterfaceTemperature, & - outFreshwaterFlux, & - outOceanHeatFlux, & - outIceHeatFlux, & - gsize) - - use shr_const_mod, only: SHR_CONST_CPICE, & - SHR_CONST_CPSW, & - SHR_CONST_LATICE, & - SHR_CONST_RHOICE, & - SHR_CONST_RHOSW, & - SHR_CONST_DTF_DP, & - SHR_CONST_DTF_DS, & - SHR_CONST_DTF_DPDS, & - SHR_CONST_TF0, & - SHR_CONST_KAPPA_LAND_ICE - - !----------------------------------------------------------------- - ! - ! input variables - ! - !----------------------------------------------------------------- - - real (kind=r8), dimension(:), intent(in) :: & - oceanTemperature, & !< Input: ocean temperature in top layer - oceanSalinity, & !< Input: ocean salinity in top layer - oceanHeatTransferVelocity, & !< Input: ocean heat transfer velocity - oceanSaltTransferVelocity, & !< Input: ocean salt transfer velocity - interfacePressure, & !< Input: pressure at the ice-ocean interface - iceTemperature, & !< Input: ice temperature in bottom layer - iceTemperatureDistance !< Input: distance to ice temperature from ice-ocean interface - integer, dimension(:), intent(in) :: & - iceFloatingMask !< Input: mask of cells that contain floating ice - - integer, intent(in) :: gsize !< Input: number of values in each array - - !----------------------------------------------------------------- - ! - ! output variables - ! - !----------------------------------------------------------------- - - real (kind=r8), dimension(:), intent(out) :: & - outInterfaceSalinity, & !< Output: ocean salinity at the interface - outInterfaceTemperature, & !< Output: ice/ocean temperature at the interface - outFreshwaterFlux, & !< Output: ocean thickness flux (melt rate) - outOceanHeatFlux, & !< Output: the temperature flux into the ocean - outIceHeatFlux !< Output: the temperature flux into the ice - - !----------------------------------------------------------------- - ! - ! local variables - ! - !----------------------------------------------------------------- - - real (kind=r8) :: T0, transferVelocityRatio, Tlatent, nu, a, b, c, eta, & - iceHeatFluxCoeff, iceDeltaT, dTf_dS - integer :: n - character(*), parameter :: subname = '(compute_melt_fluxes)' - - real (kind=r8), parameter :: minInterfaceSalinity = 0.001_r8 - - real (kind=r8), parameter :: referencePressure = 0.0_r8 ! Using reference pressure of 0 - - real (kind=r8) :: pressureOffset - - Tlatent = SHR_CONST_LATICE/SHR_CONST_CPSW - do n = 1, gsize - if (iceFloatingMask(n) == 0) cycle ! Only calculate on floating cells - - if (oceanHeatTransferVelocity(n) == 0.0_r8) then - write(logunit,*) 'compute_melt_fluxes ERROR: oceanHeatTransferVelocity value of 0 causes divide by 0 at index ', n - call shr_sys_abort('compute_melt_fluxes ERROR: oceanHeatTransferVelocity value of 0 causes divide by 0') - end if - - iceHeatFluxCoeff = SHR_CONST_RHOICE*SHR_CONST_CPICE*SHR_CONST_KAPPA_LAND_ICE/iceTemperatureDistance(n) - nu = iceHeatFluxCoeff/(SHR_CONST_RHOSW*SHR_CONST_CPSW*oceanHeatTransferVelocity(n)) - pressureOffset = max(interfacePressure(n) - referencePressure, 0.0_r8) - T0 = SHR_CONST_TF0 + SHR_CONST_DTF_DP * pressureOffset - !Note: These two terms for T0 are not needed because we are evaluating at salinity=0: - !+ SHR_CONST_DTF_DS * oceanSalinity(n) + SHR_CONST_DTF_DPDS * pressureOffset * oceanSalinity(n) - iceDeltaT = T0 - iceTemperature(n) - dTf_dS = SHR_CONST_DTF_DS + SHR_CONST_DTF_DPDS * pressureOffset - - transferVelocityRatio = oceanSaltTransferVelocity(n)/oceanHeatTransferVelocity(n) - - a = -1.0_r8 * dTf_dS * (1.0_r8 + nu) - b = transferVelocityRatio*Tlatent - nu*iceDeltaT + oceanTemperature(n) - T0 - c = -transferVelocityRatio*Tlatent*max(oceanSalinity(n), 0.0_r8) - ! a is non-negative; c is strictly non-positive so we never get imaginary roots. - ! Since a can be zero, we need a solution of the quadratic equation for 1/Si instead of Si. - ! Following: https://people.csail.mit.edu/bkph/articles/Quadratics.pdf - ! Since a and -c are are non-negative, the term in the square root is also always >= |b|. - ! In all reasonable cases, b will be strictly positive, since transferVelocityRatio*Tlatent ~ 2 C, - ! T0 ~ -1.8 C and oceanTemperature should never be able to get below about -3 C - ! As long as either b or both a and c are greater than zero, the strictly non-negative root is - outInterfaceSalinity(n) = max(-(2.0_r8*c)/(b + sqrt(b**2 - 4.0_r8*a*c)), minInterfaceSalinity) - - outInterfaceTemperature(n) = dTf_dS*outInterfaceSalinity(n)+T0 - - outFreshwaterFlux(n) = SHR_CONST_RHOSW*oceanSaltTransferVelocity(n) & - * (oceanSalinity(n)/outInterfaceSalinity(n) - 1.0_r8) - - ! According to Jenkins et al. (2001), the temperature fluxes into the ocean are: - ! 1. the advection of meltwater into the top layer (or removal for freezing) - ! 2. the turbulent transfer of heat across the boundary layer, based on the termal driving - outOceanHeatFlux(n) = SHR_CONST_CPSW*(outFreshwaterFlux(n)*outInterfaceTemperature(n) & - - SHR_CONST_RHOSW*oceanHeatTransferVelocity(n)*(oceanTemperature(n)-outInterfaceTemperature(n))) - - ! the temperature fluxes into the ice are: - ! 1. the advection of ice at the interface temperature out of the domain due to melting - ! (or in due to freezing) - ! 2. the diffusion (if any) of heat into the ice, based on temperature difference between - ! the reference point in the ice (either the surface or the middle of the bottom layer) - ! and the interface - outIceHeatFlux(n) = -SHR_CONST_CPICE*outFreshwaterFlux(n)*outInterfaceTemperature(n) - - outIceHeatFlux(n) = outIceHeatFlux(n) & - - iceHeatFluxCoeff*(iceTemperature(n) - outInterfaceTemperature(n)) - - end do - - !-------------------------------------------------------------------- - end subroutine compute_melt_fluxes - end module prep_glc_mod diff --git a/driver-moab/main/prep_ocn_mod.F90 b/driver-moab/main/prep_ocn_mod.F90 index 4583805ed38a..f31c8ba07169 100644 --- a/driver-moab/main/prep_ocn_mod.F90 +++ b/driver-moab/main/prep_ocn_mod.F90 @@ -917,14 +917,14 @@ subroutine prep_ocn_init(infodata, atm_c2_ocn, atm_c2_ice, ice_c2_ocn, rof_c2_oc if (samegrid_ro) then ! this creates a parallel communication graph between mbrxid and mboxid, ! with ids rof(1)%cplcompid, ocn(1)%cplcompid - ! this will be used in send/receive mappers + ! this will be used in send/receive mappers ierr = iMOAB_ComputeCommGraph( mbrxid, mboxid, mpicom_CPLID, mpigrp_CPLID, mpigrp_CPLID, & type_grid, type_grid, rof(1)%cplcompid, ocn(1)%cplcompid ) if (ierr .ne. 0) then write(logunit,*) subname,' error in compute graph ROF - ocean ' call shr_sys_abort(subname//' ERROR in compute graph ROF - ocean ') endif - mapper_Fr2o%intx_context = ocn(1)%cplcompid + mapper_Fr2o%intx_context = ocn(1)%cplcompid mapper_Rr2o_ice%intx_context = ocn(1)%cplcompid mapper_Rr2o_liq%intx_context = ocn(1)%cplcompid @@ -2014,6 +2014,756 @@ subroutine prep_ocn_mrg_moab(infodata, xao_ox, timer_mrg) end subroutine prep_ocn_mrg_moab !================================================================================================ + subroutine prep_ocn_merge( flux_epbalfact, a2x_o, i2x_o, r2x_o, w2x_o, g2x_o, xao_o, & + fractions_o, x2o_o ) + + use seq_flds_mod, only: wav_ocn_coup + + !----------------------------------------------------------------------- + ! + ! Arguments + real(R8) , intent(in) :: flux_epbalfact + type(mct_aVect), intent(in) :: a2x_o + type(mct_aVect), intent(in) :: i2x_o + type(mct_aVect), intent(in) :: r2x_o + type(mct_aVect), intent(in) :: w2x_o + type(mct_aVect), intent(in) :: g2x_o + type(mct_aVect), intent(in) :: xao_o + type(mct_aVect), intent(in) :: fractions_o + type(mct_aVect), intent(inout) :: x2o_o + ! + ! Local variables + integer :: n,ka,ki,ko,kr,kw,kx,kir,kor,i,i1,o1 + integer :: kof,kif + integer :: lsize + integer :: noflds,naflds,niflds,nrflds,nwflds,nxflds,ngflds + real(R8) :: ifrac,ifracr + real(R8) :: afrac,afracr + real(R8) :: frac_sum + real(R8) :: avsdr, anidr, avsdf, anidf ! albedos + real(R8) :: fswabsv, fswabsi ! sw + character(CL),allocatable :: field_ocn(:) ! string converted to char + character(CL),allocatable :: field_atm(:) ! string converted to char + character(CL),allocatable :: field_ice(:) ! string converted to char + character(CL),allocatable :: field_rof(:) ! string converted to char + character(CL),allocatable :: field_wav(:) ! string converted to char + character(CL),allocatable :: field_xao(:) ! string converted to char + character(CL),allocatable :: field_glc(:) ! string converted to char + character(CL),allocatable :: itemc_ocn(:) ! string converted to char + character(CL),allocatable :: itemc_atm(:) ! string converted to char + character(CL),allocatable :: itemc_ice(:) ! string converted to char + character(CL),allocatable :: itemc_rof(:) ! string converted to char + character(CL),allocatable :: itemc_wav(:) ! string converted to char + character(CL),allocatable :: itemc_xao(:) ! string converted to char + character(CL),allocatable :: itemc_g2x(:) ! string converted to char + integer, save :: index_a2x_Faxa_swvdr + integer, save :: index_a2x_Faxa_swvdf + integer, save :: index_a2x_Faxa_swndr + integer, save :: index_a2x_Faxa_swndf + integer, save :: index_i2x_Fioi_swpen + integer, save :: index_xao_So_avsdr + integer, save :: index_xao_So_anidr + integer, save :: index_xao_So_avsdf + integer, save :: index_xao_So_anidf + integer, save :: index_a2x_Faxa_snowc + integer, save :: index_a2x_Faxa_snowl + integer, save :: index_a2x_Faxa_rainc + integer, save :: index_a2x_Faxa_rainl + integer, save :: index_r2x_Forr_rofl + integer, save :: index_r2x_Forr_rofi + integer, save :: index_r2x_Forr_rofDIN + integer, save :: index_r2x_Forr_rofDIP + integer, save :: index_r2x_Forr_rofDON + integer, save :: index_r2x_Forr_rofDOP + integer, save :: index_r2x_Forr_rofDOC + integer, save :: index_r2x_Forr_rofPP + integer, save :: index_r2x_Forr_rofDSi + integer, save :: index_r2x_Forr_rofPOC + integer, save :: index_r2x_Forr_rofPN + integer, save :: index_r2x_Forr_rofDIC + integer, save :: index_r2x_Forr_rofFe + integer, save :: index_r2x_Forr_rofl_16O + integer, save :: index_r2x_Forr_rofi_16O + integer, save :: index_r2x_Forr_rofl_18O + integer, save :: index_r2x_Forr_rofi_18O + integer, save :: index_r2x_Forr_rofl_HDO + integer, save :: index_r2x_Forr_rofi_HDO + integer, save :: index_r2x_Flrr_flood + integer, save :: index_g2x_Fogg_rofl + integer, save :: index_g2x_Fogg_rofi + integer, save :: index_x2o_Foxx_swnet + integer, save :: index_x2o_Faxa_snow + integer, save :: index_x2o_Faxa_rain + integer, save :: index_x2o_Faxa_prec + integer, save :: index_x2o_Foxx_rofl + integer, save :: index_x2o_Foxx_rofi + integer, save :: index_x2o_Foxx_rofDIN + integer, save :: index_x2o_Foxx_rofDIP + integer, save :: index_x2o_Foxx_rofDON + integer, save :: index_x2o_Foxx_rofDOP + integer, save :: index_x2o_Foxx_rofDOC + integer, save :: index_x2o_Foxx_rofPP + integer, save :: index_x2o_Foxx_rofDSi + integer, save :: index_x2o_Foxx_rofPOC + integer, save :: index_x2o_Foxx_rofPN + integer, save :: index_x2o_Foxx_rofDIC + integer, save :: index_x2o_Foxx_rofFe + integer, save :: index_x2o_Sf_afrac + integer, save :: index_x2o_Sf_afracr + integer, save :: index_x2o_Foxx_swnet_afracr + integer, save :: index_x2o_Foxx_rofl_16O + integer, save :: index_x2o_Foxx_rofi_16O + integer, save :: index_x2o_Foxx_rofl_18O + integer, save :: index_x2o_Foxx_rofi_18O + integer, save :: index_x2o_Foxx_rofl_HDO + integer, save :: index_x2o_Foxx_rofi_HDO + integer, save :: index_a2x_Faxa_snowc_16O + integer, save :: index_a2x_Faxa_snowl_16O + integer, save :: index_a2x_Faxa_rainc_16O + integer, save :: index_a2x_Faxa_rainl_16O + integer, save :: index_x2o_Faxa_rain_16O + integer, save :: index_x2o_Faxa_snow_16O + integer, save :: index_x2o_Faxa_prec_16O + integer, save :: index_a2x_Faxa_snowc_18O + integer, save :: index_a2x_Faxa_snowl_18O + integer, save :: index_a2x_Faxa_rainc_18O + integer, save :: index_a2x_Faxa_rainl_18O + integer, save :: index_x2o_Faxa_rain_18O + integer, save :: index_x2o_Faxa_snow_18O + integer, save :: index_x2o_Faxa_prec_18O + integer, save :: index_a2x_Faxa_snowc_HDO + integer, save :: index_a2x_Faxa_snowl_HDO + integer, save :: index_a2x_Faxa_rainc_HDO + integer, save :: index_a2x_Faxa_rainl_HDO + integer, save :: index_x2o_Faxa_rain_HDO + integer, save :: index_x2o_Faxa_snow_HDO + integer, save :: index_x2o_Faxa_prec_HDO + logical :: iamroot + logical, save, pointer :: amerge(:),imerge(:),xmerge(:) + integer, save, pointer :: aindx(:), iindx(:), xindx(:) + character(CL),allocatable :: mrgstr(:) ! temporary string + type(mct_aVect_sharedindices),save :: a2x_sharedindices + type(mct_aVect_sharedindices),save :: i2x_sharedindices + type(mct_aVect_sharedindices),save :: r2x_sharedindices + type(mct_aVect_sharedindices),save :: w2x_sharedindices + type(mct_aVect_sharedindices),save :: xao_sharedindices + type(mct_aVect_sharedindices),save :: g2x_sharedindices + logical, save :: first_time = .true. + character(*),parameter :: subName = '(prep_ocn_merge) ' + + !----------------------------------------------------------------------- + + call seq_comm_setptrs(CPLID, iamroot=iamroot) + + noflds = mct_aVect_nRattr(x2o_o) + naflds = mct_aVect_nRattr(a2x_o) + niflds = mct_aVect_nRattr(i2x_o) + nrflds = mct_aVect_nRattr(r2x_o) + nwflds = mct_aVect_nRattr(w2x_o) + nxflds = mct_aVect_nRattr(xao_o) + ngflds = mct_aVect_nRattr(g2x_o) + + if (first_time) then + index_a2x_Faxa_swvdr = mct_aVect_indexRA(a2x_o,'Faxa_swvdr') + index_a2x_Faxa_swvdf = mct_aVect_indexRA(a2x_o,'Faxa_swvdf') + index_a2x_Faxa_swndr = mct_aVect_indexRA(a2x_o,'Faxa_swndr') + index_a2x_Faxa_swndf = mct_aVect_indexRA(a2x_o,'Faxa_swndf') + index_i2x_Fioi_swpen = mct_aVect_indexRA(i2x_o,'Fioi_swpen') + index_xao_So_avsdr = mct_aVect_indexRA(xao_o,'So_avsdr') + index_xao_So_anidr = mct_aVect_indexRA(xao_o,'So_anidr') + index_xao_So_avsdf = mct_aVect_indexRA(xao_o,'So_avsdf') + index_xao_So_anidf = mct_aVect_indexRA(xao_o,'So_anidf') + index_x2o_Foxx_swnet = mct_aVect_indexRA(x2o_o,'Foxx_swnet') + + index_a2x_Faxa_snowc = mct_aVect_indexRA(a2x_o,'Faxa_snowc') + index_a2x_Faxa_snowl = mct_aVect_indexRA(a2x_o,'Faxa_snowl') + index_a2x_Faxa_rainc = mct_aVect_indexRA(a2x_o,'Faxa_rainc') + index_a2x_Faxa_rainl = mct_aVect_indexRA(a2x_o,'Faxa_rainl') + index_r2x_Forr_rofl = mct_aVect_indexRA(r2x_o,'Forr_rofl') + index_r2x_Forr_rofi = mct_aVect_indexRA(r2x_o,'Forr_rofi') + if (rof2ocn_nutrients) then + index_r2x_Forr_rofDIN = mct_aVect_indexRA(r2x_o,'Forr_rofDIN') + index_r2x_Forr_rofDIP = mct_aVect_indexRA(r2x_o,'Forr_rofDIP') + index_r2x_Forr_rofDON = mct_aVect_indexRA(r2x_o,'Forr_rofDON') + index_r2x_Forr_rofDOP = mct_aVect_indexRA(r2x_o,'Forr_rofDOP') + index_r2x_Forr_rofDOC = mct_aVect_indexRA(r2x_o,'Forr_rofDOC') + index_r2x_Forr_rofPP = mct_aVect_indexRA(r2x_o,'Forr_rofPP') + index_r2x_Forr_rofDSi = mct_aVect_indexRA(r2x_o,'Forr_rofDSi') + index_r2x_Forr_rofPOC = mct_aVect_indexRA(r2x_o,'Forr_rofPOC') + index_r2x_Forr_rofPN = mct_aVect_indexRA(r2x_o,'Forr_rofPN') + index_r2x_Forr_rofDIC = mct_aVect_indexRA(r2x_o,'Forr_rofDIC') + index_r2x_Forr_rofFe = mct_aVect_indexRA(r2x_o,'Forr_rofFe') + endif + index_r2x_Flrr_flood = mct_aVect_indexRA(r2x_o,'Flrr_flood') + index_g2x_Fogg_rofl = mct_aVect_indexRA(g2x_o,'Fogg_rofl') + index_g2x_Fogg_rofi = mct_aVect_indexRA(g2x_o,'Fogg_rofi') + index_x2o_Faxa_snow = mct_aVect_indexRA(x2o_o,'Faxa_snow') + index_x2o_Faxa_rain = mct_aVect_indexRA(x2o_o,'Faxa_rain') + index_x2o_Faxa_prec = mct_aVect_indexRA(x2o_o,'Faxa_prec') + index_x2o_Foxx_rofl = mct_aVect_indexRA(x2o_o,'Foxx_rofl') + index_x2o_Foxx_rofi = mct_aVect_indexRA(x2o_o,'Foxx_rofi') + if (rof2ocn_nutrients) then + index_x2o_Foxx_rofDIN = mct_aVect_indexRA(x2o_o,'Foxx_rofDIN') + index_x2o_Foxx_rofDIP = mct_aVect_indexRA(x2o_o,'Foxx_rofDIP') + index_x2o_Foxx_rofDON = mct_aVect_indexRA(x2o_o,'Foxx_rofDON') + index_x2o_Foxx_rofDOP = mct_aVect_indexRA(x2o_o,'Foxx_rofDOP') + index_x2o_Foxx_rofDOC = mct_aVect_indexRA(x2o_o,'Foxx_rofDOC') + index_x2o_Foxx_rofPP = mct_aVect_indexRA(x2o_o,'Foxx_rofPP') + index_x2o_Foxx_rofDSi = mct_aVect_indexRA(x2o_o,'Foxx_rofDSi') + index_x2o_Foxx_rofPOC = mct_aVect_indexRA(x2o_o,'Foxx_rofPOC') + index_x2o_Foxx_rofPN = mct_aVect_indexRA(x2o_o,'Foxx_rofPN') + index_x2o_Foxx_rofDIC = mct_aVect_indexRA(x2o_o,'Foxx_rofDIC') + index_x2o_Foxx_rofFe = mct_aVect_indexRA(x2o_o,'Foxx_rofFe') + endif + + if (seq_flds_i2o_per_cat) then + index_x2o_Sf_afrac = mct_aVect_indexRA(x2o_o,'Sf_afrac') + index_x2o_Sf_afracr = mct_aVect_indexRA(x2o_o,'Sf_afracr') + index_x2o_Foxx_swnet_afracr = mct_aVect_indexRA(x2o_o,'Foxx_swnet_afracr') + endif + + !wiso: + ! H2_16O + index_a2x_Faxa_snowc_16O = mct_aVect_indexRA(a2x_o,'Faxa_snowc_16O', perrWith='quiet') + index_a2x_Faxa_snowl_16O = mct_aVect_indexRA(a2x_o,'Faxa_snowl_16O', perrWith='quiet') + index_a2x_Faxa_rainc_16O = mct_aVect_indexRA(a2x_o,'Faxa_rainc_16O', perrWith='quiet') + index_a2x_Faxa_rainl_16O = mct_aVect_indexRA(a2x_o,'Faxa_rainl_16O', perrWith='quiet') + index_r2x_Forr_rofl_16O = mct_aVect_indexRA(r2x_o,'Forr_rofl_16O' , perrWith='quiet') + index_r2x_Forr_rofi_16O = mct_aVect_indexRA(r2x_o,'Forr_rofi_16O' , perrWith='quiet') + index_x2o_Faxa_rain_16O = mct_aVect_indexRA(x2o_o,'Faxa_rain_16O' , perrWith='quiet') + index_x2o_Faxa_snow_16O = mct_aVect_indexRA(x2o_o,'Faxa_snow_16O' , perrWith='quiet') + index_x2o_Faxa_prec_16O = mct_aVect_indexRA(x2o_o,'Faxa_prec_16O' , perrWith='quiet') + index_x2o_Foxx_rofl_16O = mct_aVect_indexRA(x2o_o,'Foxx_rofl_16O' , perrWith='quiet') + index_x2o_Foxx_rofi_16O = mct_aVect_indexRA(x2o_o,'Foxx_rofi_16O' , perrWith='quiet') + ! H2_18O + index_a2x_Faxa_snowc_18O = mct_aVect_indexRA(a2x_o,'Faxa_snowc_18O', perrWith='quiet') + index_a2x_Faxa_snowl_18O = mct_aVect_indexRA(a2x_o,'Faxa_snowl_18O', perrWith='quiet') + index_a2x_Faxa_rainc_18O = mct_aVect_indexRA(a2x_o,'Faxa_rainc_18O', perrWith='quiet') + index_a2x_Faxa_rainl_18O = mct_aVect_indexRA(a2x_o,'Faxa_rainl_18O', perrWith='quiet') + index_r2x_Forr_rofl_18O = mct_aVect_indexRA(r2x_o,'Forr_rofl_18O' , perrWith='quiet') + index_r2x_Forr_rofi_18O = mct_aVect_indexRA(r2x_o,'Forr_rofi_18O' , perrWith='quiet') + index_x2o_Faxa_rain_18O = mct_aVect_indexRA(x2o_o,'Faxa_rain_18O' , perrWith='quiet') + index_x2o_Faxa_snow_18O = mct_aVect_indexRA(x2o_o,'Faxa_snow_18O' , perrWith='quiet') + index_x2o_Faxa_prec_18O = mct_aVect_indexRA(x2o_o,'Faxa_prec_18O' , perrWith='quiet') + index_x2o_Foxx_rofl_18O = mct_aVect_indexRA(x2o_o,'Foxx_rofl_18O' , perrWith='quiet') + index_x2o_Foxx_rofi_18O = mct_aVect_indexRA(x2o_o,'Foxx_rofi_18O' , perrWith='quiet') + ! HDO + index_a2x_Faxa_snowc_HDO = mct_aVect_indexRA(a2x_o,'Faxa_snowc_HDO', perrWith='quiet') + index_a2x_Faxa_snowl_HDO = mct_aVect_indexRA(a2x_o,'Faxa_snowl_HDO', perrWith='quiet') + index_a2x_Faxa_rainc_HDO = mct_aVect_indexRA(a2x_o,'Faxa_rainc_HDO', perrWith='quiet') + index_a2x_Faxa_rainl_HDO = mct_aVect_indexRA(a2x_o,'Faxa_rainl_HDO', perrWith='quiet') + index_r2x_Forr_rofl_HDO = mct_aVect_indexRA(r2x_o,'Forr_rofl_HDO' , perrWith='quiet') + index_r2x_Forr_rofi_HDO = mct_aVect_indexRA(r2x_o,'Forr_rofi_HDO' , perrWith='quiet') + index_x2o_Faxa_rain_HDO = mct_aVect_indexRA(x2o_o,'Faxa_rain_HDO' , perrWith='quiet') + index_x2o_Faxa_snow_HDO = mct_aVect_indexRA(x2o_o,'Faxa_snow_HDO' , perrWith='quiet') + index_x2o_Faxa_prec_HDO = mct_aVect_indexRA(x2o_o,'Faxa_prec_HDO' , perrWith='quiet') + index_x2o_Foxx_rofl_HDO = mct_aVect_indexRA(x2o_o,'Foxx_rofl_HDO' , perrWith='quiet') + index_x2o_Foxx_rofi_HDO = mct_aVect_indexRA(x2o_o,'Foxx_rofi_HDO' , perrWith='quiet') + + ! Compute all other quantities based on standardized naming convention (see below) + ! Only ocn field states that have the name-prefix Sx_ will be merged + ! Only field names have the same name-suffix (after the "_") will be merged + ! (e.g. Si_fldname, Sa_fldname => merged to => Sx_fldname) + ! All fluxes will be scaled by the corresponding afrac or ifrac + ! EXCEPT for + ! -- Faxa_snnet, Faxa_snow, Faxa_rain, Faxa_prec (derived) + ! All i2x_o fluxes that have the name-suffix "Faii" (atm/ice fluxes) will be ignored + ! - only ice fluxes that are Fioi_... will be used in the ocean merges + + allocate(aindx(noflds), amerge(noflds)) + allocate(iindx(noflds), imerge(noflds)) + allocate(xindx(noflds), xmerge(noflds)) + allocate(field_atm(naflds), itemc_atm(naflds)) + allocate(field_ice(niflds), itemc_ice(niflds)) + allocate(field_ocn(noflds), itemc_ocn(noflds)) + allocate(field_rof(nrflds), itemc_rof(nrflds)) + allocate(field_wav(nwflds), itemc_wav(nwflds)) + allocate(field_xao(nxflds), itemc_xao(nxflds)) + allocate(field_glc(ngflds), itemc_g2x(ngflds)) + allocate(mrgstr(noflds)) + aindx(:) = 0 + iindx(:) = 0 + xindx(:) = 0 + amerge(:) = .true. + imerge(:) = .true. + xmerge(:) = .true. + + do ko = 1,noflds + field_ocn(ko) = mct_aVect_getRList2c(ko, x2o_o) + itemc_ocn(ko) = trim(field_ocn(ko)(scan(field_ocn(ko),'_'):)) + enddo + do ka = 1,naflds + field_atm(ka) = mct_aVect_getRList2c(ka, a2x_o) + itemc_atm(ka) = trim(field_atm(ka)(scan(field_atm(ka),'_'):)) + enddo + do ki = 1,niflds + field_ice(ki) = mct_aVect_getRList2c(ki, i2x_o) + itemc_ice(ki) = trim(field_ice(ki)(scan(field_ice(ki),'_'):)) + enddo + do kr = 1,nrflds + field_rof(kr) = mct_aVect_getRList2c(kr, r2x_o) + itemc_rof(kr) = trim(field_rof(kr)(scan(field_rof(kr),'_'):)) + enddo + do kw = 1,nwflds + field_wav(kw) = mct_aVect_getRList2c(kw, w2x_o) + itemc_wav(kw) = trim(field_wav(kw)(scan(field_wav(kw),'_'):)) + enddo + do kx = 1,nxflds + field_xao(kx) = mct_aVect_getRList2c(kx, xao_o) + itemc_xao(kx) = trim(field_xao(kx)(scan(field_xao(kx),'_'):)) + enddo + do kx = 1,ngflds + field_glc(kx) = mct_aVect_getRList2c(kx, g2x_o) + itemc_g2x(kx) = trim(field_glc(kx)(scan(field_glc(kx),'_'):)) + enddo + + call mct_aVect_setSharedIndices(a2x_o, x2o_o, a2x_SharedIndices) + call mct_aVect_setSharedIndices(i2x_o, x2o_o, i2x_SharedIndices) + call mct_aVect_setSharedIndices(r2x_o, x2o_o, r2x_SharedIndices) + call mct_aVect_setSharedIndices(w2x_o, x2o_o, w2x_SharedIndices) + call mct_aVect_setSharedIndices(xao_o, x2o_o, xao_SharedIndices) + call mct_aVect_setSharedIndices(g2x_o, x2o_o, g2x_SharedIndices) + + do ko = 1,noflds + !--- document merge --- + mrgstr(ko) = subname//'x2o%'//trim(field_ocn(ko))//' =' + if (field_ocn(ko)(1:2) == 'PF') then + cycle ! if flux has first character as P, pass straight through + end if + if (field_ocn(ko)(1:1) == 'S' .and. field_ocn(ko)(2:2) /= 'x') then + cycle ! ignore all ocn states that do not have a Sx_ prefix + end if + if (trim(field_ocn(ko)) == 'Foxx_swnet' .or. & + trim(field_ocn(ko)) == 'Faxa_snow' .or. & + trim(field_ocn(ko)) == 'Faxa_rain' .or. & + trim(field_ocn(ko)) == 'Faxa_prec' )then + cycle ! ignore swnet, snow, rain, prec - treated explicitly above + end if + if (index(field_ocn(ko), 'Faxa_snow_' ) == 1 .or. & + index(field_ocn(ko), 'Faxa_rain_' ) == 1 .or. & + index(field_ocn(ko), 'Faxa_prec_' ) == 1 )then + cycle ! ignore isotope snow, rain, prec - treated explicitly above + end if + ! if (trim(field_ocn(ko)(1:5)) == 'Foxx_') then + ! cycle ! ignore runoff fields from land - treated in coupler + ! end if + + do ka = 1,naflds + if (trim(itemc_ocn(ko)) == trim(itemc_atm(ka))) then + if ((trim(field_ocn(ko)) == trim(field_atm(ka)))) then + if (field_atm(ka)(1:1) == 'F') amerge(ko) = .false. + end if + ! --- make sure only one field matches --- + if (aindx(ko) /= 0) then + write(logunit,*) subname,' ERROR: found multiple ka field matches for ',trim(itemc_atm(ka)) + call shr_sys_abort(subname//' ERROR multiple ka field matches') + endif + aindx(ko) = ka + end if + end do + do ki = 1,niflds + if (field_ice(ki)(1:1) == 'F' .and. field_ice(ki)(2:4) == 'aii') then + cycle ! ignore all i2x_o fluxes that are ice/atm fluxes + end if + if (trim(itemc_ocn(ko)) == trim(itemc_ice(ki))) then + if ((trim(field_ocn(ko)) == trim(field_ice(ki)))) then + if (field_ice(ki)(1:1) == 'F') imerge(ko) = .false. + end if + ! --- make sure only one field matches --- + if (iindx(ko) /= 0) then + write(logunit,*) subname,' ERROR: found multiple ki field matches for ',trim(itemc_ice(ki)) + call shr_sys_abort(subname//' ERROR multiple ki field matches') + endif + iindx(ko) = ki + end if + end do + do kx = 1,nxflds + if (trim(itemc_ocn(ko)) == trim(itemc_xao(kx))) then + if ((trim(field_ocn(ko)) == trim(field_xao(kx)))) then + if (field_xao(kx)(1:1) == 'F') xmerge(ko) = .false. + end if + ! --- make sure only one field matches --- + if (xindx(ko) /= 0) then + write(logunit,*) subname,' ERROR: found multiple kx field matches for ',trim(itemc_xao(kx)) + call shr_sys_abort(subname//' ERROR multiple kx field matches') + endif + xindx(ko) = kx + end if + end do + + ! --- add some checks --- + + ! --- make sure no merge of BOTH atm and xao --- + if (aindx(ko) > 0 .and. xindx(ko) > 0) then + write(logunit,*) subname,' ERROR: aindx and xindx both non-zero, not allowed' + call shr_sys_abort(subname//' ERROR aindx and xindx both non-zero') + endif + + ! --- make sure all terms agree on merge or non-merge aspect --- + if (aindx(ko) > 0 .and. iindx(ko) > 0 .and. (amerge(ko) .neqv. imerge(ko))) then + write(logunit,*) subname,' ERROR: aindx and iindx merge logic error' + call shr_sys_abort(subname//' ERROR aindx and iindx merge logic error') + endif + if (aindx(ko) > 0 .and. xindx(ko) > 0 .and. (amerge(ko) .neqv. xmerge(ko))) then + write(logunit,*) subname,' ERROR: aindx and xindx merge logic error' + call shr_sys_abort(subname//' ERROR aindx and xindx merge logic error') + endif + if (xindx(ko) > 0 .and. iindx(ko) > 0 .and. (xmerge(ko) .neqv. imerge(ko))) then + write(logunit,*) subname,' ERROR: xindx and iindx merge logic error' + call shr_sys_abort(subname//' ERROR xindx and iindx merge logic error') + endif + + end do + + end if + + call mct_aVect_zero(x2o_o) + + !--- document copy operations --- + if (first_time) then + !--- document merge --- + do i=1,a2x_SharedIndices%shared_real%num_indices + i1=a2x_SharedIndices%shared_real%aVindices1(i) + o1=a2x_SharedIndices%shared_real%aVindices2(i) + mrgstr(o1) = trim(mrgstr(o1))//' = a2x%'//trim(field_atm(i1)) + enddo + do i=1,i2x_SharedIndices%shared_real%num_indices + i1=i2x_SharedIndices%shared_real%aVindices1(i) + o1=i2x_SharedIndices%shared_real%aVindices2(i) + mrgstr(o1) = trim(mrgstr(o1))//' = i2x%'//trim(field_ice(i1)) + enddo + do i=1,r2x_SharedIndices%shared_real%num_indices + i1=r2x_SharedIndices%shared_real%aVindices1(i) + o1=r2x_SharedIndices%shared_real%aVindices2(i) + mrgstr(o1) = trim(mrgstr(o1))//' = r2x%'//trim(field_rof(i1)) + enddo + do i=1,w2x_SharedIndices%shared_real%num_indices + i1=w2x_SharedIndices%shared_real%aVindices1(i) + o1=w2x_SharedIndices%shared_real%aVindices2(i) + mrgstr(o1) = trim(mrgstr(o1))//' = w2x%'//trim(field_wav(i1)) + enddo + do i=1,xao_SharedIndices%shared_real%num_indices + i1=xao_SharedIndices%shared_real%aVindices1(i) + o1=xao_SharedIndices%shared_real%aVindices2(i) + mrgstr(o1) = trim(mrgstr(o1))//' = xao%'//trim(field_xao(i1)) + enddo + do i=1,g2x_SharedIndices%shared_real%num_indices + i1=g2x_SharedIndices%shared_real%aVindices1(i) + o1=g2x_SharedIndices%shared_real%aVindices2(i) + mrgstr(o1) = trim(mrgstr(o1))//' = g2x%'//trim(field_glc(i1)) + enddo + endif + + ! call mct_aVect_copy(aVin=a2x_o, aVout=x2o_o, vector=mct_usevector) + ! call mct_aVect_copy(aVin=i2x_o, aVout=x2o_o, vector=mct_usevector) + ! call mct_aVect_copy(aVin=r2x_o, aVout=x2o_o, vector=mct_usevector) + ! call mct_aVect_copy(aVin=w2x_o, aVout=x2o_o, vector=mct_usevector) + ! call mct_aVect_copy(aVin=xao_o, aVout=x2o_o, vector=mct_usevector) + call mct_aVect_copy(aVin=a2x_o, aVout=x2o_o, vector=mct_usevector, sharedIndices=a2x_SharedIndices) + call mct_aVect_copy(aVin=i2x_o, aVout=x2o_o, vector=mct_usevector, sharedIndices=i2x_SharedIndices) + call mct_aVect_copy(aVin=r2x_o, aVout=x2o_o, vector=mct_usevector, sharedIndices=r2x_SharedIndices) + if(wav_ocn_coup == 'twoway') call mct_aVect_copy(aVin=w2x_o, aVout=x2o_o, vector=mct_usevector, sharedIndices=w2x_SharedIndices) + call mct_aVect_copy(aVin=xao_o, aVout=x2o_o, vector=mct_usevector, sharedIndices=xao_SharedIndices) + call mct_aVect_copy(aVin=g2x_o, aVout=x2o_o, vector=mct_usevector, sharedIndices=g2x_SharedIndices) + + !--- document manual merges --- + if (first_time) then + mrgstr(index_x2o_Foxx_swnet) = trim(mrgstr(index_x2o_Foxx_swnet))//' = '// & + 'afracr*(a2x%Faxa_swvdr*(1.0-xao%So_avsdr) + '// & + 'a2x%Faxa_swvdf*(1.0-xao%So_avsdf) + '// & + 'a2x%Faxa_swndr*(1.0-xao%So_anidr) + '// & + 'a2x%Faxa_swndf*(1.0-xao%So_anidf)) + '// & + 'ifrac*i2x%Fioi_swpen' + if (seq_flds_i2o_per_cat) then + mrgstr(index_x2o_Foxx_swnet_afracr) = trim(mrgstr(index_x2o_Foxx_swnet_afracr))//' = '// & + 'afracr*(a2x%Faxa_swvdr*(1.0-xao%So_avsdr) + '// & + 'a2x%Faxa_swvdf*(1.0-xao%So_avsdf) + '// & + 'a2x%Faxa_swndr*(1.0-xao%So_anidr) + '// & + 'a2x%Faxa_swndf*(1.0-xao%So_anidf))' + end if + mrgstr(index_x2o_Faxa_snow) = trim(mrgstr(index_x2o_Faxa_snow))//' = '// & + 'afrac*(a2x%Faxa_snowc + a2x%Faxa_snowl)*flux_epbalfact' + mrgstr(index_x2o_Faxa_rain) = trim(mrgstr(index_x2o_Faxa_rain))//' = '// & + 'afrac*(a2x%Faxa_rainc + a2x%Faxa_rainl)*flux_epbalfact' + mrgstr(index_x2o_Faxa_prec) = trim(mrgstr(index_x2o_Faxa_prec))//' = '// & + 'afrac*(a2x%Faxa_snowc + a2x%Faxa_snowl + a2x%Faxa_rainc + a2x%Faxa_rainl)*flux_epbalfact' + mrgstr(index_x2o_Foxx_rofl) = trim(mrgstr(index_x2o_Foxx_rofl))//' = '// & + '(r2x%Forr_rofl + r2x%Flrr_flood + g2x%Fogg_rofl)*flux_epbalfact' + mrgstr(index_x2o_Foxx_rofi) = trim(mrgstr(index_x2o_Foxx_rofi))//' = '// & + '(r2x%Forr_rofi + g2x%Fogg_rofi)*flux_epbalfact' + + ! river nutrients + if (rof2ocn_nutrients) then + mrgstr(index_x2o_Foxx_rofDIN) = trim(mrgstr(index_x2o_Foxx_rofDIN))//' = '// & + 'r2x%Forr_rofDIN' + mrgstr(index_x2o_Foxx_rofDIP) = trim(mrgstr(index_x2o_Foxx_rofDIP))//' = '// & + 'r2x%Forr_rofDIP' + mrgstr(index_x2o_Foxx_rofDON) = trim(mrgstr(index_x2o_Foxx_rofDON))//' = '// & + 'r2x%Forr_rofDON' + mrgstr(index_x2o_Foxx_rofDOP) = trim(mrgstr(index_x2o_Foxx_rofDOP))//' = '// & + 'r2x%Forr_rofDOP' + mrgstr(index_x2o_Foxx_rofDOC) = trim(mrgstr(index_x2o_Foxx_rofDOC))//' = '// & + 'r2x%Forr_rofDOC' + mrgstr(index_x2o_Foxx_rofPP) = trim(mrgstr(index_x2o_Foxx_rofPP))//' = '// & + 'r2x%Forr_rofPP' + mrgstr(index_x2o_Foxx_rofDSi) = trim(mrgstr(index_x2o_Foxx_rofDSi))//' = '// & + 'r2x%Forr_rofDSi' + mrgstr(index_x2o_Foxx_rofPOC) = trim(mrgstr(index_x2o_Foxx_rofPOC))//' = '// & + 'r2x%Forr_rofPOC' + mrgstr(index_x2o_Foxx_rofPN) = trim(mrgstr(index_x2o_Foxx_rofPN))//' = '// & + 'r2x%Forr_rofPN' + mrgstr(index_x2o_Foxx_rofDIC) = trim(mrgstr(index_x2o_Foxx_rofDIC))//' = '// & + 'r2x%Forr_rofDIC' + mrgstr(index_x2o_Foxx_rofFe) = trim(mrgstr(index_x2o_Foxx_rofFe))//' = '// & + 'r2x%Forr_rofFe' + endif + + ! water isotope snow, rain prec + if ( index_x2o_Faxa_snow_16O /= 0 )then + mrgstr(index_x2o_Faxa_snow_16O) = trim(mrgstr(index_x2o_Faxa_snow_16O))//' = '// & + 'afrac*(a2x%Faxa_snowc_16O + a2x%Faxa_snowl_16O)*flux_epbalfact' + mrgstr(index_x2o_Faxa_rain_16O) = trim(mrgstr(index_x2o_Faxa_rain_16O))//' = '// & + 'afrac*(a2x%Faxa_rainc_16O + a2x%Faxa_rainl_16O)*flux_epbalfact' + mrgstr(index_x2o_Faxa_prec_16O) = trim(mrgstr(index_x2o_Faxa_prec_16O))//' = '// & + 'afrac*(a2x%Faxa_snowc_16O + a2x%Faxa_snowl_16O + a2x%Faxa_rainc_16O + '// & + 'a2x%Faxa_rainl_16O)*flux_epbalfact' + end if + if ( index_x2o_Faxa_snow_18O /= 0 )then + mrgstr(index_x2o_Faxa_snow_18O) = trim(mrgstr(index_x2o_Faxa_snow_18O))//' = '// & + 'afrac*(a2x%Faxa_snowc_18O + a2x%Faxa_snowl_18O)*flux_epbalfact' + mrgstr(index_x2o_Faxa_rain_18O) = trim(mrgstr(index_x2o_Faxa_rain_18O))//' = '// & + 'afrac*(a2x%Faxa_rainc_18O + a2x%Faxa_rainl_18O)*flux_epbalfact' + mrgstr(index_x2o_Faxa_prec_18O) = trim(mrgstr(index_x2o_Faxa_prec_18O))//' = '// & + 'afrac*(a2x%Faxa_snowc_18O + a2x%Faxa_snowl_18O + a2x%Faxa_rainc_18O + '// & + 'a2x%Faxa_rainl_18O)*flux_epbalfact' + end if + if ( index_x2o_Faxa_snow_HDO /= 0 )then + mrgstr(index_x2o_Faxa_snow_HDO) = trim(mrgstr(index_x2o_Faxa_snow_HDO))//' = '// & + 'afrac*(a2x%Faxa_snowc_HDO + a2x%Faxa_snowl_HDO)*flux_epbalfact' + mrgstr(index_x2o_Faxa_rain_HDO) = trim(mrgstr(index_x2o_Faxa_rain_HDO))//' = '// & + 'afrac*(a2x%Faxa_rainc_HDO + a2x%Faxa_rainl_HDO)*flux_epbalfact' + mrgstr(index_x2o_Faxa_prec_HDO) = trim(mrgstr(index_x2o_Faxa_prec_HDO))//' = '// & + 'afrac*(a2x%Faxa_snowc_HDO + a2x%Faxa_snowl_HDO + a2x%Faxa_rainc_HDO + '// & + 'a2x%Faxa_rainl_HDO)*flux_epbalfact' + end if + endif + + ! Compute input ocn state (note that this only applies to non-land portion of gridcell) + + kif = mct_aVect_indexRa(fractions_o,"ifrac",perrWith=subName) + kof = mct_aVect_indexRa(fractions_o,"ofrac",perrWith=subName) + kir = mct_aVect_indexRa(fractions_o,"ifrad",perrWith=subName) + kor = mct_aVect_indexRa(fractions_o,"ofrad",perrWith=subName) + lsize = mct_aVect_lsize(x2o_o) + do n = 1,lsize + + ifrac = fractions_o%rAttr(kif,n) + afrac = fractions_o%rAttr(kof,n) + frac_sum = ifrac + afrac + if ((frac_sum) /= 0._R8) then + ifrac = ifrac / (frac_sum) + afrac = afrac / (frac_sum) + endif + + ifracr = fractions_o%rAttr(kir,n) + afracr = fractions_o%rAttr(kor,n) + frac_sum = ifracr + afracr + if ((frac_sum) /= 0._R8) then + ifracr = ifracr / (frac_sum) + afracr = afracr / (frac_sum) + endif + + ! Derived: compute net short-wave + avsdr = xao_o%rAttr(index_xao_So_avsdr,n) + anidr = xao_o%rAttr(index_xao_So_anidr,n) + avsdf = xao_o%rAttr(index_xao_So_avsdf,n) + anidf = xao_o%rAttr(index_xao_So_anidf,n) + fswabsv = a2x_o%rAttr(index_a2x_Faxa_swvdr,n) * (1.0_R8 - avsdr) & + + a2x_o%rAttr(index_a2x_Faxa_swvdf,n) * (1.0_R8 - avsdf) + fswabsi = a2x_o%rAttr(index_a2x_Faxa_swndr,n) * (1.0_R8 - anidr) & + + a2x_o%rAttr(index_a2x_Faxa_swndf,n) * (1.0_R8 - anidf) + x2o_o%rAttr(index_x2o_Foxx_swnet,n) = (fswabsv + fswabsi) * afracr + & + i2x_o%rAttr(index_i2x_Fioi_swpen,n) * ifrac + + if (seq_flds_i2o_per_cat) then + x2o_o%rAttr(index_x2o_Sf_afrac,n) = afrac + x2o_o%rAttr(index_x2o_Sf_afracr,n) = afracr + x2o_o%rAttr(index_x2o_Foxx_swnet_afracr,n) = (fswabsv + fswabsi) * afracr + end if + + ! Derived: compute total precipitation - scale total precip and runoff + + x2o_o%rAttr(index_x2o_Faxa_snow ,n) = a2x_o%rAttr(index_a2x_Faxa_snowc,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_snowl,n) * afrac + x2o_o%rAttr(index_x2o_Faxa_rain ,n) = a2x_o%rAttr(index_a2x_Faxa_rainc,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_rainl,n) * afrac + + x2o_o%rAttr(index_x2o_Faxa_snow ,n) = x2o_o%rAttr(index_x2o_Faxa_snow ,n) * flux_epbalfact + x2o_o%rAttr(index_x2o_Faxa_rain ,n) = x2o_o%rAttr(index_x2o_Faxa_rain ,n) * flux_epbalfact + + x2o_o%rAttr(index_x2o_Faxa_prec ,n) = x2o_o%rAttr(index_x2o_Faxa_rain ,n) + & + x2o_o%rAttr(index_x2o_Faxa_snow ,n) + + x2o_o%rAttr(index_x2o_Foxx_rofl, n) = (r2x_o%rAttr(index_r2x_Forr_rofl , n) + & + r2x_o%rAttr(index_r2x_Flrr_flood, n) + & + g2x_o%rAttr(index_g2x_Fogg_rofl , n)) * flux_epbalfact + x2o_o%rAttr(index_x2o_Foxx_rofi, n) = (r2x_o%rAttr(index_r2x_Forr_rofi , n) + & + g2x_o%rAttr(index_g2x_Fogg_rofi , n)) * flux_epbalfact + + if (rof2ocn_nutrients) then + x2o_o%rAttr(index_x2o_Foxx_rofDIN, n) = r2x_o%rAttr(index_r2x_Forr_rofDIN , n) + x2o_o%rAttr(index_x2o_Foxx_rofDIP, n) = r2x_o%rAttr(index_r2x_Forr_rofDIP , n) + x2o_o%rAttr(index_x2o_Foxx_rofDON, n) = r2x_o%rAttr(index_r2x_Forr_rofDON , n) + x2o_o%rAttr(index_x2o_Foxx_rofDOP, n) = r2x_o%rAttr(index_r2x_Forr_rofDOP , n) + x2o_o%rAttr(index_x2o_Foxx_rofDOC, n) = r2x_o%rAttr(index_r2x_Forr_rofDOC , n) + x2o_o%rAttr(index_x2o_Foxx_rofPP , n) = r2x_o%rAttr(index_r2x_Forr_rofPP , n) + x2o_o%rAttr(index_x2o_Foxx_rofDSi, n) = r2x_o%rAttr(index_r2x_Forr_rofDSi , n) + x2o_o%rAttr(index_x2o_Foxx_rofPOC, n) = r2x_o%rAttr(index_r2x_Forr_rofPOC , n) + x2o_o%rAttr(index_x2o_Foxx_rofPN , n) = r2x_o%rAttr(index_r2x_Forr_rofPN , n) + x2o_o%rAttr(index_x2o_Foxx_rofDIC, n) = r2x_o%rAttr(index_r2x_Forr_rofDIC , n) + x2o_o%rAttr(index_x2o_Foxx_rofFe, n) = r2x_o%rAttr(index_r2x_Forr_rofFe , n) + endif + + if ( index_x2o_Foxx_rofl_16O /= 0 ) then + x2o_o%rAttr(index_x2o_Foxx_rofl_16O, n) = (r2x_o%rAttr(index_r2x_Forr_rofl_16O, n) + & + r2x_o%rAttr(index_r2x_Flrr_flood, n) + & + g2x_o%rAttr(index_g2x_Fogg_rofl , n)) * flux_epbalfact + x2o_o%rAttr(index_x2o_Foxx_rofi_16O, n) = (r2x_o%rAttr(index_r2x_Forr_rofi_16O , n) + & + g2x_o%rAttr(index_g2x_Fogg_rofi , n)) * flux_epbalfact + x2o_o%rAttr(index_x2o_Foxx_rofl_18O, n) = (r2x_o%rAttr(index_r2x_Forr_rofl_18O, n) + & + r2x_o%rAttr(index_r2x_Flrr_flood, n) + & + g2x_o%rAttr(index_g2x_Fogg_rofl , n)) * flux_epbalfact + x2o_o%rAttr(index_x2o_Foxx_rofi_18O, n) = (r2x_o%rAttr(index_r2x_Forr_rofi_18O , n) + & + g2x_o%rAttr(index_g2x_Fogg_rofi , n)) * flux_epbalfact + x2o_o%rAttr(index_x2o_Foxx_rofl_HDO, n) = (r2x_o%rAttr(index_r2x_Forr_rofl_HDO, n) + & + r2x_o%rAttr(index_r2x_Flrr_flood, n) + & + g2x_o%rAttr(index_g2x_Fogg_rofl , n)) * flux_epbalfact + x2o_o%rAttr(index_x2o_Foxx_rofi_HDO, n) = (r2x_o%rAttr(index_r2x_Forr_rofi_HDO , n) + & + g2x_o%rAttr(index_g2x_Fogg_rofi , n)) * flux_epbalfact + end if + + ! Derived: water isotopes total preciptiation and scaling + + if ( index_x2o_Faxa_snow_16O /= 0 )then + x2o_o%rAttr(index_x2o_Faxa_snow_16O ,n) = a2x_o%rAttr(index_a2x_Faxa_snowc_16O,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_snowl_16O,n) * afrac + x2o_o%rAttr(index_x2o_Faxa_rain_16O ,n) = a2x_o%rAttr(index_a2x_Faxa_rainc_16O,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_rainl_16O,n) * afrac + + x2o_o%rAttr(index_x2o_Faxa_snow_16O ,n) = x2o_o%rAttr(index_x2o_Faxa_snow_16O ,n) * flux_epbalfact + x2o_o%rAttr(index_x2o_Faxa_rain_16O ,n) = x2o_o%rAttr(index_x2o_Faxa_rain_16O ,n) * flux_epbalfact + + x2o_o%rAttr(index_x2o_Faxa_prec_16O ,n) = x2o_o%rAttr(index_x2o_Faxa_rain_16O ,n) + & + x2o_o%rAttr(index_x2o_Faxa_snow_16O ,n) + end if + + if ( index_x2o_Faxa_snow_18O /= 0 )then + x2o_o%rAttr(index_x2o_Faxa_snow_18O ,n) = a2x_o%rAttr(index_a2x_Faxa_snowc_18O,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_snowl_18O,n) * afrac + x2o_o%rAttr(index_x2o_Faxa_rain_18O ,n) = a2x_o%rAttr(index_a2x_Faxa_rainc_18O,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_rainl_18O,n) * afrac + + x2o_o%rAttr(index_x2o_Faxa_snow_18O ,n) = x2o_o%rAttr(index_x2o_Faxa_snow_18O ,n) * flux_epbalfact + x2o_o%rAttr(index_x2o_Faxa_rain_18O ,n) = x2o_o%rAttr(index_x2o_Faxa_rain_18O ,n) * flux_epbalfact + + x2o_o%rAttr(index_x2o_Faxa_prec_18O ,n) = x2o_o%rAttr(index_x2o_Faxa_rain_18O ,n) + & + x2o_o%rAttr(index_x2o_Faxa_snow_18O ,n) + end if + + if ( index_x2o_Faxa_snow_HDO /= 0 )then + x2o_o%rAttr(index_x2o_Faxa_snow_HDO ,n) = a2x_o%rAttr(index_a2x_Faxa_snowc_HDO,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_snowl_HDO,n) * afrac + x2o_o%rAttr(index_x2o_Faxa_rain_HDO ,n) = a2x_o%rAttr(index_a2x_Faxa_rainc_HDO,n) * afrac + & + a2x_o%rAttr(index_a2x_Faxa_rainl_HDO,n) * afrac + + x2o_o%rAttr(index_x2o_Faxa_snow_HDO ,n) = x2o_o%rAttr(index_x2o_Faxa_snow_HDO ,n) * flux_epbalfact + x2o_o%rAttr(index_x2o_Faxa_rain_HDO ,n) = x2o_o%rAttr(index_x2o_Faxa_rain_HDO ,n) * flux_epbalfact + + x2o_o%rAttr(index_x2o_Faxa_prec_HDO ,n) = x2o_o%rAttr(index_x2o_Faxa_rain_HDO ,n) + & + x2o_o%rAttr(index_x2o_Faxa_snow_HDO ,n) + end if + end do + + do ko = 1,noflds + !--- document merge --- + if (first_time) then + if (iindx(ko) > 0) then + if (imerge(ko)) then + mrgstr(ko) = trim(mrgstr(ko))//' + ifrac*i2x%'//trim(field_ice(iindx(ko))) + else + mrgstr(ko) = trim(mrgstr(ko))//' = ifrac*i2x%'//trim(field_ice(iindx(ko))) + end if + end if + if (aindx(ko) > 0) then + if (amerge(ko)) then + mrgstr(ko) = trim(mrgstr(ko))//' + afrac*a2x%'//trim(field_atm(aindx(ko))) + else + mrgstr(ko) = trim(mrgstr(ko))//' = afrac*a2x%'//trim(field_atm(aindx(ko))) + end if + end if + if (xindx(ko) > 0) then + if (xmerge(ko)) then + mrgstr(ko) = trim(mrgstr(ko))//' + afrac*xao%'//trim(field_xao(xindx(ko))) + else + mrgstr(ko) = trim(mrgstr(ko))//' = afrac*xao%'//trim(field_xao(xindx(ko))) + end if + end if + endif + + do n = 1,lsize + ifrac = fractions_o%rAttr(kif,n) + afrac = fractions_o%rAttr(kof,n) + frac_sum = ifrac + afrac + if ((frac_sum) /= 0._R8) then + ifrac = ifrac / (frac_sum) + afrac = afrac / (frac_sum) + endif + if (iindx(ko) > 0) then + if (imerge(ko)) then + x2o_o%rAttr(ko,n) = x2o_o%rAttr(ko,n) + i2x_o%rAttr(iindx(ko),n) * ifrac + else + x2o_o%rAttr(ko,n) = i2x_o%rAttr(iindx(ko),n) * ifrac + end if + end if + if (aindx(ko) > 0) then + if (amerge(ko)) then + x2o_o%rAttr(ko,n) = x2o_o%rAttr(ko,n) + a2x_o%rAttr(aindx(ko),n) * afrac + else + x2o_o%rAttr(ko,n) = a2x_o%rAttr(aindx(ko),n) * afrac + end if + end if + if (xindx(ko) > 0) then + if (xmerge(ko)) then + x2o_o%rAttr(ko,n) = x2o_o%rAttr(ko,n) + xao_o%rAttr(xindx(ko),n) * afrac + else + x2o_o%rAttr(ko,n) = xao_o%rAttr(xindx(ko),n) * afrac + end if + end if + end do + end do + + if (first_time) then + if (iamroot) then + write(logunit,'(A)') subname//' Summary:' + do ko = 1,noflds + write(logunit,'(A)') trim(mrgstr(ko)) + enddo + endif + deallocate(mrgstr) + deallocate(field_atm,itemc_atm) + deallocate(field_ocn,itemc_ocn) + deallocate(field_ice,itemc_ice) + deallocate(field_rof,itemc_rof) + deallocate(field_wav,itemc_wav) + deallocate(field_xao,itemc_xao) + endif + + first_time = .false. + + end subroutine prep_ocn_merge + + !================================================================================================ + subroutine prep_ocn_calc_a2x_ox(timer) !--------------------------------------------------------------- ! Arguments diff --git a/driver-moab/main/seq_diag_mct.F90 b/driver-moab/main/seq_diag_mct.F90 index 8534008f9be7..9938119f4b0f 100644 --- a/driver-moab/main/seq_diag_mct.F90 +++ b/driver-moab/main/seq_diag_mct.F90 @@ -1873,7 +1873,7 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, & if (plev > 0) then ! ---- doprint ---- doprint ---- doprint ---- - + if (.not.sumdone) then if (do_bgc_budg) then @@ -2240,7 +2240,7 @@ SUBROUTINE seq_diag_print_mct(EClock, stop_alarm, do_bgc_budg, & ! ---- doprint ---- doprint ---- doprint ---- if (do_bgc_budg) then - call seq_diagBGC_print_mct(EClock, ip, plev) + call seq_diagBGC_print_mct(EClock, ip, plev) endif endif ! plev > 0 diff --git a/driver-moab/shr/seq_flds_mod.F90 b/driver-moab/shr/seq_flds_mod.F90 index 99d6580183e9..c807fb6ef663 100644 --- a/driver-moab/shr/seq_flds_mod.F90 +++ b/driver-moab/shr/seq_flds_mod.F90 @@ -836,7 +836,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'kg m-3' attname = 'Sa_dens' call metadata_set(attname, longname, stdname, units) - + ! UoverN for use by topounits call seq_flds_add(a2x_states,"Sa_uovern") call seq_flds_add(x2l_states,"Sa_uovern") @@ -1668,6 +1668,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) ! Ocean land ice frazil production call seq_flds_add(o2x_fluxes,"Foxo_frazil_li") + call seq_flds_add(x2g_fluxes,"Foxo_frazil_li") longname = 'Ocean land ice frazil production' stdname = 'ocean_land_ice_frazil_ice_production' units = 'kg m-2 s-1' @@ -1676,6 +1677,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) ! Water flux from ice shelf melt call seq_flds_add(o2x_fluxes,"Foxo_ismw") + call seq_flds_add(x2g_fluxes,"Foxo_ismw") longname = 'Water flux due to basal melting of ice shelves' stdname = 'basal_iceshelf_melt_flux' units = 'kg m-2 s-1' @@ -2464,7 +2466,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'kg m-2 s-1' attname = 'Flrr_supply' call metadata_set(attname, longname, stdname, units) - + call seq_flds_add(r2x_fluxes,'Flrr_deficit') call seq_flds_add(x2l_fluxes,'Flrr_deficit') longname = 'River model supply deficit' @@ -2818,12 +2820,12 @@ subroutine seq_flds_set(nmlfile, ID, infodata) !---------------------------- ! lnd->iac, iac->lnd, iac->atm !---------------------------- - + ! lnd/iac coupling needs one field in each class per pft ! Note that this ends up as 17*4=68 coupled fields... ! also send harvest fraction from iac to land ! use the same loop and index string - + ! these two variables are hardcoded above because there is not an appropriate ! namelist to put them in: iac_npft and iac_nharvest @@ -2842,7 +2844,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'Sl_hr_pft' // pftstr attname = trim(attname) call metadata_set(attname, longname, stdname, units) - + if(add_iac_to_cplstate)call seq_flds_add(l2x_states,'Sl_npp_pft' // pftstr) call seq_flds_add(x2z_states,'Sl_npp_pft' // pftstr) longname = 'Net primary production for pft ' // pftstr @@ -2850,7 +2852,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) units = 'gC/m^2/s' attname = 'Sl_npp_pft' // pftstr call metadata_set(attname, longname, stdname, units) - + ! Review if(add_iac_to_cplstate)call seq_flds_add(l2x_states,'Sl_pftwgt_pft' //pftstr) call seq_flds_add(x2z_states,'Sl_pftwgt_pft' //pftstr) @@ -2885,7 +2887,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'Sz_pct_pft_prev' //pftstr attname = trim(attname) call metadata_set(attname, longname, stdname, units) - + ! send the harvest data also, these are for model year if (i <= iac_nharvest) then call seq_flds_add(z2x_states,trim('Sz_harvest_frac' //pftstr)) @@ -3230,54 +3232,24 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call set_glc_elevclass_field(name, attname, longname, stdname, units, x2l_fluxes_from_glc, & additional_list = .true.) - name = 'So_blt' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf boundary layer ocean temperature' - stdname = 'Ice_shelf_boundary_layer_ocean_temperature' - units = 'C' - attname = 'So_blt' - call metadata_set(attname, longname, stdname, units) - - name = 'So_bls' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf boundary layer ocean salinity' - stdname = 'Ice_shelf_boundary_layer_ocean_salinity' - units = 'psu' - attname = 'So_bls' - call metadata_set(attname, longname, stdname, units) - - name = 'So_htv' - call seq_flds_add(o2x_states,trim(name)) - call seq_flds_add(x2g_states,trim(name)) - call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf ocean heat transfer velocity' - stdname = 'Ice_shelf_ocean_heat_transfer_velocity' - units = 'm/s' - attname = 'So_htv' - call metadata_set(attname, longname, stdname, units) - - name = 'So_stv' + name = 'So_tfrz_isf' call seq_flds_add(o2x_states,trim(name)) call seq_flds_add(x2g_states,trim(name)) call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ice shelf ocean salinity transfer velocity' - stdname = 'Ice_shelf_ocean_salinity_transfer_velocity' - units = 'm/s' - attname = 'So_stv' + longname = 'Ice shelf-ocean freezing temperature' + stdname = 'Ice_shelf_ocean_freezing_temperature' + units = 'K' + attname = 'So_tfrz_isf' call metadata_set(attname, longname, stdname, units) - name = 'So_rhoeff' + name = 'So_liflfrac' call seq_flds_add(o2x_states,trim(name)) call seq_flds_add(x2g_states,trim(name)) call seq_flds_add(x2g_shelf_states_from_ocn,trim(name)) - longname = 'Ocean effective pressure' - stdname = 'Ocean_effective_pressure' - units = 'Pa' - attname = 'So_rhoeff' + longname = 'Floating ice shelf area fraction' + stdname = 'floating_ice_shelf_area_fraction' + units = '1' + attname = name call metadata_set(attname, longname, stdname, units) if ((flds_tf) .and. (glc_nzoc > 0)) then @@ -3309,42 +3281,6 @@ subroutine seq_flds_set(nmlfile, ID, infodata) additional_list = .true.) end if - name = 'Fogx_qicelo' - call seq_flds_add(g2x_fluxes,trim(name)) - call seq_flds_add(x2o_fluxes,trim(name)) - longname = 'Subshelf liquid flux for ocean' - stdname = 'Subshelf_liquid_flux_for_ocean' - units = 'kg m-2 s-1' - attname = 'Fogx_qicelo' - call metadata_set(attname, longname, stdname, units) - - name = 'Fogx_qiceho' - call seq_flds_add(g2x_fluxes,trim(name)) - call seq_flds_add(x2o_fluxes,trim(name)) - longname = 'Subshelf heat flux for the ocean' - stdname = 'Subshelf_heat_flux_for_the_ocean' - units = 'W m-2' - attname = 'Fogx_qiceho' - call metadata_set(attname, longname, stdname, units) - - name = 'Sg_blit' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Boundary layer interface temperature for ocean' - stdname = 'Boundary_layer_interface_temperature_for_ocean' - units = 'C' - attname = 'Sg_blit' - call metadata_set(attname, longname, stdname, units) - - name = 'Sg_blis' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Boundary layer interface salinity for ocean' - stdname = 'Boundary_layer_interface_salinity_for_ocean' - units = 'psu' - attname = 'Sg_blis' - call metadata_set(attname, longname, stdname, units) - name = 'Sg_lithop' call seq_flds_add(g2x_states,trim(name)) call seq_flds_add(x2o_states,trim(name)) @@ -3372,32 +3308,6 @@ subroutine seq_flds_set(nmlfile, ID, infodata) attname = 'Sg_icemask_floating' call metadata_set(attname, longname, stdname, units) - name = 'Sg_tbot' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Bottom layer ice temperature' - stdname = 'Bottom_layer_ice_temperature' - units = 'C' - attname = 'Sg_tbot' - call metadata_set(attname, longname, stdname, units) - - name = 'Sg_dztbot' - call seq_flds_add(g2x_states,trim(name)) - call seq_flds_add(x2o_states,trim(name)) - longname = 'Bottom layer ice layer half thickness' - stdname = 'Bottom_layer_ice_layer_half_thickness' - units = 'm' - attname = 'Sg_dztbot' - call metadata_set(attname, longname, stdname, units) - - name = 'Fogx_qiceli' - call seq_flds_add(x2g_fluxes,trim(name)) - longname = 'Subshelf mass flux for ice sheet' - stdname = 'Subshelf_mass_flux_for_ice_sheet' - units = 'kg m-2 s-1' - attname = 'Fogx_qiceli' - call metadata_set(attname, longname, stdname, units) - name = 'Fogx_qicehi' call seq_flds_add(x2g_fluxes,trim(name)) longname = 'Subshelf heat flux for ice sheet' @@ -4386,7 +4296,7 @@ subroutine seq_flds_set(nmlfile, ID, infodata) call catFields(seq_flds_z2x_fields, seq_flds_z2x_states, seq_flds_z2x_fluxes) call catFields(seq_flds_x2z_fields, seq_flds_x2z_states, seq_flds_x2z_fluxes) call catFields(seq_flds_o2x_fields_to_lnd, seq_flds_o2x_states_to_lnd, seq_flds_o2x_fluxes_to_lnd) - + if (seq_comm_iamroot(ID)) then write(logunit,*) subname//': seq_flds_dom_fields= ',trim(seq_flds_dom_fields) write(logunit,*) subname//': seq_flds_xao_fields= ',trim(seq_flds_xao_fields) @@ -4685,7 +4595,7 @@ subroutine set_glc_zocnclass_field(name, attname, longname, stdname, units, fiel integer :: num character(len= 16) :: cnum logical :: l_additional_list ! local version of the optional additional_list argument - + l_additional_list = .false. if (present(additional_list)) then l_additional_list = additional_list