Skip to content

Commit

Permalink
Update temperature and waterFrac in RK loop, but not enthalpy
Browse files Browse the repository at this point in the history 
Update temperature and waterFrac in RK loop, but not enthalpy. Updating
enthalpy resulted in negative temperatures, which caused an error in li_thermal.
Updating temperature and waterFrac instead after enthalpy is advected allows
enthalpy tests in full_integration suite to pass execution and validation.
  • Loading branch information
@trhille
trhille committed Oct 9, 2023
1 parent 35977fc commit 4005761
Showing 1 changed file with 15 additions and 33 deletions.
48 changes: 15 additions & 33 deletions components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe_rk.F
View file
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
integer :: rkStage, iCell, iTracer, k
real (kind=RKIND), dimension(:,:), pointer :: layerThickness, &
temperature, &
enthalpy, &
waterFrac
real (kind=RKIND), dimension(:), pointer :: thickness, damage, passiveTracer2d
real (kind=RKIND), dimension(:), pointer :: sfcMassBalApplied, &
Expand All @@ -141,7 +140,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
real (kind=RKIND), dimension(:,:), allocatable :: layerThicknessPrev, &
layerThicknessTmp, &
temperaturePrev, &
enthalpyPrev, &
waterFracPrev, &
passiveTracer3d, &
passiveTracer3dPrev, &
damage3d, &
Expand Down Expand Up @@ -186,7 +185,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
call mpas_pool_get_array(geometryPool, 'groundedToFloatingThickness', groundedToFloatingThickness)

allocate(temperaturePrev(nVertLevels, nCells+1))
allocate(enthalpyPrev(nVertLevels, nCells+1))
allocate(waterFracPrev(nVertLevels, nCells+1))
allocate(layerThicknessPrev(nVertLevels, nCells+1))
allocate(layerThicknessTmp(nVertLevels, nCells+1))
allocate(damage3dPrev(nVertLevels, nCells+1))
Expand All @@ -203,7 +202,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
allocate(fluxAcrossGroundingLineAccum(nEdges+1))

temperaturePrev(:,:) = 0.0_RKIND
enthalpyPrev(:,:) = 0.0_RKIND
layerThicknessPrev(:,:) = 0.0_RKIND
layerThicknessTmp(:,:) = 0.0_RKIND
damage3dPrev(:,:) = 0.0_RKIND
Expand Down Expand Up @@ -272,11 +270,10 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
call mpas_pool_get_array(thermalPool, 'temperature', temperature)
call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy)
call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac)
! Save relevant fields before RK loop, to be used in update at the end
temperaturePrev(:,:) = temperature(:,:)
enthalpyPrev(:,:) = enthalpy(:,:)
waterFracPrev(:,:) = waterFrac(:,:)
layerThicknessPrev(:,:) = layerThickness(:,:)
cellMaskPrev(:) = cellMask(:)
do k = 1, nVertLevels
Expand Down Expand Up @@ -384,40 +381,26 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
call mpas_pool_get_array(geometryPool, 'damage', damage)
call mpas_pool_get_array(thermalPool, 'temperature', temperature)
call mpas_pool_get_array(thermalPool, 'enthalpy', enthalpy)
call mpas_pool_get_array(thermalPool, 'waterFrac', waterFrac)
layerThicknessTmp(:,:) = layerThickness(:,:)
layerThickness(:,:) = rkSSPweights(rkStage) * layerThicknessPrev(:,:) + &
(1.0_RKIND - rkSSPweights(rkStage)) * layerThickness(:,:)
thickness = sum(layerThickness, 1)
if (trim(config_thermal_solver) == 'enthalpy') then
where (layerThickness(:,:) > 0.0_RKIND)
enthalpy(:,:) = ( rkSSPweights(rkStage) * enthalpyPrev(:,:) * layerThicknessPrev(:,:) + &
(1.0_RKIND - rkSSPweights(rkStage)) * enthalpy(:,:) * &
layerThicknessTmp(:,:) ) / layerThickness(:,:)
! elsewhere, what?
end where
! given the enthalpy, compute the temperature and water fraction
if (trim(config_thermal_solver) .ne. 'none') then
do iCell = 1, nCells
call li_enthalpy_to_temperature_kelvin(&
layerCenterSigma, &
thickness(iCell), &
enthalpy(:,iCell), &
temperature(:,iCell), &
waterFrac(:,iCell))
do k = 1, nVertLevels
if (layerThickness(k,iCell) > 0.0_RKIND) then
temperature(k,iCell) = ( rkSSPweights(rkStage) * temperaturePrev(k,iCell) * layerThicknessPrev(k,iCell) + &
(1.0_RKIND - rkSSPweights(rkStage)) * temperature(k,iCell) * &
layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell)
waterFrac(k,iCell) = ( rkSSPweights(rkStage) * waterFracPrev(k,iCell) * layerThicknessPrev(k,iCell) + &
(1.0_RKIND - rkSSPweights(rkStage)) * waterFrac(k,iCell) * &
layerThicknessTmp(k,iCell) ) / layerThickness(k,iCell)
endif
enddo
enddo
elseif (trim(config_thermal_solver) == 'temperature') then
where (layerThickness(:,:) > 0.0_RKIND)
temperature(:,:) = ( rkSSPweights(rkStage) * temperaturePrev(:,:) * layerThicknessPrev(:,:) + &
(1.0_RKIND - rkSSPweights(rkStage)) * temperature(:,:) * &
layerThicknessTmp(:,:) ) / layerThickness(:,:)
! elsewhere, what?
end where
endif
if (config_calculate_damage) then
Expand Down Expand Up @@ -486,7 +469,6 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
call mpas_dmpar_field_halo_exch(domain, 'thickness')
call mpas_dmpar_field_halo_exch(domain, 'temperature')
call mpas_dmpar_field_halo_exch(domain, 'waterFrac')
call mpas_dmpar_field_halo_exch(domain, 'enthalpy')
call mpas_dmpar_field_halo_exch(domain, 'damage')
call mpas_dmpar_field_halo_exch(domain, 'passiveTracer2d')
Expand Down Expand Up @@ -580,7 +562,7 @@ subroutine li_time_integrator_forwardeuler_rungekutta(domain, err)
endif

deallocate(temperaturePrev)
deallocate(enthalpyPrev)
deallocate(waterFracPrev)
deallocate(layerThicknessPrev)
deallocate(layerThicknessTmp)
deallocate(damage3dPrev)
Expand Down

0 comments on commit 4005761

Please sign in to comment.