diff --git a/components/mpas-albany-landice/src/Registry.xml b/components/mpas-albany-landice/src/Registry.xml index 14b8527e6225..f716f398f994 100644 --- a/components/mpas-albany-landice/src/Registry.xml +++ b/components/mpas-albany-landice/src/Registry.xml @@ -42,7 +42,7 @@ - + possible_values="'fo', 'fct', 'none'" + /> + + + @@ -741,6 +753,7 @@ + @@ -850,6 +863,7 @@ + @@ -1270,6 +1284,9 @@ is the value of that variable from the *previous* time level! + diff --git a/components/mpas-albany-landice/src/mode_forward/Makefile b/components/mpas-albany-landice/src/mode_forward/Makefile index 71cae3b4d1d1..c2c9d880ca6b 100644 --- a/components/mpas-albany-landice/src/mode_forward/Makefile +++ b/components/mpas-albany-landice/src/mode_forward/Makefile @@ -7,6 +7,8 @@ OBJS = mpas_li_core.o \ mpas_li_time_integration_fe.o \ mpas_li_diagnostic_vars.o \ mpas_li_advection.o \ + mpas_li_advection_fct_shared.o \ + mpas_li_advection_fct.o \ mpas_li_calving.o \ mpas_li_statistics.o \ mpas_li_velocity.o \ @@ -45,8 +47,12 @@ mpas_li_time_integration_fe.o: mpas_li_advection.o \ mpas_li_bedtopo.o mpas_li_advection.o: mpas_li_thermal.o \ + mpas_li_advection_fct.o \ + mpas_li_advection_fct_shared.o \ mpas_li_diagnostic_vars.o +mpas_li_advection_fct.o: mpas_li_advection_fct_shared.o + mpas_li_calving.o: mpas_li_thermal.o \ mpas_li_advection.o diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F index 4b84ab7206fd..df8f1b95ba8a 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F @@ -27,6 +27,8 @@ module li_advection use li_setup use li_mask use li_constants + use li_config + use li_mesh implicit none private @@ -82,7 +84,8 @@ subroutine li_advection_thickness_tracers(& advectTracersIn) use li_setup, only: li_calculate_layerThickness - + use li_tracer_advection_fct, only: li_tracer_advection_fct_tend + use li_tracer_advection_fct_shared !----------------------------------------------------------------- ! ! input variables @@ -162,7 +165,8 @@ subroutine li_advection_thickness_tracers(& layerNormalVelocity, & ! normal component of velocity on cell edges at layer midpoints layerThickness, & ! thickness of each layer layerThicknessOld, & ! old layer thickness - layerThicknessEdge ! layer thickness on upstream edge of cell + layerThicknessEdge, & ! layer thickness on upstream edge of cell + normalVelocity ! horizontal velocity on interfaces real (kind=RKIND), dimension(:,:,:), pointer :: & advectedTracers, & ! values of advected tracers @@ -187,6 +191,10 @@ subroutine li_advection_thickness_tracers(& cellMaskTemporaryField, & ! scratch field containing old values of cellMask thermalCellMaskField + ! Allocatable arrays need for flux-corrected transport advection + real (kind=RKIND), dimension(:,:,:), allocatable :: tend + real (kind=RKIND), dimension(:,:,:), allocatable :: activeTracerHorizontalAdvectionEdgeFlux + integer, dimension(:), pointer :: & cellMask, & ! integer bitmask for cells edgeMask, & ! integer bitmask for edges @@ -208,15 +216,22 @@ subroutine li_advection_thickness_tracers(& config_ice_density, & ! rhoi config_ocean_density, & ! rhoo config_sea_level, & ! sea level relative to z = 0 - config_thermal_thickness + config_thermal_thickness, & + config_advection_coef_3rd_order + + integer, pointer :: config_num_halos logical :: advectTracers ! if true, then advect tracers as well as thickness integer :: iCell1, iCell2, theGroundedCell + integer, parameter :: horizAdvOrder = 4 + integer, parameter :: vertAdvOrder = 1 real(kind=RKIND) :: GLfluxSign, thicknessFluxEdge integer :: err_tmp + integer :: nTracers + integer, pointer :: maxEdges2 !WHL - debug integer, dimension(:), pointer :: indexToCellID, indexToEdgeID @@ -277,6 +292,7 @@ subroutine li_advection_thickness_tracers(& ! get arrays from the velocity pool call mpas_pool_get_array(velocityPool, 'layerNormalVelocity', layerNormalVelocity) + call mpas_pool_get_array(velocityPool, 'normalVelocity', normalVelocity) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLine', fluxAcrossGroundingLine) call mpas_pool_get_array(velocityPool, 'fluxAcrossGroundingLineOnCells', fluxAcrossGroundingLineOnCells) @@ -295,11 +311,12 @@ subroutine li_advection_thickness_tracers(& call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) call mpas_pool_get_config(liConfigs, 'config_ocean_density', config_ocean_density) call mpas_pool_get_config(liConfigs, 'config_sea_level', config_sea_level) + call mpas_pool_get_config(liConfigs, 'config_num_halos', config_num_halos) - !WHL - debug call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_dimension(meshPool, 'nEdges', nEdges) call mpas_pool_get_dimension(meshPool, 'nVertLevels', nVertLevels) + call mpas_pool_get_dimension(meshPool, 'maxEdges2', maxEdges2) call mpas_pool_get_array(meshPool, 'edgesOnCell', edgesOnCell) call mpas_pool_get_array(meshPool, 'nEdgesOnCell', nEdgesOnCell) call mpas_pool_get_array(meshPool, 'indexToCellID', indexToCellID) @@ -339,6 +356,9 @@ subroutine li_advection_thickness_tracers(& call mpas_allocate_scratch_field(thermalCellMaskField, .true.) thermalCellMask => thermalCellMaskField % array + ! define max number of advection cells as max number of edges*2 + nAdvCellsMax = maxEdges2 + ! given the old thickness, compute the thickness in each layer call li_calculate_layerThickness(meshPool, thickness, layerThickness) @@ -353,9 +373,8 @@ subroutine li_advection_thickness_tracers(& ! Horizontal transport of thickness and tracers !----------------------------------------------------------------- - select case (trim(config_thickness_advection)) - - case ('fo') + if ( (trim(config_thickness_advection) .eq. 'fo') .or. & + (trim(config_thickness_advection) .eq. 'fct') ) then if (config_print_thickness_advection_info) then call mpas_log_write('Using first-order upwind for thickness advection') @@ -367,7 +386,7 @@ subroutine li_advection_thickness_tracers(& endif - if (advectTracers) then + if ( (advectTracers) .or. (trim(config_thickness_advection) .eq. 'fct') ) then ! Copy the old tracer values into the advectedTracers array. ! This requires setting a tracer pointer to either temperature or enthalpy, @@ -381,7 +400,8 @@ subroutine li_advection_thickness_tracers(& thermalPool, & advectedTracers, & surfaceTracers, & - basalTracers) + basalTracers, & + nTracers) else ! no tracer advection; pass a tracer array full of zeroes @@ -389,6 +409,14 @@ subroutine li_advection_thickness_tracers(& endif + if ((trim(config_thickness_advection) .eq. 'fct') .or. & + (trim(config_tracer_advection) .eq. 'fct' ) ) then + allocate(tend(nTracers,nVertLevels,nCells+1)) + tend(:,:,:) = 0.0_RKIND + allocate(activeTracerHorizontalAdvectionEdgeFlux(nTracers,nVertLevels+1,nEdges+1)) + activeTracerHorizontalAdvectionEdgeFlux(:,:,:) = 0.0_RKIND + endif + ! Transport thickness and tracers using a first-order upwind scheme ! Note: For the enthalpy scheme, temperature and waterFrac are the primary prognostic ! variables to be updated, but enthalpy is the advected tracer (for reasons of @@ -402,8 +430,10 @@ subroutine li_advection_thickness_tracers(& intArgs=(/iCell, indexToCellID(iCell)/), realArgs=(/thickness(iCell), advectedTracers(1,1,iCell)/)) do iEdgeOnCell = 1, nEdgesOnCell(iCell) iEdge = edgesOnCell(iEdgeOnCell,iCell) - call mpas_log_write('iEdge (local=$i, global=$i), layerNormalVelocity=$r:', & - intArgs=(/iEdge, indexToEdgeID(iEdge)/), realArgs=(/layerNormalVelocity(1,iEdge)/)) + call mpas_log_write('iEdge (local=$i, global=$i), & + layerNormalVelocity=$r:', & + intArgs=(/iEdge, indexToEdgeID(iEdge)/), & + realArgs=(/layerNormalVelocity(1,iEdge)/)) enddo endif enddo @@ -414,18 +444,78 @@ subroutine li_advection_thickness_tracers(& advectedTracersOld(:,:,:) = advectedTracers(:,:,:) ! compute new values of layer thickness and tracers - call advect_thickness_tracers_upwind(& - dt, & - meshPool, & - layerNormalVelocity, & - layerThicknessEdge, & - layerThicknessOld, & - advectedTracersOld, & - cellMask, & - edgeMask, & - layerThickness, & - advectedTracers, & - err) + if ((trim(config_thickness_advection) .eq. 'fo') .and. & + ( (trim(config_tracer_advection) .eq. 'fo') .or. & + (trim(config_tracer_advection) .eq. 'none') ) ) then + call advect_thickness_tracers_upwind(& + dt, & + meshPool, & + layerNormalVelocity, & + layerThicknessEdge, & + layerThicknessOld, & + advectedTracersOld, & + cellMask, & + edgeMask, & + layerThickness, & + advectedTracers, & + err) + elseif ((trim(config_thickness_advection) .eq. 'fo') .and. & + trim(config_tracer_advection) .eq. 'fct') then + ! TODO: Decide what to do about w (vertical vel). + ! For now, just setting it to 0 to try to get this to compile. + call advect_thickness_tracers_upwind(& + dt, & + meshPool, & + layerNormalVelocity, & + layerThicknessEdge, & + layerThicknessOld, & + advectedTracersOld, & + cellMask, & + edgeMask, & + layerThickness, & + advectedTracers, & + err) + ! Reset tracers after fo advection for fct advection + advectedTracers(:,:,:) = advectedTracersOld(:,:,:) + + ! Pass FO upwind normalThicknessFlux (layerThicknessEdge * layerNormalVelocity) + ! to fct tracer routine + call li_tracer_advection_fct_tend(& + tend, advectedTracers, layerThicknessOld, & + layerThicknessEdge * layerNormalVelocity, 0 * normalVelocity, dt, & + nTracers, activeTracerHorizontalAdvectionEdgeFlux, computeBudgets=.false.)!{{{ + elseif (trim(config_thickness_advection) .eq. 'fct') then + ! Call fct routine for thickness first, and use activeTracerHorizontalAdvectionEdgeFlux + ! returned by that call as normalThicknessFlux for call to tracer fct + call li_tracer_advection_fct_tend(& + tend(nTracers:,:,:), advectedTracers(nTracers:,:,:), layerThicknessOld * 0.0_RKIND + 1.0_RKIND, & + layerNormalVelocity, 0.0_RKIND * normalVelocity, dt, & + 1, activeTracerHorizontalAdvectionEdgeFlux(nTracers:,:,:), computeBudgets=.false.) + ! layerThickness is last tracer. However, for some reason + ! this: layerThickness(:,:) = advectedTracers(nTracers,:,:) does not conserve mass! + ! This does conserve mass: + layerThickness(:,:) = layerThickness(:,:) + tend(nTracers,:,:) * dt + + if (trim(config_tracer_advection) .eq. 'fct') then + ! Call fct for tracers, using activeTracerHorizontalAdvectionEdgeFlux + ! from fct thickness advection as normalThicknessFlux + call li_tracer_advection_fct_tend(& + tend(1:nTracers-1,:,:), advectedTracers(1:nTracers-1,:,:), layerThicknessOld, & + activeTracerHorizontalAdvectionEdgeFlux(nTracers,:,:), 0.0_RKIND * normalVelocity, dt, & + nTracers-1, computeBudgets=.false.) + elseif (trim(config_tracer_advection) .eq. 'none') then + ! do nothing + else + err_tmp = 1 + call mpas_log_write(trim(config_tracer_advection) // & + ' tracer advection is not currently supported with fct thickness advection.', MPAS_LOG_ERR) + endif + else + err_tmp = 1 + call mpas_log_write("config_thickness_advection = " // trim(config_thickness_advection) // & + ", config_tracer_advection = " // trim(config_tracer_advection) // & + " is not a supported combination.", MPAS_LOG_ERR) + endif if (config_print_thickness_advection_info) then do iCell = 1, nCells @@ -601,23 +691,35 @@ subroutine li_advection_thickness_tracers(& geometryPool, & thermalPool, & advectedTracers) - endif - case ('none') + elseif (trim(config_thickness_advection) .eq. 'none') then ! Do nothing - case default + else call mpas_log_write(trim(config_thickness_advection) // ' is not a valid option for thickness/tracer advection.', & MPAS_LOG_ERR) err_tmp = 1 - end select + end if ! config_thickness_advection err = ior(err,err_tmp) + ! Deallocate arrays for fct + if ( (trim(config_thickness_advection) .eq. 'fct') .or. & + (trim(config_tracer_advection) .eq. 'fct') ) then + deallocate( nAdvCellsForEdge, & + advCellsForEdge, & + advCoefs, & + advCoefs3rd, & + advMaskHighOrder, & + advMask2ndOrder, & + tend, & + activeTracerHorizontalAdvectionEdgeFlux) + endif + ! clean up call mpas_deallocate_scratch_field(advectedTracersField, .true.) call mpas_deallocate_scratch_field(advectedTracersOldField, .true.) @@ -981,10 +1083,12 @@ subroutine tracer_setup(& thermalPool, & advectedTracers, & surfaceTracers, & - basalTracers) + basalTracers, & + nTracers) use li_thermal, only: li_temperature_to_enthalpy_kelvin - + use li_tracer_advection_fct_shared + use li_tracer_advection_fct !----------------------------------------------------------------- ! ! input variables @@ -1015,6 +1119,7 @@ subroutine tracer_setup(& basalTracers ! tracer values for new ice at lower surface ! dimension 1 = maxTracers, 2 = nCells + integer, intent(out) :: nTracers !----------------------------------------------------------------- ! ! local variables @@ -1022,7 +1127,8 @@ subroutine tracer_setup(& !----------------------------------------------------------------- character (len=StrKIND), pointer :: & - config_thermal_solver + config_thermal_solver, & + config_thickness_advection logical, pointer :: & config_calculate_damage @@ -1036,7 +1142,8 @@ subroutine tracer_setup(& real (kind=RKIND), dimension(:,:), pointer :: & temperature, & ! interior ice temperature waterFrac, & ! interior water fraction - enthalpy ! interior ice enthalpy + enthalpy, & ! interior ice enthalpy + layerThickness real (kind=RKIND), dimension(:), pointer :: & surfaceAirTemperature, & ! surface air temperature @@ -1049,12 +1156,17 @@ subroutine tracer_setup(& thickness ! ice thickness real (kind=RKIND), dimension(:), pointer :: & - damage ! damage + damage, & ! damage + passiveTracer2d real (kind=RKIND), pointer :: & config_ice_density ! ice density - integer :: iCell, iTracer, k + integer :: iCell, iTracer, k, err, err1, err2 + + err = 0 + err1 = 0 + err2 = 0 ! get dimensions call mpas_pool_get_dimension(meshPool, 'nCells', nCells) @@ -1065,7 +1177,9 @@ subroutine tracer_setup(& ! get arrays from geometry pool call mpas_pool_get_array(geometryPool, 'thickness', thickness) + call mpas_pool_get_array(geometryPool, 'layerThickness', layerThickness) call mpas_pool_get_array(geometryPool, 'damage', damage) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) ! get arrays from thermal pool call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -1078,6 +1192,7 @@ subroutine tracer_setup(& call mpas_pool_get_config(liConfigs, 'config_thermal_solver', config_thermal_solver) call mpas_pool_get_config(liConfigs, 'config_calculate_damage', config_calculate_damage) call mpas_pool_get_config(liConfigs, 'config_ice_density', config_ice_density) + call mpas_pool_get_config(liConfigs, 'config_thickness_advection', config_thickness_advection) ! initialize advectedTracers(:,:,:) = 0.0_RKIND @@ -1155,8 +1270,41 @@ subroutine tracer_setup(& endif + ! always advect passiveTracer2d + iTracer = iTracer + 1 + ! copy 2d passiveTracer2d to all vertical levels in the tracer array + do k = 1, nVertLevels + advectedTracers(iTracer,k,:) = passiveTracer2d(:) + enddo + surfaceTracers(iTracer,:) = 0.0_RKIND + basalTracers(iTracer,:) = 0.0_RKIND + + if (trim(config_thickness_advection) == 'fct') then + ! increment the tracer index + iTracer = iTracer + 1 + + ! copy 2d thickness to all vertical levels in the tracer array + advectedTracers(iTracer,:,:) = layerThickness(:,:) + surfaceTracers(iTracer,:) = 0.0_RKIND + basalTracers(iTracer,:) = 0.0_RKIND + + endif + + nTracers = iTracer ! Note: Other tracers (e.g., ice age) can be added here as needed. ! May need to increase maxTracers in the Registry. + if ( (trim(config_tracer_advection) == 'fct') .or. & + (trim(config_thickness_advection) == 'fct') ) then + call li_tracer_advection_fct_shared_init(geometryPool, err1) + call li_tracer_advection_fct_init(err2) + + if (err1 /= 0 .or. err2 /= 0) then + err = 1 + call mpas_log_write( & + 'Error encountered during fct tracer advection init', & + MPAS_LOG_ERR, masterOnly=.true.) + endif + endif end subroutine tracer_setup @@ -1228,7 +1376,8 @@ subroutine tracer_finish(& thickness ! ice thickness real (kind=RKIND), dimension(:), pointer :: & - damage ! damage + damage, & ! damage + passiveTracer2d real (kind=RKIND), dimension(:,:), pointer :: & temperature, & ! interior ice temperature @@ -1249,6 +1398,7 @@ subroutine tracer_finish(& ! get arrays from geometry pool call mpas_pool_get_array(geometryPool, 'thickness', thickness) call mpas_pool_get_array(geometryPool, 'damage', damage) + call mpas_pool_get_array(geometryPool, 'passiveTracer2d', passiveTracer2d) ! get arrays from thermal pool call mpas_pool_get_array(thermalPool, 'temperature', temperature) @@ -1299,6 +1449,12 @@ subroutine tracer_finish(& enddo endif + ! Always advect passiveTracer2d + iTracer = iTracer + 1 + do iCell = 1, nCellsSolve + passiveTracer2d(iCell) = sum(advectedTracers(iTracer,:,iCell) * layerThicknessFractions) ! thickness-weighted average + enddo + ! Note: Other tracers (e.g., ice age) can be added here as needed. ! May need to increase maxTracers in the Registry. diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F new file mode 100644 index 000000000000..7dccfdaa8b88 --- /dev/null +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct.F @@ -0,0 +1,1108 @@ +! and the University Corporation for Atmospheric Research (UCAR). +! Copyright (c) 2013, Los Alamos National Security, LLC (LANS) +! +! 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.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! tracer_advection_mono +! +!> \brief MPAS monotonic tracer advection with FCT +!> \author Mark Petersen, David Lee, Doug Jacobsen, Phil Jones +!> \date October 2017, updated May 2019 +!> \details +!> This module contains routines for monotonic advection of tracers +!> using a Flux Corrected Transport (FCT) algorithm +! +!----------------------------------------------------------------------- + +module li_tracer_advection_fct + + ! module includes +#ifdef _ADV_TIMERS + use mpas_timer +#endif + use mpas_kind_types + + use li_config + use li_mesh +! use ocn_diagnostics_variables + use li_tracer_advection_fct_shared +! use tracer_advection_vert + + implicit none + private + save + + ! module private variables + real (kind=RKIND) :: & + coef3rdOrder !< high-order horizontal coefficient + logical :: & + monotonicityCheck !< flag to check monotonicity + + ! public method interfaces + public :: li_tracer_advection_fct_tend, & + li_tracer_advection_fct_init + +!********************************************************************** + + contains + +!********************************************************************** +! +! routine li_tracer_advection_fct_tend +! +!> \brief MPAS monotonic tracer horizontal advection tendency with FCT +!> \author Mark Petersen, David Lee, Doug Jacobsen, Phil Jones +!> \date October 2017, updated May 2019 +!> \details +!> This routine computes the monotonic tracer horizontal advection +!> tendency using a flux-corrected transport (FCT) algorithm. +! +!----------------------------------------------------------------------- + + subroutine li_tracer_advection_fct_tend( & + tend, tracers, layerThickness, & + normalThicknessFlux, w, dt, & + nTracers, & + activeTracerHorizontalAdvectionEdgeFlux, & + computeBudgets)!{{{ + use li_mesh + !----------------------------------------------------------------- + ! Input/Output parameters + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + tend !< [inout] Tracer tendency to which advection added + real (kind=RKIND), dimension(:,:,:), intent(inout), optional :: & + activeTracerHorizontalAdvectionEdgeFlux !< [inout] used to compute higher order normalThicknessFlux + !----------------------------------------------------------------- + ! Input parameters + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(inout) :: & + tracers !< [in] Current tracer values + + real (kind=RKIND), dimension(:,:), intent(in) :: & + layerThickness, &!< [in] Thickness + normalThicknessFlux, &!< [in] Thichness weighted velocitiy + w !< [in] Vertical velocity + + real (kind=RKIND), intent(in) :: & + dt !< [in] Timestep + + integer, intent(in) :: nTracers + + logical, intent(in) :: & + computeBudgets !< [in] Flag to compute active tracer budgets + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer :: & + i, iCell, iEdge, &! horz indices + cell1, cell2, &! neighbor cell indices + k,kmin, kmax, &! vert index variants + kmin1, kmax1, &! vert index variants + iTracer ! tracer index + + real (kind=RKIND) :: & + signedFactor, &! temp factor including flux sign + tracerMinNew, &! updated tracer minimum + tracerMaxNew, &! updated tracer maximum + tracerUpwindNew, &! tracer updated with upwind flx + scaleFactor, &! factor for normalizing fluxes + flux, &! flux temporary + tracerWeight, &! tracer weighting temporary + invAreaCell1, &! inverse cell area + coef1, coef3 ! temporary coefficients + + real (kind=RKIND), dimension(:), allocatable :: & + wgtTmp, &! vertical temporaries for + flxTmp, sgnTmp ! high-order flux computation + + real (kind=RKIND), dimension(:,:), allocatable :: & + tracerCur, &! reordered current tracer + tracerMax, &! max tracer in neighbors for limiting + tracerMin, &! min tracer in neighbors for limiting + hNewInv, &! inverse of new layer thickness + hProv, &! provisional layer thickness + hProvInv, &! inverse of provisional layer thickness + flxIn, &! flux coming into each cell + flxOut, &! flux going out of each cell + workTend, &! temp for holding some tendency values + lowOrderFlx, &! low order flux for FCT + highOrderFlx ! high order flux for FCT + + real (kind=RKIND), parameter :: & + eps = 1.e-10_RKIND ! small number to avoid numerical difficulties + + ! end of preamble + !---------------- + ! begin code + +#ifdef _ADV_TIMERS + call mpas_timer_start('allocates') +#endif + + ! allocate temporary arrays + allocate(wgtTmp (nVertLevels), & + flxTmp (nVertLevels), & + sgnTmp (nVertLevels), & + tracerCur (nVertLevels ,nCells+1), & + tracerMin (nVertLevels ,nCells), & + tracerMax (nVertLevels ,nCells), & + hNewInv (nVertLevels ,nCells), & + hProv (nVertLevels ,nCells), & + hProvInv (nVertLevels ,nCells), & + flxIn (nVertLevels ,nCells+1), & + flxOut (nVertLevels ,nCells+1), & + workTend (nVertLevels ,nCells+1), & + lowOrderFlx (nVertLevels+1,max(nCells,nEdges)+1), & + highOrderFlx(nVertLevels+1,max(nCells,nEdges)+1)) + + ! Initialize variables so you don't get floating point exceptions + wgtTmp(:) = 0.0_RKIND + flxTmp(:) = 0.0_RKIND + sgnTmp(:) = 0.0_RKIND + tracerCur(:,:) = 0.0_RKIND + tracerMin(:,:) = 0.0_RKIND + tracerMax(:,:) = 0.0_RKIND + hNewInv(:,:) = 0.0_RKIND + hProv(:,:) = 0.0_RKIND + hProvInv(:,:) = 0.0_RKIND + flxIn(:,:) = 0.0_RKIND + flxOut(:,:) = 0.0_RKIND + workTend(:,:) = 0.0_RKIND + lowOrderFlx(:,:) = 0.0_RKIND + highOrderFlx(:,:) = 0.0_RKIND + !$acc enter data create(wgtTmp, flxTmp, sgnTmp, & + !$acc tracerCur, tracerMin, tracerMax, & + !$acc hNewInv, hProv, hProvInv, flxIn, flxOut, & + !$acc workTend, lowOrderFlx, highOrderFlx) + +#ifdef _ADV_TIMERS + call mpas_timer_stop('allocates') + call mpas_timer_start('prov thickness') +#endif + + ! Compute some provisional layer thicknesses + ! Note: This assumes we are in the first part of the horizontal/ + ! vertical operator splitting, which is true because currently + ! we dont flip order and horizontal is always first. + ! See notes in commit 2cd4a89d. + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(areaCell, dvEdge, minLevelCell, maxLevelCell, & + !$acc nEdgesOnCell, edgesOnCell, edgeSignOnCell, & + !$acc layerThickness, normalThicknessFlux, w, & + !$acc hProv, hProvInv, hNewInv) & + !$acc private(k,kmin,kmax,i,iEdge, invAreaCell1, signedFactor) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k,kmin,kmax,i,iEdge, invAreaCell1, signedFactor) +#endif + do iCell = 1, nCells + invAreaCell1 = dt/areaCell(iCell) + kmin = 1 !minLevelCell(iCell) + kmax = nVertLevels !maxLevelCell(iCell) + do k = kmin, kmax + hProv(k, iCell) = layerThickness(k, iCell) + end do + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i,iCell) + signedFactor = invAreaCell1*dvEdge(iEdge)* & + edgeSignOnCell(i,iCell) + ! Provisional layer thickness is after horizontal + ! thickness flux only + do k = kmin, kmax + hProv(k,iCell) = hProv(k,iCell) & + + signedFactor*normalThicknessFlux(k,iEdge) + end do + end do + ! New layer thickness is after horizontal and vertical + ! thickness flux + do k = kmin, kmax + if (hProv(k, iCell) > 0.0_RKIND) then + hProvInv(k,iCell) = 1.0_RKIND/ hProv(k,iCell) + hNewInv (k,iCell) = 1.0_RKIND/(hProv(k,iCell) - & + dt*w(k,iCell) + dt*w(k+1, iCell)) + else + hProvInv(k, iCell) = 0.0_RKIND + hNewInv(k, iCell) = 0.0_RKIND + endif + end do + end do +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef _ADV_TIMERS + call mpas_timer_stop('prov thickness') +#endif + + ! Loop over tracers. One tracer is advected at a time. + do iTracer = 1, nTracers + +#ifdef _ADV_TIMERS + call mpas_timer_start('cell init') +#endif + + ! Extract current tracer and change index order to improve locality +#ifdef MPAS_OPENACC + !$acc parallel loop collapse(2) & + !$acc present(tracerCur, tracers) +#else + !$omp parallel + !$omp do schedule(runtime) private(k) +#endif + do iCell = 1, nCells+1 + do k=1, nVertLevels + tracerCur(k,iCell) = tracers(iTracer,k,iCell) + end do ! k loop + end do ! iCell loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + ! Compute the high and low order horizontal fluxes. +#ifdef _ADV_TIMERS + call mpas_timer_stop('cell init') + call mpas_timer_start('tracer bounds') +#endif + + ! set nCells to first halo level + ! nCells = nCellsHalo(1) + + ! Determine bounds on tracer (tracerMin and tracerMax) from + ! surrounding cells for later limiting. + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(minLevelCell, maxLevelCell, nEdgesOnCell, & + !$acc cellsOnCell, tracerCur, tracerMin, tracerMax) & + !$acc private(k, kmin, kmax, kmin1, kmax1, i, cell2) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(k, kmin, kmax, kmin1, kmax1, i, cell2) +#endif + do iCell = 1, nCellsHalo(1) + kmin = 1 !minLevelCell(iCell) + kmax = nVertLevels !maxLevelCell(iCell) + do k=kmin,kmax + tracerMin(k,iCell) = tracerCur(k,iCell) + tracerMax(k,iCell) = tracerCur(k,iCell) + end do + ! TODO: Determine how/if this translates to MALI + do i = 1, nEdgesOnCell(iCell) + cell2 = cellsOnCell(i,iCell) + !kmin1 = max(kmin, minLevelCell(cell2)) + !kmax1 = min(kmax, maxLevelCell(cell2)) + do k=kmin, kmax !kmin1,kmax1 + tracerMax(k,iCell) = max(tracerMax(k,iCell), & + tracerCur(k,cell2)) + tracerMin(k,iCell) = min(tracerMin(k,iCell), & + tracerCur(k,cell2)) + end do ! k loop + end do ! i loop over nEdgesOnCell + end do ! cell loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef _ADV_TIMERS + call mpas_timer_stop('tracer bounds') + call mpas_timer_start('horiz flux') +#endif + ! Need all the edges around the 1 halo cells and owned cells + ! nEdges = nEdgesHalo(2) + + ! Compute the high order horizontal flux + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(minLevelCell, maxLevelCell, minLevelEdgeBot, & + !$acc maxLevelEdgeTop, nAdvCellsForEdge, & + !$acc advCellsForEdge, cellsOnEdge, dvEdge, & + !$acc advMaskHighOrder, advCoefs, advCoefs3rd, & + !$acc tracerCur, normalThicknessFlux, & + !$acc highOrderFlx, lowOrderFlx) & + !$acc private(i, k, icell, cell1, cell2, coef1, coef3, & + !$acc wgtTmp, sgnTmp, flxTmp, tracerWeight) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, k, icell, cell1, cell2, coef1, coef3, & + !$omp wgtTmp, sgnTmp, flxTmp, tracerWeight) +#endif + do iEdge = 1, nEdgesHalo(2) + cell1 = cellsOnEdge(1, iEdge) + cell2 = cellsOnEdge(2, iEdge) + + ! compute some common intermediate factors + do k = 1, nVertLevels + wgtTmp(k) = normalThicknessFlux(k,iEdge)* & + advMaskHighOrder(k,iEdge) + sgnTmp(k) = sign(1.0_RKIND, & + normalThicknessFlux(k,iEdge)) + flxTmp(k) = 0.0_RKIND + end do + + ! Compute 3rd or 4th fluxes where requested. + do i = 1, nAdvCellsForEdge(iEdge) + iCell = advCellsForEdge(i,iEdge) + coef1 = advCoefs (i,iEdge) + coef3 = advCoefs3rd (i,iEdge)*coef3rdOrder + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + flxTmp(k) = flxTmp(k) + tracerCur(k,iCell)* & + wgtTmp(k)*(coef1 + coef3*sgnTmp(k)) + end do ! k loop + end do ! i loop over nAdvCellsForEdge + + do k=1,nVertLevels + highOrderFlx(k,iEdge) = flxTmp(k) + end do + + ! Compute 2nd order fluxes where needed. + ! Also compute low order upwind horizontal flux (monotonic) + ! Remove low order flux from the high order flux + ! Store left over high order flux in highOrderFlx array + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + tracerWeight = advMask2ndOrder(k,iEdge) & + * (1.0_RKIND - advMaskHighOrder(k, iEdge)) & + * (dvEdge(iEdge) * 0.5_RKIND) & + * normalThicknessFlux(k, iEdge) + + lowOrderFlx(k,iEdge) = dvEdge(iEdge) * & + (max(0.0_RKIND,normalThicknessFlux(k,iEdge))*tracerCur(k,cell1) & + + min(0.0_RKIND,normalThicknessFlux(k,iEdge))*tracerCur(k,cell2)) + + highOrderFlx(k,iEdge) = highOrderFlx(k,iedge) & + + tracerWeight*(tracerCur(k,cell1) & + + tracerCur(k,cell2)) + + ! only remove low order flux where high order flux is valid to + ! avoid introducing nonzero values to highOrderFlx where it is invalid + highOrderFlx(k,iEdge) = highOrderFlx(k,iEdge) & + - lowOrderFlx(k,iEdge) * & + (advMaskHighOrder(k,iEdge) + advMask2ndOrder(k,iEdge) & + * (1.0_RKIND - advMaskHighOrder(k, iEdge)) ) + end do ! k loop + end do ! iEdge loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef _ADV_TIMERS + call mpas_timer_stop('horiz flux') + call mpas_timer_start('scale factor build') +#endif + + ! Initialize flux arrays for all cells + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present (workTend, flxIn, flxOut) +#else + !$omp parallel + !$omp do schedule(runtime) +#endif + do iCell = 1, nCells+1 + do k=1, nVertLevels + workTend(k, iCell) = 0.0_RKIND + flxIn (k, iCell) = 0.0_RKIND + flxOut (k, iCell) = 0.0_RKIND + end do ! k loop + end do ! iCell loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + + ! Need one halo of cells around owned cells + ! nCells = nCellsHalo(1) + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(minLevelEdgeBot, maxLevelEdgeTop, & + !$acc minLevelCell, maxLevelCell, areaCell, & + !$acc nEdgesOnCell, edgesOnCell, cellsOnEdge, & + !$acc edgeSignOnCell, layerThickness, hProvInv, & + !$acc tracerCur, tracerMax, tracerMin, workTend, & + !$acc flxOut, flxIn, lowOrderFlx, highOrderFlx) & + !$acc private(i, k, iEdge, cell1, cell2, & + !$acc invAreaCell1, signedFactor, scaleFactor, & + !$acc tracerUpwindNew, tracerMinNew, tracerMaxNew) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, k, iEdge, cell1, cell2, & + !$omp invAreaCell1, signedFactor, scaleFactor, & + !$omp tracerUpwindNew, tracerMinNew, tracerMaxNew) +#endif + do iCell = 1, nCellsHalo(1) + invAreaCell1 = 1.0_RKIND / areaCell(iCell) + + ! Finish computing the low order horizontal fluxes + ! Upwind fluxes are accumulated in workTend + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + signedFactor = edgeSignOnCell(i, iCell) * invAreaCell1 + + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + + ! Here workTend is the advection tendency due to the + ! upwind (low order) fluxes. + workTend(k,iCell) = workTend(k,iCell) & + + signedFactor*lowOrderFlx(k,iEdge) + + ! Accumulate remaining high order fluxes + flxOut(k,iCell) = flxOut(k,iCell) + min(0.0_RKIND, & + signedFactor*highOrderFlx(k,iEdge)) + flxIn (k,iCell) = flxIn (k,iCell) + max(0.0_RKIND, & + signedFactor*highOrderFlx(k,iEdge)) + + end do + end do + + ! Build the factors for the FCT + ! Computed using the bounds that were computed previously, + ! and the bounds on the newly updated value + ! Factors are placed in the flxIn and flxOut arrays + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + ! Here workTend is the upwind tendency + tracerUpwindNew = (tracerCur(k,iCell)*layerThickness(k,iCell) & + + dt*workTend(k,iCell))*hProvInv(k,iCell) + tracerMinNew = tracerUpwindNew & + + dt*flxOut(k,iCell)*hProvInv(k,iCell) + tracerMaxNew = tracerUpwindNew & + + dt*flxIn (k,iCell)*hProvInv(k,iCell) + + scaleFactor = (tracerMax(k,iCell) - tracerUpwindNew)/ & + (tracerMaxNew - tracerUpwindNew + eps) + flxIn (k,iCell) = min(1.0_RKIND, & + max(0.0_RKIND, scaleFactor)) + scaleFactor = (tracerUpwindNew - tracerMin(k,iCell))/ & + (tracerUpwindNew - tracerMinNew + eps) + flxOut(k,iCell) = min(1.0_RKIND, & + max(0.0_RKIND, scaleFactor)) + end do ! k loop + end do ! iCell loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef _ADV_TIMERS + call mpas_timer_stop('scale factor build') + call mpas_timer_start('rescale horiz fluxes') +#endif + ! Need all of the edges around owned cells + ! nEdges = nEdgesHalo(1) + ! rescale the high order horizontal fluxes + +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(minLevelEdgeBot, maxLevelEdgeTop, & + !$acc cellsOnEdge, highOrderFlx, flxIn, flxOut) & + !$acc private(k, cell1, cell2) +#else + !$omp parallel + !$omp do schedule(runtime) private(k, cell1, cell2) +#endif + do iEdge = 1, nEdgesHalo(1) + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + highOrderFlx(k,iEdge) = max(0.0_RKIND,highOrderFlx(k,iEdge))* & + min(flxOut(k,cell1), flxIn (k,cell2)) & + + min(0.0_RKIND,highOrderFlx(k,iEdge))* & + min(flxIn (k,cell1), flxOut(k,cell2)) + end do ! k loop + end do ! iEdge loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef _ADV_TIMERS + call mpas_timer_stop('rescale horiz fluxes') + call mpas_timer_start('flux accumulate') +#endif + + ! Accumulate the scaled high order vertical tendencies + ! and the upwind tendencies +#ifdef MPAS_OPENACC + !$acc parallel loop & + !$acc present(minLevelEdgeBot, maxLevelEdgeTop, & + !$acc minLevelCell, maxLevelCell, areaCell, & + !$acc nEdgesOnCell, edgesOnCell, edgeSignOnCell, & + !$acc workTend, highOrderFlx, layerThickness, & + !$acc tracerCur, hProvInv, tend) & + !$acc private(i, k, iEdge, invAreaCell1, signedFactor) +#else + !$omp parallel + !$omp do schedule(runtime) & + !$omp private(i, k, iEdge, invAreaCell1, signedFactor) +#endif + do iCell = 1, nCellsSolve + invAreaCell1 = 1.0_RKIND / areaCell(iCell) + + ! Accumulate the scaled high order horizontal tendencies + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + signedFactor = invAreaCell1*edgeSignOnCell(i,iCell) + do k = 1, nVertLevels !minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge) + ! workTend on RHS is upwind tendency + ! workTend on LHS is total horiz advect tendency + workTend(k,iCell) = workTend(k,iCell) & + + signedFactor*highOrderFlx(k,iEdge) + end do + end do + + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + ! workTend on RHS is total horiz advection tendency + ! tracerCur on LHS is provisional tracer after + ! horizontal fluxes only. + tracerCur(k,iCell) = (tracerCur(k,iCell)* & + layerThickness(k,iCell) & + + dt*workTend(k,iCell)) & + * hProvInv(k,iCell) + tend(iTracer,k,iCell) = tend(iTracer,k,iCell) & + + workTend(k,iCell) + end do + + end do ! iCell loop +#ifndef MPAS_OPENACC + !$omp end do + !$omp end parallel +#endif + +#ifdef _ADV_TIMERS + call mpas_timer_stop('flux accumulate') + call mpas_timer_start('advect diags horiz') +#endif + ! Compute budget and monotonicity diagnostics if needed +! if (computeBudgets) then + + !nEdges = nEdgesHalo(2) +#ifdef MPAS_OPENACC + !$acc parallel loop private(k) & + !$acc present(activeTracerHorizontalAdvectionEdgeFlux, & + !$acc minLevelEdgeBot, maxLevelEdgeTop, & + !$acc lowOrderFlx, highOrderFlx, dvEdge) +#endif + +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(k) +!#endif + + ! Use activeTracerHorizontalAdvectionEdgeFlux from the call to fct for thickness + ! advection as the higher-order normalThicknessFlux for fct tracer advection. + if (present(activeTracerHorizontalAdvectionEdgeFlux)) then + do iEdge = 1,nEdges + do k = 1, nVertLevels + ! Save u*h*T flux on edge for analysis. This variable will be + ! divided by h at the end of the time step. + if (dvEdge(iEdge) > 0.0_RKIND) then + activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = & + (lowOrderFlx(k,iEdge) + highOrderFlx(k,iEdge))/dvEdge(iEdge) + else + activeTracerHorizontalAdvectionEdgeFlux(iTracer,k,iEdge) = 0.0_RKIND + endif + enddo + enddo + endif +!#ifndef MPAS_OPENACC +! !$omp end do +!#endif + +#ifdef MPAS_OPENACC + !$acc parallel loop private(k) & + !$acc present(activeTracerHorizontalAdvectionTendency, & + !$acc minLevelCell, maxLevelCell, workTend) + +#endif + +!#else +! !$omp do schedule(runtime) private(k) +!#endif +! ! TODO: Determine if it is necessary to define activeTracerHorizontalAdvectionTendency +! ! for MALI +! !do iCell = 1, nCellsSolve +! !do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) +! ! activeTracerHorizontalAdvectionTendency(iTracer,k,iCell) = & +! ! workTend(k,iCell) +! !end do +! !end do ! iCell loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif + +! end if ! computeBudgets + + if (monotonicityCheck) then + ! Check tracer values against local min,max to detect + ! non-monotone values and write warning if found + + ! Perform check on host since print involved + !$acc update host(tracerCur, tracerMin, tracerMax) + + !$omp parallel + !$omp do schedule(runtime) private(k) + do iCell = 1, nCellsSolve + do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) + if(tracerCur(k,iCell) < tracerMin(k, iCell)-eps) then + call mpas_log_write( & + 'Horizontal minimum out of bounds on tracer: $i $r $r ', & + MPAS_LOG_WARN, intArgs=(/iTracer/), & + realArgs=(/ tracerMin(k, iCell), tracerCur(k,iCell) /) ) + end if + + if(tracerCur(k,iCell) > tracerMax(k,iCell)+eps) then + call mpas_log_write( & + 'Horizontal maximum out of bounds on tracer: $i $r $r ', & + MPAS_LOG_WARN, intArgs=(/iTracer/), & + realArgs=(/ tracerMax(k, iCell), tracerCur(k,iCell) /) ) + end if + end do + end do + !$omp end do + !$omp end parallel + end if ! monotonicity check +#ifdef _ADV_TIMERS + call mpas_timer_stop('advect diags horiz') +#endif + +!TODO: implement vertical advection? + !----------------------------------------------------------------------- + ! + ! Horizontal advection complete + ! Begin vertical advection + ! + !----------------------------------------------------------------------- + +!#ifdef _ADV_TIMERS +! call mpas_timer_start('cell init vert') +!#endif +! ! Initialize variables for use in this iTracer iteration +! +!#ifdef MPAS_OPENACC +! !$acc parallel loop collapse(2) present(workTend) +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(k) +!#endif +! do iCell = 1, nCells +! do k=1, nVertLevels +! workTend(k, iCell) = 0.0_RKIND +! end do ! k loop +! end do ! iCell loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif +! +!#ifdef _ADV_TIMERS +! call mpas_timer_stop('cell init vert') +! call mpas_timer_start('vertical bounds') +!#endif +! +! ! Need all owned and 1 halo cells +! ! nCells = nCellsHalo(1) +! +! ! Determine bounds on tracerCur from neighbor values for limiting +!#ifdef MPAS_OPENACC +! !$acc parallel loop & +! !$acc present(tracerCur, tracerMin, tracerMax, & +! !$acc minLevelCell, maxLevelCell) & +! !$acc private(k, kmin, kmax) +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(k, kmin, kmax) +!#endif +! do iCell = 1, nCellsHalo(1) +! kmin = 1 !minLevelCell(iCell) +! kmax = nVertLevels !maxLevelCell(iCell) +! +! ! take care of top cell +! tracerMax(kmin,iCell) = max(tracerCur(1,iCell), & +! tracerCur(2,iCell)) +! tracerMin(kmin,iCell) = min(tracerCur(1,iCell), & +! tracerCur(2,iCell)) +! do k=kmin+1,kmax-1 +! tracerMax(k,iCell) = max(tracerCur(k-1,iCell), & +! tracerCur(k ,iCell), & +! tracerCur(k+1,iCell)) +! tracerMin(k,iCell) = min(tracerCur(k-1,iCell), & +! tracerCur(k ,iCell), & +! tracerCur(k+1,iCell)) +! end do +! ! finish with bottom cell +! tracerMax(kmax,iCell) = max(tracerCur(kmax ,iCell), & +! tracerCur(kmax-1,iCell)) +! tracerMin(kmax,iCell) = min(tracerCur(kmax ,iCell), & +! tracerCur(kmax-1,iCell)) +! end do ! cell loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif +! +!!#ifdef _ADV_TIMERS +!! call mpas_timer_stop('vertical bounds') +!! call mpas_timer_start('vertical flux high') +!!#endif +!! +!! ! Compute the high order vertical fluxes based on order selected. +!! +!! call tracer_advection_vert_flx(tracerCur, w, hProv, & +!! highOrderFlx) +!! +!!#ifdef _ADV_TIMERS +!! call mpas_timer_stop('vertical flux high') +!! call mpas_timer_start('vertical flux') +!!#endif +! +! ! Compute low order upwind vertical flux (monotonic) +! ! Remove low order flux from the high order flux. +! ! Store left over high order flux in highOrderFlx array. +! +!#ifdef MPAS_OPENACC +! !$acc parallel loop & +! !$acc present(minLevelCell, maxLevelCell, w, tracerCur, & +! !$acc lowOrderFlx, highOrderFlx, workTend, & +! !$acc flxIn, flxOut) & +! !$acc private(k, kmin, kmax) +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(k, kmin, kmax) +!#endif +! do iCell = 1, nCellsHalo(1) +! kmin = 1 !minLevelCell(iCell) +! kmax = nVertLevels !maxLevelCell(iCell) +! +! lowOrderFlx(kmin,iCell) = 0.0_RKIND +! do k = kmin+1, kmax +! lowOrderFlx(k,iCell) = & +! min(0.0_RKIND,w(k,iCell))*tracerCur(k-1,iCell) + & +! max(0.0_RKIND,w(k,iCell))*tracerCur(k ,iCell) +! highOrderFlx(k,iCell) = highOrderFlx(k,iCell) & +! - lowOrderFlx(k,iCell) +! end do ! k loop +! lowOrderFlx(kmax+1,iCell) = 0.0_RKIND +! +! ! Upwind fluxes are accumulated in workTend +! ! flxIn contains total remaining high order flux into iCell +! ! it is positive. +! ! flxOut contains total remaining high order flux out of iCell +! ! it is negative +! do k = kmin,kmax +! workTend(k,iCell) = lowOrderFlx(k+1,iCell) & +! - lowOrderFlx(k ,iCell) +! flxIn (k,iCell) = max(0.0_RKIND, highOrderFlx(k+1,iCell)) & +! - min(0.0_RKIND, highOrderFlx(k ,iCell)) +! flxOut(k,iCell) = min(0.0_RKIND, highOrderFlx(k+1,iCell)) & +! - max(0.0_RKIND, highOrderFlx(k ,iCell)) +! end do ! k Loop +! end do ! iCell Loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif +! +!#ifdef _ADV_TIMERS +! call mpas_timer_stop('vertical flux') +! call mpas_timer_start('scale factor build') +!#endif +! +! ! Build the scale factors to limit flux for FCT +! ! Computed using the bounds that were computed previously, +! ! and the bounds on the newly updated value +! ! Factors are placed in the flxIn and flxOut arrays +! +! !nCells = nCellsHalo(1) ! Need one halo around owned cells +! +!#ifdef MPAS_OPENACC +! !$acc parallel loop & +! !$acc present(minLevelCell, maxLevelCell, hProv, hNewInv, & +! !$acc tracerCur, tracerMax, tracerMin, & +! !$acc workTend, flxIn, flxOut) & +! !$acc private(k, tracerMinNew, tracerMaxNew, & +! !$acc tracerUpwindNew, scaleFactor) +!#else +! !$omp parallel +! !$omp do schedule(runtime) & +! !$omp private(k, tracerMinNew, tracerMaxNew, & +! !$omp tracerUpwindNew, scaleFactor) +!#endif +! do iCell = 1, nCellsHalo(1) +! +! do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) +! ! workTend on the RHS is the upwind tendency +! tracerMinNew = (tracerCur(k,iCell)*hProv(k,iCell) & +! + dt*(workTend(k,iCell)+flxOut(k,iCell))) & +! * hNewInv(k,iCell) +! tracerMaxNew = (tracerCur(k,iCell)*hProv(k,iCell) & +! + dt*(workTend(k,iCell)+flxIn(k,iCell))) & +! * hNewInv(k,iCell) +! tracerUpwindNew = (tracerCur(k,iCell)*hProv(k,iCell) & +! + dt*workTend(k,iCell)) * hNewInv(k,iCell) +! +! scaleFactor = max(0.0_RKIND, (tracerMax(k,iCell)-tracerUpwindNew)/ & +! (tracerMaxNew-tracerUpwindNew+eps) ) +! flxIn (k,iCell) = min(1.0_RKIND, scaleFactor) +! +! +! scaleFactor = max(0.0_RKIND, (tracerUpwindNew-tracerMin(k,iCell))/ & +! (tracerUpwindNew-tracerMinNew+eps) ) +! flxOut(k,iCell) = min(1.0_RKIND, scaleFactor) +! end do ! k loop +! end do ! iCell loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif +! +!#ifdef _ADV_TIMERS +! call mpas_timer_stop('scale factor build') +! call mpas_timer_start('flux accumulate') +!#endif +! +! ! Accumulate the scaled high order vertical tendencies +! ! and the upwind tendencies +! +!#ifdef MPAS_OPENACC +! !$acc parallel loop & +! !$acc present(minLevelCell, maxLevelCell, highOrderFlx, & +! !$acc flxIn, flxOut, workTend, tend) & +! !$acc private(k, kmin, kmax, flux) +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(k, kmin, kmax, flux) +!#endif +! do iCell = 1, nCellsSolve +! kmin = 1 !minLevelCell(iCell) +! kmax = nVertLevels !maxLevelCell(iCell) +! ! rescale the high order vertical flux +! do k = kmin+1, kmax +! flux = highOrderFlx(k,iCell) +! highOrderFlx(k,iCell) = max(0.0_RKIND,flux)* & +! min(flxOut(k ,iCell), & +! flxIn (k-1,iCell)) & +! + min(0.0_RKIND,flux)* & +! min(flxOut(k-1,iCell), & +! flxIn (k ,iCell)) +! end do ! k loop +! +! do k = kmin,kmax +! ! workTend on the RHS is upwind tendency +! ! workTend on the LHS is total vertical advection tendency +! workTend(k, iCell) = workTend(k, iCell) & +! + (highOrderFlx(k+1,iCell) & +! - highOrderFlx(k ,iCell)) +! tend(iTracer,k,iCell) = tend(iTracer,k,iCell) & +! + workTend(k,iCell) +! end do ! k loop +! +! end do ! iCell loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif +! +! ! Compute advection diagnostics and monotonicity checks if requested +!#ifdef _ADV_TIMERS +! call mpas_timer_stop('flux accumulate') +! call mpas_timer_start('advect diags vert') +!#endif +! if (computeBudgets) then +! +!#ifdef MPAS_OPENACC +! !$acc parallel loop & +! !$acc present(lowOrderFlx, highOrderFlx, workTend, & +! !$acc activeTracerVerticalAdvectionTopFlux, & +! !$acc activeTracerVerticalAdvectionTendency, & +! !$acc minLevelCell, maxLevelCell) & +! !$acc private(k, kmin, kmax) +!#endif +! +!!#else +!! !$omp parallel +!! !$omp do schedule(runtime) private(k,kmin,kmax) +!!#endif +!! do iCell = 1, nCellsSolve +!! kmin = minLevelCell(iCell) +!! kmax = maxLevelCell(iCell) +!! do k = kmin+1,kmax +!! activeTracerVerticalAdvectionTopFlux(iTracer,k,iCell) = & +!! lowOrderFlx(k,iCell) + highOrderFlx(k,iCell) +!! end do +!! do k = kmin,kmax +!! activeTracerVerticalAdvectionTendency(iTracer,k,iCell) = & +!! workTend(k,iCell) +!! end do +!! end do ! iCell loop +!!#ifndef MPAS_OPENACC +!! !$omp end do +!! !$omp end parallel +!!#endif +! end if ! computeBudgets +! +! if (monotonicityCheck) then +! +! ! Check for monotonicity of new tracer value +! ! Use flxIn as a temp for new tracer value +! +!#ifdef MPAS_OPENACC +! !$acc parallel loop collapse(2) & +! !$acc present(flxIn, tracerCur, hProv, hNewInv, & +! !$acc workTend, cellMask) +!#else +! !$omp parallel +! !$omp do schedule(runtime) private(k) +!#endif +! do iCell = 1,nCellsSolve +! do k = 1,nVertLevels +! ! workTend on the RHS is total vertical advection tendency +! flxIn(k,iCell) = (tracerCur(k, iCell)*hProv(k, iCell) & +! + dt*workTend(k, iCell)) & +! * hNewInv(k,iCell)*cellMask(iCell) +! end do ! vert +! end do ! cell loop +!#ifndef MPAS_OPENACC +! !$omp end do +! !$omp end parallel +!#endif +! +! ! Transfer new tracer values to host for diagnostic check +! !$acc update host(flxIn) +! +! ! Print out info for cells that are out of bounds +! ! flxIn holds new tracer values +! do iCell = 1, nCellsSolve +! do k = 1, nVertLevels !minLevelCell(iCell), maxLevelCell(iCell) +! if (flxIn(k,iCell) < tracerMin(k, iCell)-eps) then +! call mpas_log_write( & +! 'Vertical minimum out of bounds on tracer: $i $i $i $r $r ',& +! MPAS_LOG_WARN, intArgs=(/iTracer, k, iCell/), & +! realArgs=(/ tracerMin(k,iCell), flxIn(k,iCell) /) ) +! end if +! if (flxIn(k,iCell) > tracerMax(k,iCell)+eps) then +! call mpas_log_write( & +! 'Vertical maximum out of bounds on tracer: $i $i $i $r $r ',& +! MPAS_LOG_WARN, intArgs=(/iTracer, k, iCell/), & +! realArgs=(/ tracerMax(k,iCell), flxIn(k,iCell) /) ) +! end if +! end do +! end do +! +! end if ! monotonicity check +!#ifdef _ADV_TIMERS +! call mpas_timer_stop('advect diags vert') +!#endif + ! Update tracer array + tracers(iTracer,:,:) = tracerCur(:,:) + end do ! iTracer loop +! +#ifdef _ADV_TIMERS + call mpas_timer_start('deallocates') +#endif + !$acc exit data delete(wgtTmp, flxTmp, sgnTmp, & + !$acc tracerCur, tracerMin, tracerMax, & + !$acc hNewInv, hProv, hProvInv, flxIn, flxOut, & + !$acc workTend, lowOrderFlx, highOrderFlx) + + deallocate(wgtTmp, & + flxTmp, & + sgnTmp, & + tracerCur, & + tracerMin, & + tracerMax, & + hNewInv, & + hProv, & + hProvInv, & + flxIn, & + flxOut, & + workTend, & + lowOrderFlx, & + highOrderFlx) + +#ifdef _ADV_TIMERS + call mpas_timer_stop('deallocates') +#endif + + end subroutine li_tracer_advection_fct_tend!}}} + +!********************************************************************** +! +! routine li_tracer_advection_fct_init +! +!> \brief MPAS initialize monotonic tracer advection tendency with FCT +!> \author Mark Petersen, David Lee, Doug Jacobsen, Phil Jones +!> \date October 2017, updated May 2019 +!> \details +!> This routine initializes monotonic tracer advection quantities for +!> the flux-corrected transport (FCT) algorithm. +! +!----------------------------------------------------------------------- + + subroutine li_tracer_advection_fct_init(err)!{{{ + + !*** output parameters + + integer, intent(out) :: & + err !< [out] error flag + + ! end of preamble + !---------------- + ! begin code + + err = 0 ! initialize error code to success + + ! Check that the halo is wide enough for FCT + if (config_num_halos < 3) then + call mpas_log_write( & + 'Monotonic advection cannot be used with less than 3 halos.', & + MPAS_LOG_CRIT) + err = -1 + end if + + ! Set blending coefficient if 3rd order horizontal advection chosen + select case (config_horiz_tracer_adv_order) + case (2) + coef3rdOrder = 0.0_RKIND + case (3) + coef3rdOrder = config_advection_coef_3rd_order + case (4) + coef3rdOrder = 0.0_RKIND + case default + coef3rdOrder = 0.0_RKIND + call mpas_log_write( & + 'Invalid value for horizontal advection order, defaulting to 2',& + MPAS_LOG_WARN) + end select ! horizontal advection order + + ! Set flag for checking monotonicity + monotonicityCheck = config_check_tracer_monotonicity + + end subroutine li_tracer_advection_fct_init!}}} + +!*********************************************************************** + +end module li_tracer_advection_fct + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct_shared.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct_shared.F new file mode 100644 index 000000000000..ce6748e9ea37 --- /dev/null +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_advection_fct_shared.F @@ -0,0 +1,781 @@ +! 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.com/license.html +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! mpas_tracer_advection_shared +! +!> \brief MPAS ocean tracer advection common coefficients and variables +!> \author Doug Jacobsen +!> \date 03/09/12 +!> \details +!> This module contains the routines and arrays for some common +!> tracer advection quantities. +! +!----------------------------------------------------------------------- + +module li_tracer_advection_fct_shared + + use mpas_kind_types + use mpas_derived_types + use mpas_hash + use mpas_sort + use mpas_geometry_utils + use mpas_log + + use li_config + use li_mesh + use li_mask + + implicit none + private + save + + !-------------------------------------------------------------------- + ! Public variables (should only be used by other advection modules) + !-------------------------------------------------------------------- + + integer, public :: & + nAdvCellsMax ! largest number of advection cells for any edge + + integer, public, dimension(:), allocatable :: & + nAdvCellsForEdge ! number of cells contrib to advection at edge + + integer, public, dimension(:,:), allocatable :: & + advCellsForEdge ! index of cells contributing to advection at edge + + real (kind=RKIND), public, dimension(:,:), allocatable :: & + advCoefs, &! common advection coefficients + advCoefs3rd, &! common advection coeffs for high order + advMaskHighOrder, & ! mask where high order advection terms are valid + advMask2ndOrder ! mask where 2nd order advection terms are valid + + !-------------------------------------------------------------------- + ! Public member functions + !-------------------------------------------------------------------- + + public :: li_tracer_advection_fct_shared_init + +!*********************************************************************** + + contains + +!*********************************************************************** +! +! routine li_tracer_advection_fct_shared_init +! +!> \brief MPAS tracer advection coefficients +!> \author Doug Jacobsen, Bill Skamarock, +!> adapted for MALI by Trevor Hillebrand and Matt Hoffman +!> \date 03/09/12, adapted for MALI 03/2023 +!> \details +!> This routine precomputes advection coefficients for horizontal +!> advection of tracers. +! +!----------------------------------------------------------------------- + + subroutine li_tracer_advection_fct_shared_init(geometryPool, err)!{{{ + use li_config + use li_mesh + !----------------------------------------------------------------- + ! Output variables + !----------------------------------------------------------------- + + integer, intent(out) :: err !< [out] Error flag + + !----------------------------------------------------------------- + ! Input variables + !----------------------------------------------------------------- + + type (mpas_pool_type), intent(in) :: geometryPool + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer :: & + n, i, k, &! loop indices for neighbor, vertical loops + iEdge, iCell, &! loop indices for edge, cell loops + cell1, cell2, & ! neighbor cell indices across edge + iNeighbor, jCell + + integer, dimension(:), allocatable :: & + cellIndx ! cell indices gathered from neighbors + + integer, dimension(:), pointer :: cellMask + + integer, dimension(:,:), allocatable :: & + cellIndxSorted ! nbr cell indices sorted + + real (kind=RKIND), dimension(:,:,:), allocatable :: & + derivTwo !< 2nd derivative values for polynomial fit to tracers + + type (hashtable) :: cell_hash + + ! End preamble + !------------- + ! Begin code + + call mpas_pool_get_array(geometryPool, 'cellMask', cellMask) + + err = 0 ! initialize error code + + ! define max number of advection cells as max number of edges*2 + nAdvCellsMax = maxEdges2 + + ! Allocate common variables + allocate( nAdvCellsForEdge ( nEdges), & + advCellsForEdge (nAdvCellsMax,nEdges), & + advCoefs (nAdvCellsMax,nEdges), & + advCoefs3rd (nAdvCellsMax,nEdges), & + advMaskHighOrder(nVertLevels, nEdges), & + advMask2ndOrder(nVertLevels, nEdges), & + derivTwo (nAdvCellsMax,2,nEdges)) + + ! Initialize all these allocatable arrays + nAdvCellsForEdge(:) = 0 + advCellsForEdge(:,:) = 0 + advCoefs(:,:) = 0.0_RKIND + advCoefs3rd(:,:) = 0.0_RKIND + advMaskHighOrder(:,:) = 0.0_RKIND + advMask2ndOrder(:,:) = 0.0_RKIND + derivTwo(:,:,:) = 0.0_RKIND + + ! Compute derivTwo array + call computeDerivTwo(derivTwo, err) + if (err /= 0) then + call mpas_log_write( & + 'Error computing derivTwo in ocn advect shared init ', & + MPAS_LOG_CRIT) + deallocate (derivTwo) + return + endif + + allocate(cellIndx ( maxEdges2 + 2), & + cellIndxSorted(2, maxEdges2 + 2)) + + cellIndx(:) = 0 + cellIndxSorted(:,:) = 0 + + ! calculate boundaryCell. A boundary cell is one that + ! has at least one empty or non-dynamic neighbor + boundaryCell(:) = 0 + do iCell = 1, nCells + do iNeighbor = 1, nEdgesOnCell(iCell) + jCell = cellsOnCell(iNeighbor, iCell) + if (.not. li_mask_is_dynamic_ice(cellMask(jCell))) then + boundaryCell(iCell) = 1 + exit ! no need to loop over more neigbors + endif + enddo + enddo + + call mpas_log_write('Marked $i boundary cells out of $i cells belonging to this processor', & + intArgs=(/sum(boundaryCell(1:nCellsSolve)), nCellsSolve/)) + + ! Compute masks and coefficients + do iEdge = 1, nEdges + nAdvCellsForEdge(iEdge) = 0 + cell1 = cellsOnEdge(1,iEdge) + cell2 = cellsOnEdge(2,iEdge) + + ! at boundaries, must stay at low order + if (boundaryCell(cell1) == 1 .or. & + boundaryCell(cell2) == 1) then + advMaskHighOrder(:,iEdge) = 0.0_RKIND + else + advMaskHighOrder(:,iEdge) = 1.0_RKIND + end if + + ! second order mask must have ice on both sides + if ( li_mask_is_dynamic_ice(cellMask(cell1)) .and. li_mask_is_dynamic_ice(cellMask(cell2)) ) then + advMask2ndOrder(:,iEdge) = 1.0_RKIND + else + advMask2ndOrder(:,iEdge) = 0.0_RKIND + end if + + ! do only if this edge flux is needed to update owned cells + if (cell1 <= nCells .and. cell2 <= nCells) then + ! Insert cellsOnEdge to list of advection cells + call mpas_hash_init(cell_hash) + call mpas_hash_insert(cell_hash, cell1) + call mpas_hash_insert(cell_hash, cell2) + cellIndx(1) = cell1 + cellIndx(2) = cell2 + cellIndxSorted(1,1) = indexToCellID(cell1) + cellIndxSorted(2,1) = cell1 + cellIndxSorted(1,2) = indexToCellID(cell2) + cellIndxSorted(2,2) = cell2 + n = 2 + + ! Build unique list of cells used for advection on edge + ! by expanding to the extended neighbor cells + do i = 1, nEdgesOnCell(cell1) + if (.not. mpas_hash_search(cell_hash, & + cellsOnCell(i,cell1))) then + n = n + 1 + cellIndx(n) = cellsOnCell(i, cell1) + cellIndxSorted(1,n) = indexToCellID(cellsOnCell(i,cell1)) + cellIndxSorted(2,n) = cellsOnCell(i,cell1) + call mpas_hash_insert(cell_hash,cellsOnCell(i,cell1)) + end if + end do + + do i = 1, nEdgesOnCell(cell2) + if(.not. mpas_hash_search(cell_hash, & + cellsOnCell(i, cell2))) then + n = n + 1 + cellIndx(n) = cellsOnCell(i, cell2) + cellIndxSorted(1,n) = indexToCellID(cellsOnCell(i,cell2)) + cellIndxSorted(2,n) = cellsOnCell(i,cell2) + call mpas_hash_insert(cell_hash, cellsOnCell(i,cell2)) + end if + end do + + call mpas_hash_destroy(cell_hash) + + ! sort the cell indices by cellID + call mpas_quicksort(n, cellIndxSorted) + + ! store local cell indices for high-order calculations + nAdvCellsForEdge(iEdge) = n + do iCell = 1, nAdvCellsForEdge(iEdge) + advCellsForEdge(iCell,iEdge) = cellIndxSorted(2,iCell) + end do + + ! equation 7 in Skamarock, W. C., & Gassmann, A. (2011): + ! F(u,psi)_{i+1/2} = u_{i+1/2} * + ! [ 1/2 (psi_{i+1} + psi_i) term 1 + ! - 1/12(dx^2psi_{i+1} + dx^2psi_i) term 2 + ! + sign(u) beta/12 (dx^2psi_{i+1} - dx^2psi_i)] term 3 + ! (note minus sign) + ! + ! advCoefs accounts for terms 1 and 2 in SG11 equation 7. + ! Term 1 is the 2nd-order flux-function term. advCoefs + ! accounts for this with the "+ 0.5" lines below. In the + ! advection routines that use these coefficients, the + ! 2nd-order flux loop is then skipped. Term 2 is the + ! 4th-order flux-function term. advCoefs accounts for + ! term 3, the beta term. beta > 0 corresponds to the + ! third-order flux function. The - sign in the derivTwo + ! accumulation is for the i+1 part of term 3, while + ! the + sign is for the i part. + + do i=1,nAdvCellsMax + advCoefs (i,iEdge) = 0.0_RKIND + advCoefs3rd(i,iEdge) = 0.0_RKIND + end do + + ! pull together third and fourth order contributions to the + ! flux first from cell1 + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), indexToCellID(cell1)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(1,1,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + + derivTwo(1,1,iEdge) + end if + + do iCell = 1, nEdgesOnCell(cell1) + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cellsOnCell(iCell,cell1))) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(iCell+1,1,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + + derivTwo(iCell+1,1,iEdge) + end if + end do + + ! pull together third and fourth order contributions to the + ! flux now from cell2 + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), indexToCellID(cell2)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(1,2,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + - derivTwo(1,2,iEdge) + end if + + do iCell = 1, nEdgesOnCell(cell2) + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cellsOnCell(iCell,cell2))) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs (i,iEdge) = advCoefs (i,iEdge) & + + derivTwo(iCell+1,2,iEdge) + advCoefs3rd(i,iEdge) = advCoefs3rd(i,iEdge) & + - derivTwo(iCell+1,2,iEdge) + end if + end do + + do iCell = 1,nAdvCellsForEdge(iEdge) + advCoefs (iCell,iEdge) = - (dcEdge(iEdge)**2)* & + advCoefs(iCell,iEdge) / 12. + advCoefs3rd(iCell,iEdge) = - (dcEdge(iEdge)**2)* & + advCoefs3rd(iCell,iEdge) / 12. + end do + + ! 2nd order centered contribution + ! place this in the main flux weights + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cell1)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs(i,iEdge) = advCoefs(i, iEdge) + 0.5 + end if + + i = mpas_binary_search(cellIndxSorted, 2, 1, & + nAdvCellsForEdge(iEdge), & + indexToCellID(cell2)) + if (i <= nAdvCellsForEdge(iEdge)) then + advCoefs(i,iEdge) = advCoefs(i, iEdge) + 0.5 + end if + + ! multiply by edge length - thus the flux is just dt*ru + ! times the results of the vector-vector multiply + do iCell=1,nAdvCellsForEdge(iEdge) + advCoefs (iCell,iEdge) = dvEdge(iEdge)* & + advCoefs (iCell,iEdge) + advCoefs3rd(iCell,iEdge) = dvEdge(iEdge)* & + advCoefs3rd(iCell,iEdge) + end do + end if ! only do for edges of owned-cells + end do ! end loop over edges + + deallocate(cellIndx, & + cellIndxSorted, & + derivTwo) + + ! If 2nd order advection, disable high-order terms by + ! setting mask to zero. + if (config_horiz_tracer_adv_order == 2) then + advMaskHighOrder(:,:) = 0.0_RKIND + endif + + ! Copy module variables to device + !$acc enter data copyin(nAdvCellsForEdge, advCellsForEdge, & + !$acc advCoefs, advCoefs3rd, advMaskHighOrder) + + !-------------------------------------------------------------------- + + end subroutine li_tracer_advection_fct_shared_init !}}} + +!*********************************************************************** +! +! routine computeDerivTwo +! +!> \brief MPAS deriv two computation +!> \author Doug Jacobsen, Bill Skamarock +!> \date 03/09/12 +!> \details +!> This routine precomputes the second derivative values for tracer +!> advection. It computes cell coefficients for the polynomial fit +!> as described in: +!> Skamarock, W. C., & Gassmann, A. (2011). +!> Conservative Transport Schemes for Spherical Geodesic Meshs: +!> High-Order Flux Operators for ODE-Based Time Integration. +!> Monthly Weather Review, 139(9), 2962-2975. +!> doi:10.1175/MWR-D-10-05056.1 +!> This is performed during model initialization. +! +!----------------------------------------------------------------------- + + subroutine computeDerivTwo(derivTwo, err)!{{{ + + !----------------------------------------------------------------- + ! Output variables + !----------------------------------------------------------------- + + real (kind=RKIND), dimension(:,:,:), intent(out) :: & + derivTwo !< [out] 2nd deriv values of polynomial for tracer fit + + integer, intent(out) :: err !< [out] returned error flag + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + logical :: & + addCell, &! flag for adding cell add to list + doCell ! flag for whether to process cell + + integer, parameter :: & + polynomial_order = 2 ! option for polynomial order, forced to 2 + + integer :: & + i, j, n, &! loop counters + ip1, ip2, &! temps for i+1, i+2 + iCell, iEdge, &! loop indices for cell, edge loops + ma, na, cell_add, mw + + real (kind=RKIND) :: & + xec, yec, zec, &! arc bisection coords + thetae_tmp, &! angle + xv1, xv2, yv1, yv2, zv1, zv2, &! vertex cart coords + pii, &! pi math constant + length_scale, &! length scale + cos2t, costsint, sin2t ! trig function temps + + real (kind=RKIND), dimension(:,:), allocatable :: & + thetae ! sphere angle between vectors + + real (kind=RKIND), dimension(:), allocatable :: & + theta_abs, &! sphere angle + angle_2d ! 2d angle + + integer, dimension(25) :: cell_list + + real (kind=RKIND), dimension(25) :: & + xc, yc, zc, &! cell center coordinates + xp, yp, &! + thetav, thetat, dl_sphere + + real (kind=RKIND) :: & + amatrix(25,25), bmatrix(25,25), wmatrix(25,25) + + logical, parameter :: onSphere=.false. + + ! End preamble + !------------- + ! Begin code + + ! set return code and check proper poly order + err = 0 + + if (polynomial_order > 2) then + call mpas_log_write( & + 'Polynomial for second derivitave can only be 2', & + MPAS_LOG_ERR) + err = 1 + return + end if + + ! allocate local arrays + allocate(thetae(2, nEdges)) + allocate(theta_abs(nCells)) + allocate(angle_2d(maxEdges)) + + ! Initialize derivTwo and pi + + pii = 2.*asin(1.0) + derivTwo(:,:,:) = 0.0_RKIND + + do iCell = 1, nCells + + cell_list(1) = iCell + do i=2, nEdgesOnCell(iCell)+1 + cell_list(i) = cellsOnCell(i-1,iCell) + end do + n = nEdgesOnCell(iCell) + 1 + + !if ( polynomial_order > 2 ) then + ! do i=2, nEdgesOnCell(iCell) + 1 + ! do j=1, nEdgesOnCell( cell_list(i) ) + ! cell_add = cellsOnCell(j,cell_list(i)) + ! addCell = .true. + ! do k=1,n + ! if ( cell_add == cell_list(k) ) addCell = .false. + ! end do + ! if (addCell) then + ! n = n+1 + ! cell_list(n) = cell_add + ! end if + ! end do + ! end do + !end if + + ! check to see if we are reaching outside the halo + + doCell = .true. + do i=1,n + if (cell_list(i) > nCells) doCell = .false. + end do + + if ( .not. doCell ) cycle + + ! compute poynomial fit for this cell if all needed + ! neighbors exist + if ( onSphere ) then + + do i=1,n + j = cell_list(i) + xc(i) = xCell(j) / sphereRadius + yc(i) = yCell(j) / sphereRadius + zc(i) = zCell(j) / sphereRadius + end do + + if ( zc(1) == 1.0_RKIND) then + theta_abs(iCell) = pii/2.0_RKIND + else + theta_abs(iCell) = pii/2.0_RKIND & + - mpas_sphere_angle( xc(1), yc(1), zc(1), & + xc(2), yc(2), zc(2), & + 0.0_RKIND, 0.0_RKIND, 1.0_RKIND) + end if + + ! angles from cell center to neighbor centers (thetav) + + do i=1,n-1 + + ip2 = i+2 + if (ip2 > n) ip2 = 2 + + thetav(i) = mpas_sphere_angle(xc( 1), yc( 1), zc( 1),& + xc(i+1), yc(i+1), zc(i+1),& + xc(ip2), yc(ip2), zc(ip2)) + + dl_sphere(i) = sphereRadius* & + mpas_arc_length(xc( 1),yc( 1),zc( 1), & + xc(i+1),yc(i+1),zc(i+1)) + end do + + length_scale = 1.0_RKIND + do i=1,n-1 + dl_sphere(i) = dl_sphere(i)/length_scale + end do + + ! thetat(1) = 0. ! this defines the x direction, cell center 1 -> + ! this defines the x direction, longitude line + thetat(1) = theta_abs(iCell) + do i=2,n-1 + thetat(i) = thetat(i-1) + thetav(i-1) + end do + + do i=1,n-1 + xp(i) = cos(thetat(i)) * dl_sphere(i) + yp(i) = sin(thetat(i)) * dl_sphere(i) + end do + + else ! On an x-y plane + + do i=1,n-1 + + angle_2d(i) = angleEdge(edgesOnCell(i,iCell)) + iEdge = edgesOnCell(i,iCell) + if ( iCell .ne. cellsOnEdge(1,iEdge)) & + angle_2d(i) = angle_2d(i) - pii + + xp(i) = dcEdge(edgesOnCell(i,iCell)) * cos(angle_2d(i)) + yp(i) = dcEdge(edgesOnCell(i,iCell)) * sin(angle_2d(i)) + + end do + + end if + + ma = n-1 + mw = nEdgesOnCell(iCell) + + bmatrix = 0. + amatrix = 0. + wmatrix = 0. + + if (polynomial_order == 2) then + na = 6 + ma = ma+1 + + amatrix(1,1) = 1. + wmatrix(1,1) = 1. + do i=2,ma + amatrix(i,1) = 1. + amatrix(i,2) = xp(i-1) + amatrix(i,3) = yp(i-1) + amatrix(i,4) = xp(i-1)**2 + amatrix(i,5) = xp(i-1) * yp(i-1) + amatrix(i,6) = yp(i-1)**2 + + wmatrix(i,i) = 1. + end do + + else if (polynomial_order == 3) then + na = 10 + ma = ma+1 + + amatrix(1,1) = 1. + wmatrix(1,1) = 1. + do i=2,ma + amatrix(i,1) = 1. + amatrix(i,2) = xp(i-1) + amatrix(i,3) = yp(i-1) + + amatrix(i,4) = xp(i-1)**2 + amatrix(i,5) = xp(i-1) * yp(i-1) + amatrix(i,6) = yp(i-1)**2 + + amatrix(i,7) = xp(i-1)**3 + amatrix(i,8) = yp(i-1) * (xp(i-1)**2) + amatrix(i,9) = xp(i-1) * (yp(i-1)**2) + amatrix(i,10) = yp(i-1)**3 + + wmatrix(i,i) = 1. + + end do + + else + na = 15 + ma = ma+1 + + amatrix(1,1) = 1. + wmatrix(1,1) = 1. + do i=2,ma + amatrix(i,1) = 1. + amatrix(i,2) = xp(i-1) + amatrix(i,3) = yp(i-1) + + amatrix(i,4) = xp(i-1)**2 + amatrix(i,5) = xp(i-1) * yp(i-1) + amatrix(i,6) = yp(i-1)**2 + + amatrix(i,7) = xp(i-1)**3 + amatrix(i,8) = yp(i-1) * (xp(i-1)**2) + amatrix(i,9) = xp(i-1) * (yp(i-1)**2) + amatrix(i,10) = yp(i-1)**3 + + amatrix(i,11) = xp(i-1)**4 + amatrix(i,12) = yp(i-1) * (xp(i-1)**3) + amatrix(i,13) = (xp(i-1)**2)*(yp(i-1)**2) + amatrix(i,14) = xp(i-1) * (yp(i-1)**3) + amatrix(i,15) = yp(i-1)**4 + + wmatrix(i,i) = 1. + + end do + + do i=1,mw + wmatrix(i,i) = 1. + end do + + end if + + call mpas_poly_fit_2( amatrix, bmatrix, wmatrix, ma, na, 25 ) + + do i=1, nEdgesOnCell(iCell) + ip1 = i+1 + if (ip1 > n-1) ip1 = 1 + + iEdge = edgesOnCell(i,iCell) + + if ( onSphere ) then + xv1 = xVertex(verticesOnEdge(1,iedge)) / sphereRadius + yv1 = yVertex(verticesOnEdge(1,iedge)) / sphereRadius + zv1 = zVertex(verticesOnEdge(1,iedge)) / sphereRadius + xv2 = xVertex(verticesOnEdge(2,iedge)) / sphereRadius + yv2 = yVertex(verticesOnEdge(2,iedge)) / sphereRadius + zv2 = zVertex(verticesOnEdge(2,iedge)) / sphereRadius + else + xv1 = xVertex(verticesOnEdge(1,iedge)) + yv1 = yVertex(verticesOnEdge(1,iedge)) + zv1 = zVertex(verticesOnEdge(1,iedge)) + xv2 = xVertex(verticesOnEdge(2,iedge)) + yv2 = yVertex(verticesOnEdge(2,iedge)) + zv2 = zVertex(verticesOnEdge(2,iedge)) + end if + + if ( onSphere ) then + call mpas_arc_bisect( xv1, yv1, zv1, & + xv2, yv2, zv2, & + xec, yec, zec ) + + thetae_tmp = mpas_sphere_angle( & + xc( 1), yc( 1), zc( 1), & + xc(i+1), yc(i+1), zc(i+1), & + xec, yec, zec ) + thetae_tmp = thetae_tmp + thetat(i) + if (iCell == cellsOnEdge(1,iEdge)) then + thetae(1, edgesOnCell(i,iCell)) = thetae_tmp + else + thetae(2, edgesOnCell(i,iCell)) = thetae_tmp + end if + end if + + end do + + ! fill second derivative stencil for rk advection + + do i=1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i,iCell) + + if ( onSphere ) then + if (iCell == cellsOnEdge(1,iEdge)) then + + cos2t = cos(thetae(1, edgesOnCell(i,iCell))) + sin2t = sin(thetae(1, edgesOnCell(i,iCell))) + costsint = cos2t*sin2t + cos2t = cos2t**2 + sin2t = sin2t**2 + + do j=1,n + derivTwo(j,1,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + + else + + cos2t = cos(thetae(2, edgesOnCell(i,iCell))) + sin2t = sin(thetae(2, edgesOnCell(i,iCell))) + costsint = cos2t*sin2t + cos2t = cos2t**2 + sin2t = sin2t**2 + + do j=1,n + derivTwo(j,2,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + end if + + else + + cos2t = cos(angle_2d(i)) + sin2t = sin(angle_2d(i)) + costsint = cos2t*sin2t + cos2t = cos2t**2 + sin2t = sin2t**2 + +! do j=1,n +! +! derivTwo(j,1,iEdge) = 2.*xe(iEdge)*xe(iEdge)*bmatrix(4,j) & +! + 2.*xe(iEdge)*ye(iEdge)*bmatrix(5,j) & +! + 2.*ye(iEdge)*ye(iEdge)*bmatrix(6,j) +! end do + + if (iCell == cellsOnEdge(1,iEdge)) then + do j=1,n + derivTwo(j,1,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + else + do j=1,n + derivTwo(j,2,iEdge) = 2.*cos2t *bmatrix(4,j) & + + 2.*costsint*bmatrix(5,j) & + + 2.*sin2t *bmatrix(6,j) + end do + end if + + end if + end do + + end do ! end of loop over cells + + deallocate(thetae) + deallocate(theta_abs) + deallocate(angle_2d) + + !-------------------------------------------------------------------- + + end subroutine computeDerivTwo !}}} + +!*********************************************************************** + +end module li_tracer_advection_fct_shared + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F index 9975963323e2..9c48a67a745e 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_core.F @@ -15,6 +15,8 @@ module li_core use li_velocity use li_setup use li_mask + use li_mesh + use li_config implicit none private @@ -121,6 +123,7 @@ function li_core_init(domain, startTimeStamp) result(err) err_tmp = 0 globalErr = 0 + call li_config_init(domain % configs) call mpas_pool_get_config(domain % configs, 'config_do_restart', config_do_restart) ! @@ -206,6 +209,11 @@ function li_core_init(domain, startTimeStamp) result(err) call mpas_timer_stop('reset_io_alarms') call mpas_log_write("Finished processing recurring input streams at initial time.") + ! === + ! Initialize land ice mesh public paramters + call li_meshCreate(domain) + ! === + ! === ! Initialize some time stuff on each block ! === @@ -254,6 +262,11 @@ function li_core_init(domain, startTimeStamp) result(err) block => block % next end do + ! === + ! This is where li_updateMeshFIelds would be called if we wanted to be consistent with mpas-ocean. + ! We don't currently think this needs to be called, so leaving it out for now. + ! === + ! === ! === Initialize modules === ! === @@ -310,8 +323,6 @@ function li_core_init(domain, startTimeStamp) result(err) ! === call li_velocity_external_write_albany_mesh(domain) - ! check for errors and exit - call mpas_dmpar_max_int(domain % dminfo, err, globalErr) ! Find out if any blocks got an error if (globalErr > 0) then call mpas_log_write("An error has occurred in li_core_init. Aborting...", MPAS_LOG_CRIT) @@ -649,6 +660,9 @@ function li_core_finalize(domain) result(err) err = ior(err, err_tmp) + call li_meshDestroy(err_tmp) + err = ior(err, err_tmp) + ! === error check and exit call mpas_dmpar_max_int(domain % dminfo, err, globalErr) ! Find out if any blocks got an error if (globalErr > 0) then diff --git a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F index 2bbc58d3f6f8..803fbaa79d5f 100644 --- a/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F +++ b/components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration_fe.F @@ -661,7 +661,8 @@ subroutine advection_solver(domain, err) ! === ! === Calculate layerThicknessEdge, which is needed for advection ! === - if (trim(config_thickness_advection) == 'fo' .or. trim(config_tracer_advection) == 'fo') then + if (trim(config_thickness_advection) == 'fo' .or. trim(config_tracer_advection) == 'fo' .or. & + trim(config_thickness_advection) == 'fct' .or. trim(config_tracer_advection) == 'fct') then block => domain % blocklist do while (associated(block)) @@ -697,7 +698,9 @@ subroutine advection_solver(domain, err) call mpas_pool_get_subpool(block % structs, 'scratch', scratchPool) call mpas_pool_get_array(meshPool, 'deltat', deltat) - if (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fo') then + if ( (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fo') .or. & + (trim(config_thickness_advection) == 'fct') .or. & + (trim(config_thickness_advection) == 'fo' .and. trim(config_tracer_advection) == 'fct') ) then ! Note: This subroutine requires that thickness and tracers are correct in halos @@ -752,7 +755,6 @@ subroutine advection_solver(domain, err) call mpas_pool_get_dimension(meshPool, 'nCells', nCells) call mpas_pool_get_array(geometryPool, 'thickness', thickness, timeLevel=1) - allocate( masktmp(nCells + 1) ) masktmp = 0 @@ -789,6 +791,7 @@ subroutine advection_solver(domain, err) 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') call mpas_timer_stop("halo updates") diff --git a/components/mpas-albany-landice/src/shared/Makefile b/components/mpas-albany-landice/src/shared/Makefile index b9e1f6820813..e0f1828298d9 100644 --- a/components/mpas-albany-landice/src/shared/Makefile +++ b/components/mpas-albany-landice/src/shared/Makefile @@ -3,6 +3,8 @@ OBJS = mpas_li_constants.o \ mpas_li_mask.o \ + mpas_li_mesh.o \ + mpas_li_config.o \ mpas_li_setup.o all: $(OBJS) @@ -13,7 +15,9 @@ mpas_li_setup.o: mpas_li_mask.o: mpas_li_setup.o +mpas_li_mesh.o: +mpas_li_config.o: clean: $(RM) *.o *.mod *.f90 diff --git a/components/mpas-albany-landice/src/shared/mpas_li_config.F b/components/mpas-albany-landice/src/shared/mpas_li_config.F new file mode 100644 index 000000000000..771c89bc6afc --- /dev/null +++ b/components/mpas-albany-landice/src/shared/mpas_li_config.F @@ -0,0 +1,55 @@ +! 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 +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! li_config +! +!> \brief MALI specific config +!> \details +!> This module contains config specific to MALI. +! +!----------------------------------------------------------------------- + +module li_config + + use mpas_derived_types + use mpas_pool_routines + use mpas_kind_types + + implicit none + public + save + +#include "../inc/config_declare.inc" + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! routine li_config_init +! +!> \brief Initializes the MALI config +!> \details +!> This routine sets up config for use in MALI. +! +!----------------------------------------------------------------------- + subroutine li_config_init(configPool)!{{{ + type (mpas_pool_type), pointer :: configPool + +#include "../inc/config_get.inc" + + end subroutine li_config_init!}}} + +!*********************************************************************** + +end module li_config + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-albany-landice/src/shared/mpas_li_mesh.F b/components/mpas-albany-landice/src/shared/mpas_li_mesh.F new file mode 100644 index 000000000000..b59497b155e3 --- /dev/null +++ b/components/mpas-albany-landice/src/shared/mpas_li_mesh.F @@ -0,0 +1,1768 @@ +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! \file mpas_li_mesh.F +! +! 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 +! +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! +! li_mesh +! +!> \brief MALI mesh structure with GPU support +!> \author Rob Aulwes and Phil Jones, modified for MALI +!> by Trevor Hillebrand and Matthew Hoffman +!> \date 14 Jan 2020 +!> \details +!> This module creates and maintains a primary land ice mesh structure +!> and ensures all mesh variables are copied to an accelerator device +!> if needed. Currently it consists of pointers to the existing MPAS mesh pool +!> variables, but is intended to eventually replace the mesh pool later. +! +!------------------------------------------------------------------------------- + +module li_mesh + + ! module dependencies + use mpas_dmpar + use mpas_derived_types + use mpas_pool_routines + use mpas_constants + use mpas_log + +! use ocn_config + + implicit none + private + + !---------------------------------------------------------------------------- + ! + ! Public parameters + ! + !---------------------------------------------------------------------------- + !{{{ + + logical, public :: & + onSphere ! this mesh is on the sphere + + real (kind=RKIND), public :: & + sphereRadius ! radius of sphere for spherical meshes + + integer, public :: &! mesh, array sizes + nCells, &! total number of local (owned+halo) cells in primary + nEdges, &! total number of local edge midpoints + nVertices, &! total number of local cells in dual (cell vertices) + nCellsSolve, &! number of cells owned by the local domain + nEdgesSolve, &! number of edges owned by the local domain + nVerticesSolve, &! number of vertices owned by the local domain + maxEdges, &! largest number of edges any polygon has + maxEdges2, &! 2x the largest number of edges any polygon has + vertexDegree, &! number of cells or edges touching each vertex + nVertLevels ! number of vertical levels +! nVertLevelsP1 ! number of vertical interfaces (levels plus one) + + integer, public, dimension(:), allocatable :: & + nCellsHalo, &! number of owned+halo(n) cells in local domain + nEdgesHalo, &! number of owned+halo(n) edges in local domain + nVerticesHalo ! number of owned+halo(n) vertices in local domain + + integer, public, dimension(:), pointer :: & + nEdgesOnEdge, &! number of edges connected to each edge point + nEdgesOnCell, & ! number of edges associated with each cell center +! minLevelCell, &! max ocean level at bottom of cell +! minLevelEdgeTop, &! max ocean level at top of edge column +! minLevelEdgeBot, &! max ocean level at bottom of edge column +! minLevelVertexTop, &! max ocean level at top of each vertex +! minLevelVertexBot, &! max ocean level at bottom of each vertex +! maxLevelCell, &! max ocean level at bottom of cell +! maxLevelEdgeTop, &! max ocean level at top of edge column +! maxLevelEdgeBot, &! max ocean level at bottom of edge column +! maxLevelVertexTop, &! max ocean level at top of each vertex +! maxLevelVertexBot, &! max ocean level at bottom of each vertex + indexToCellID, &! global ID of each local cell + indexToEdgeID, &! global ID of each local edge + indexToVertexID ! global ID of each local vertex + + integer, public, dimension(:, :), pointer :: & + edgesOnEdge, &! index of edges connected to each edge + cellsOnEdge, &! index of cells connected to each edge + verticesOnEdge, &! index of vertices connected to each edge + cellsOnCell, &! index of cells connected to each cell + edgesOnCell, &! index of edges connected to each cell + verticesOnCell, &! index of vertices connected to each cell + cellsOnVertex, &! index of cells connected to each vertex + edgesOnVertex, & ! index of edges connected to each vertex + edgeSignOnCell, &! sign of edge contributions to a cell + edgeSignOnVertex ! sign of edge contributions to a vertex +! kiteIndexOnCell ! index of kite associated with each cell + + real(kind=RKIND), public, dimension(:), pointer :: & + latCell, &! latitude of cell centers + lonCell, &! longitude of cell centers + xCell, &! Cartesian x coord of cell center + yCell, &! Cartesian y coord of cell center + zCell, &! Cartesian z coord of cell center + latEdge, &! latitude of edge + lonEdge, &! longitude of edge + xEdge, &! Cartesian x coord of edge + yEdge, &! Cartesian y coord of edge + zEdge, &! Cartesian z coord of edge + latVertex, &! latitude of vertex + lonVertex, &! longitude of vertex + xVertex, &! Cartesian coord of vertex + yVertex, &! Cartesian y coord of vertex + zVertex, &! Cartesian z coord of vertex +! fEdge, &! Coriolus parameter at edge +! fVertex, &! Coriolus parameter at vertex +! fCell, &! Coriolus parameter at cell center + dcEdge, &! length of edge = dist between cells across edge + dvEdge, &! length of edge = dist between vertices along edge + areaCell, &! area of each cell + areaTriangle, &! area of each cell on dual grid +! bottomDepth, &! ocean bottom depth at each cell center +! refBottomDepth, &! ocean depth at bottom of cell for reference profile +! refBottomDepthTopOfCell, &! depth at top of cell for reference profile +! vertCoordMovementWeights, &! weights for distributing height perturb +! meshScalingDel2, &! mesh scaling factor for use in del2 diffusion +! meshScalingDel4, &! mesh scaling factor for use in del4 diffusion + meshDensity, &! density of mesh + angleEdge ! angle the edge normal makes with local east +! distanceToCoast ! distance to nearest coast cell + + real(kind=RKIND), public, dimension(:), allocatable :: & + invAreaCell ! inverse of area of each cell + + ! Multiplicative masks and vectors for various conditions + integer, public, dimension(:), allocatable :: & + boundaryEdge, &! mask for boundary edges at each level + boundaryCell, &! mask for boundary cells at each level + boundaryVertex ! mask for boundary vertices at each level + + real(kind=RKIND), public, dimension(:, :), pointer :: & + weightsOnEdge, &! weights on each edge + kiteAreasOnVertex, &! real (vertexDegree nVertices) + edgeTangentVectors, &! tangent unit vector at edge + edgeNormalVectors, &! normal unit vector at edge + localVerticalUnitVectors ! local unit vector iin vertical + + real(kind=RKIND), public, dimension(:, :, :), pointer :: & + cellTangentPlane, &! two vectors defining tangent plane at cell center + coeffs_reconstruct ! coeffs for reconstructing vectors at cell centers + + ! Derived mesh values needed + real(kind=RKIND), public, dimension(:, :), allocatable :: & + edgeAreaFractionOfCell + + !}}} + + !---------------------------------------------------------------------------- + ! + ! Public member functions + ! + !---------------------------------------------------------------------------- + !{{{ + + public :: & + li_meshCreate, & + li_meshUpdateFields, & + li_meshDestroy + !}}} + +!*********************************************************************** + +contains + +!*********************************************************************** +! +! li_meshCreate +! +!> \brief Creates the ocean mesh data structure on both host and device +!> \author Rob Aulwes and Phil Jones +!> \date 14 Jan 2020 +!> \details +!> This module creates and maintains public ocean mesh data +!> and ensures all mesh variables are copied to an accelerator device +!> if needed. +! +!----------------------------------------------------------------------- + + subroutine li_meshCreate(domain) !{{{ + + ! Input arguments + + type(domain_type) :: & + domain !< [in] MPAS type to describe domain + + ! Local variables + + integer :: & + blockCount ! counter for number of blocks + + type(block_type), pointer :: & + block ! variables in current subblock + + type(mpas_pool_type), pointer :: & + meshPool ! mesh variables in MPAS pool structure + + real (kind=RKIND) :: & + maxDensityLocal, maxDensityGlobal ! temps for mesh density + + ! scalar pointers for retrieval, but convert to actual scalars in struct + logical, pointer :: & + on_a_sphere + + real (kind=RKIND), pointer :: & + sphere_radius + + integer, pointer :: &! mesh dimensions + nCellsTmp, &! + nEdgesTmp, &! + nVerticesTmp, &! + maxEdgesTmp, &! + maxEdges2Tmp, &! + vertexDegreeTmp, &! + nVertLevelsTmp ! + !nVertLevelsP1Tmp ! + + ! temporary pointers for converting index arrays + integer, dimension(:), pointer :: & + nCellsArrayTmp, & + nEdgesArrayTmp, & + nVerticesArrayTmp + + ! temporary pointers for converting masks + integer i, k, n ! loop indices + integer, dimension(:, :), pointer :: & + edgeSignOnCellTmp, & + edgeSignOnVertexTmp + + !*** + !*** end of preamble, begin code + !*** + + ! We only support one block so test for condition here + blockCount = 0 + block => domain%blocklist + do while (associated(block)) + blockCount = blockCount + 1 + if (blockCount > 1) then + call mpas_log_write( & + 'li_meshCreate: more than one block no longer supported', & + MPAS_LOG_CRIT) + endif + block => block%next + end do + + ! Reset to original block + block => domain%blocklist + + ! retrieve the mpas mesh pool + call mpas_pool_get_subpool(block%structs, 'mesh', meshPool) + + !----------------------------------------------------------------- + ! first set pointers/values for all mesh variables + ! many variables already initialized based on read of mesh file + !----------------------------------------------------------------- + + ! set all mesh properties + call mpas_pool_get_config(meshPool, 'on_a_sphere', & + on_a_sphere) + call mpas_pool_get_config(meshPool, 'sphere_radius', & + sphere_radius) + + ! set all mesh dimensions + call mpas_pool_get_dimension(meshPool, 'nCells', & + nCellsTmp) + call mpas_pool_get_dimension(meshPool, 'nEdges', & + nEdgesTmp) + call mpas_pool_get_dimension(meshPool, 'nVertices', & + nVerticesTmp) + call mpas_pool_get_dimension(meshPool, 'maxEdges', & + maxEdgesTmp) + call mpas_pool_get_dimension(meshPool, 'maxEdges2', & + maxEdges2Tmp) + call mpas_pool_get_dimension(meshPool, 'vertexDegree', & + vertexDegreeTmp) + call mpas_pool_get_dimension(meshPool, 'nVertLevels', & + nVertLevelsTmp) +! call mpas_pool_get_dimension(meshPool, 'nVertLevelsP1', & +! nVertLevelsP1Tmp) + call mpas_pool_get_dimension(meshPool, 'nCellsArray', & + nCellsArrayTmp) + call mpas_pool_get_dimension(meshPool, 'nEdgesArray', & + nEdgesArrayTmp) + call mpas_pool_get_dimension(meshPool, 'nVerticesArray', & + nVerticesArrayTmp) + + ! translate scalar pointers to scalars in new mesh structure + onSphere = on_a_sphere + sphereRadius = sphere_radius + maxEdges = maxEdgesTmp + maxEdges2 = maxEdges2Tmp + vertexDegree = vertexDegreeTmp + nVertLevels = nVertLevelsTmp +! nVertLevelsP1 = nVertLevelsP1Tmp + + ! convert previous index limits into new halo definitions + nCells = nCellsTmp + nEdges = nEdgesTmp + nVertices = nVerticesTmp + + n = size(nCellsArrayTmp) + allocate (nCellsHalo(n - 1)) + nCellsSolve = nCellsArrayTmp(1) + do i = 2, n + nCellsHalo(i - 1) = nCellsArrayTmp(i) + end do + + n = size(nEdgesArrayTmp) + allocate (nEdgesHalo(n - 1)) + nEdgesSolve = nEdgesArrayTmp(1) + do i = 2, n + nEdgesHalo(i - 1) = nEdgesArrayTmp(i) + end do + + n = size(nVerticesArrayTmp) + allocate (nVerticesHalo(n - 1)) + nVerticesSolve = nVerticesArrayTmp(1) + do i = 2, n + nVerticesHalo(i - 1) = nVerticesArrayTmp(i) + end do + + ! set pointers for a lot of connectivity info + call mpas_pool_get_array(meshPool, 'nEdgesOnEdge', & + nEdgesOnEdge) + call mpas_pool_get_array(meshPool, 'nEdgesOnCell', & + nEdgesOnCell) +! call mpas_pool_get_array(meshPool, 'minLevelCell', & +! minLevelCell) +! call mpas_pool_get_array(meshPool, 'minLevelEdgeTop', & +! minLevelEdgeTop) +! call mpas_pool_get_array(meshPool, 'minLevelEdgeBot', & +! minLevelEdgeBot) +! call mpas_pool_get_array(meshPool, 'minLevelVertexTop', & +! minLevelVertexTop) +! call mpas_pool_get_array(meshPool, 'minLevelVertexBot', & +! minLevelVertexBot) +! call mpas_pool_get_array(meshPool, 'maxLevelCell', & +! maxLevelCell) +! call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', & +! maxLevelEdgeTop) +! call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', & +! maxLevelEdgeBot) +! call mpas_pool_get_array(meshPool, 'maxLevelVertexTop', & +! maxLevelVertexTop) +! call mpas_pool_get_array(meshPool, 'maxLevelVertexBot', & +! maxLevelVertexBot) + call mpas_pool_get_array(meshPool, 'indexToCellID', & + indexToCellID) + call mpas_pool_get_array(meshPool, 'indexToEdgeID', & + indexToEdgeID) + call mpas_pool_get_array(meshPool, 'indexToVertexID', & + indexToVertexID) + call mpas_pool_get_array(meshPool, 'edgesOnEdge', & + edgesOnEdge) + call mpas_pool_get_array(meshPool, 'cellsOnEdge', & + cellsOnEdge) + call mpas_pool_get_array(meshPool, 'verticesOnEdge', & + verticesOnEdge) + call mpas_pool_get_array(meshPool, 'cellsOnCell', & + cellsOnCell) + call mpas_pool_get_array(meshPool, 'edgesOnCell', & + edgesOnCell) + call mpas_pool_get_array(meshPool, 'verticesOnCell', & + verticesOnCell) + call mpas_pool_get_array(meshPool, 'cellsOnVertex', & + cellsOnVertex) + call mpas_pool_get_array(meshPool, 'edgesOnVertex', & + edgesOnVertex) +! call mpas_pool_get_array(meshPool, 'kiteIndexOnCell', & +! kiteIndexOnCell) + + ! now set a number of physics and numerical properties of mesh + call mpas_pool_get_array(meshPool, 'latCell', & + latCell) + call mpas_pool_get_array(meshPool, 'lonCell', & + lonCell) + call mpas_pool_get_array(meshPool, 'xCell', & + xCell) + call mpas_pool_get_array(meshPool, 'yCell', & + yCell) + call mpas_pool_get_array(meshPool, 'zCell', & + zCell) + call mpas_pool_get_array(meshPool, 'latEdge', & + latEdge) + call mpas_pool_get_array(meshPool, 'lonEdge', & + lonEdge) + call mpas_pool_get_array(meshPool, 'xEdge', & + xEdge) + call mpas_pool_get_array(meshPool, 'yEdge', & + yEdge) + call mpas_pool_get_array(meshPool, 'zEdge', & + zEdge) + call mpas_pool_get_array(meshPool, 'latVertex', & + latVertex) + call mpas_pool_get_array(meshPool, 'lonVertex', & + lonVertex) + call mpas_pool_get_array(meshPool, 'xVertex', & + xVertex) + call mpas_pool_get_array(meshPool, 'yVertex', & + yVertex) + call mpas_pool_get_array(meshPool, 'zVertex', & + zVertex) +! call mpas_pool_get_array(meshPool, 'fEdge', & +! fEdge) +! call mpas_pool_get_array(meshPool, 'fVertex', & +! fVertex) +! call mpas_pool_get_array(meshPool, 'fCell', & +! fCell) + call mpas_pool_get_array(meshPool, 'dcEdge', & + dcEdge) + call mpas_pool_get_array(meshPool, 'dvEdge', & + dvEdge) + call mpas_pool_get_array(meshPool, 'areaCell', & + areaCell) + call mpas_pool_get_array(meshPool, 'areaTriangle', & + areaTriangle) + call mpas_pool_get_array(meshPool, 'weightsOnEdge', & + weightsOnEdge) +! call mpas_pool_get_array(meshPool, 'bottomDepth', & +! bottomDepth) +! call mpas_pool_get_array(meshPool, 'refBottomDepth', & +! refBottomDepth) +! call mpas_pool_get_array(meshPool, 'refBottomDepthTopOfCell', & +! refBottomDepthTopOfCell) +! call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', & +! vertCoordMovementWeights) +! call mpas_pool_get_array(meshPool, 'meshScalingDel2', & +! meshScalingDel2) +! call mpas_pool_get_array(meshPool, 'meshScalingDel4', & +! meshScalingDel4) + call mpas_pool_get_array(meshPool, 'meshDensity', & + meshDensity) + call mpas_pool_get_array(meshPool, 'angleEdge', & + angleEdge) +! call mpas_pool_get_array(meshPool, 'distanceToCoast', & +! distanceToCoast) + + call mpas_pool_get_array(meshPool, 'weightsOnEdge', & + weightsOnEdge) + call mpas_pool_get_array(meshPool, 'kiteAreasOnVertex', & + kiteAreasOnVertex) + call mpas_pool_get_array(meshPool, 'edgeTangentVectors', & + edgeTangentVectors) + call mpas_pool_get_array(meshPool, 'edgeNormalVectors', & + edgeNormalVectors) + call mpas_pool_get_array(meshPool, 'localVerticalUnitVectors', & + localVerticalUnitVectors) + call mpas_pool_get_array(meshPool, 'cellTangentPlane', & + cellTangentPlane) + call mpas_pool_get_array(meshPool, 'coeffs_reconstruct', & + coeffs_reconstruct) + + ! For masks, we wish to convert to real multiplicative masks + ! so retrieve integer version pointers and allocate real masks + ! Once these are converted in Registry, we can eliminate this. + call mpas_pool_get_array(meshPool, 'edgeSignOnCell', & + edgeSignOnCell) + call mpas_pool_get_array(meshPool, 'edgeSignOnVertex', & + edgeSignOnVertex) + + allocate ( & + boundaryEdge(nEdges+1), & + boundaryCell(nCells+1), & + boundaryVertex(nVertices+1), & + invAreaCell(nCells+1)) + + !----------------------------------------------------------------- + ! Now that all pointers are set and mesh variables allocated + ! we initialize other mesh quantities + !----------------------------------------------------------------- + + ! Start by initializing vertical mesh, min/max cells and + ! sign/index fields + !call meshVertCoordInit(domain) + call meshMinMaxLevel() + call meshSignIndexFields() + + ! Compute max mesh density and mesh scaling +! if (config_maxMeshDensity < 0.0_RKIND) then +! maxDensityLocal = maxval(meshDensity) +! call mpas_dmpar_max_real(domain%dminfo, maxDensityLocal, & +! maxDensityGlobal) +! config_maxMeshDensity = maxDensityGlobal +! endif +! call meshScaling() + + ! Routine calls above initialized the real mask arrays + ! so convert the meshPool integer pointers here +! do n = 1, nCells+1 +! do k = 1, nVertLevels +! boundaryCellTmp(k,n) = nint(boundaryCell(k,n)) +! end do +! end do +! +! do n = 1, nCells+1 +! do k = 1, maxEdges +! edgeSignOnCellTmp(k, n) = nint(edgeSignOnCell(k,n)) +! end do +! end do +! +! do n = 1, nEdges+1 +! do k = 1, nVertLevels +! boundaryEdgeTmp(k, n) = nint(boundaryEdge(k,n)) +! end do +! end do +! +! do n = 1, nVertices+1 +! do k = 1, nVertLevels +! boundaryVertexTmp(k, n) = nint(boundaryVertex(k,n)) +! end do +! end do +! +! do n = 1, nVertices+1 +! do k = 1, vertexDegree +! edgeSignOnVertexTmp(k,n) = nint(edgeSignOnVertex(k,n)) +! end do +! end do + + ! Compute area weighting for some interpolation + call meshAreaWeights() + areaCell(nCells+1) = -1.0e34_RKIND + + ! Compute the inverse of areaCell + do n = 1, nCells + invAreaCell(n) = 1.0_RKIND / areaCell(n) + end do + invAreaCell(nCells+1) = 0.0_RKIND + + !----------------------------------------------------------------- + ! Copy mesh data to accelerator if needed. + !----------------------------------------------------------------- + + !$acc enter data copyin( & + !$acc onSphere, & + !$acc sphereRadius, & + !$acc nCells, & + !$acc nEdges, & + !$acc nVertices, & + !$acc nCellsSolve, & + !$acc nEdgesSolve, & + !$acc nVerticesSolve, & + !$acc maxEdges, & + !$acc maxEdges2, & + !$acc vertexDegree, & + !$acc nVertLevels, & + !$acc nVertLevelsP1, & + !$acc nEdgesHalo, & + !$acc nCellsHalo, & + !$acc nVerticesHalo, & + !$acc nEdgesOnEdge, & + !$acc nEdgesOnCell, & + !$acc minLevelCell, & + !$acc minLevelEdgeTop, & + !$acc minLevelEdgeBot, & + !$acc minLevelVertexTop, & + !$acc minLevelVertexBot, & + !$acc maxLevelCell, & + !$acc maxLevelEdgeTop, & + !$acc maxLevelEdgeBot, & + !$acc maxLevelVertexTop, & + !$acc maxLevelVertexBot, & + !$acc indexToCellID, & + !$acc indexToEdgeID, & + !$acc indexToVertexID, & + !$acc edgesOnEdge, & + !$acc cellsOnEdge, & + !$acc verticesOnEdge, & + !$acc cellsOnCell, & + !$acc edgesOnCell, & + !$acc verticesOnCell, & + !$acc cellsOnVertex, & + !$acc edgesOnVertex, & + !$acc kiteIndexOnCell, & + !$acc latCell, & + !$acc lonCell, & + !$acc xCell, & + !$acc yCell, & + !$acc zCell, & + !$acc latEdge, & + !$acc lonEdge, & + !$acc xEdge, & + !$acc yEdge, & + !$acc zEdge, & + !$acc latVertex, & + !$acc lonVertex, & + !$acc xVertex, & + !$acc yVertex, & + !$acc zVertex, & + !$acc fEdge, & + !$acc fVertex, & + !$acc fCell, & + !$acc dcEdge, & + !$acc dvEdge, & + !$acc areaCell, & + !$acc invAreaCell, & + !$acc areaTriangle, & + !$acc bottomDepth, & + !$acc refBottomDepth, & + !$acc refBottomDepthTopOfCell, & + !$acc vertCoordMovementWeights, & + !$acc meshScalingDel2, & + !$acc meshScalingDel4, & + !$acc meshDensity, & + !$acc angleEdge, & + !$acc distanceToCoast, & + !$acc boundaryEdge, & + !$acc boundaryCell, & + !$acc boundaryVertex, & + !$acc edgeSignOnCell, & + !$acc edgeSignOnVertex, & + !$acc weightsOnEdge, & + !$acc kiteAreasOnVertex, & + !$acc edgeTangentVectors, & + !$acc edgeNormalVectors, & + !$acc localVerticalUnitVectors, & + !$acc cellTangentPlane, & + !$acc coeffs_reconstruct) + +!------------------------------------------------------------------------------- + + end subroutine li_meshCreate !}}} + +!******************************************************************************* +! +! li_meshDestroy +! +!> \brief Destroy mesh structure and removes from device +!> \author Rob Aulwes and Phil Jones +!> \date 14 Jan 2020 +!> \details +!> This module removes the mesh variables from the device and invalidates +!> all pointers in the mesh structure. +! +!------------------------------------------------------------------------------- + + subroutine li_meshDestroy(err) !{{{ + + ! Input variables + + ! Since the ocnMesh is currently a public module variable, no inputs + ! here, but eventually may want to treat ocnMesh as a specific + ! instantiation instead and pass via args everywhere. If so, need an + ! input mesh here + + ! Output variables + + integer, intent(out) :: & + err ! returned error flag + + ! Local variables + + !*** + !*** end of preamble, begin code + !*** + + err = 0 + + ! First remove data from the device. Must remove components first, + ! then remove the mesh type itself. + + !$acc exit data delete(onSphere, & + !$acc sphereRadius, & + !$acc nCells, & + !$acc nEdges, & + !$acc nVertices, & + !$acc nCellsSolve, & + !$acc nEdgesSolve, & + !$acc nVerticesSolve, & + !$acc maxEdges, & + !$acc maxEdges2, & + !$acc vertexDegree, & + !$acc nVertLevels, & + !$acc nVertLevelsP1, & + !$acc nEdgesHalo, & + !$acc nCellsHalo, & + !$acc nVerticesHalo, & + !$acc nEdgesOnEdge, & + !$acc nEdgesOnCell, & + !$acc minLevelCell, & + !$acc minLevelEdgeTop, & + !$acc minLevelEdgeBot, & + !$acc minLevelVertexTop, & + !$acc minLevelVertexBot, & + !$acc maxLevelCell, & + !$acc maxLevelEdgeTop, & + !$acc maxLevelEdgeBot, & + !$acc maxLevelVertexTop, & + !$acc maxLevelVertexBot, & + !$acc indexToCellID, & + !$acc indexToEdgeID, & + !$acc indexToVertexID, & + !$acc edgesOnEdge, & + !$acc cellsOnEdge, & + !$acc verticesOnEdge, & + !$acc cellsOnCell, & + !$acc edgesOnCell, & + !$acc verticesOnCell, & + !$acc cellsOnVertex, & + !$acc edgesOnVertex, & + !$acc kiteIndexOnCell, & + !$acc latCell, & + !$acc lonCell, & + !$acc xCell, & + !$acc yCell, & + !$acc zCell, & + !$acc latEdge, & + !$acc lonEdge, & + !$acc xEdge, & + !$acc yEdge, & + !$acc zEdge, & + !$acc latVertex, & + !$acc lonVertex, & + !$acc xVertex, & + !$acc yVertex, & + !$acc zVertex, & + !$acc fEdge, & + !$acc fVertex, & + !$acc fCell, & + !$acc dcEdge, & + !$acc dvEdge, & + !$acc areaCell, & + !$acc invAreaCell, & + !$acc areaTriangle, & + !$acc bottomDepth, & + !$acc refBottomDepth, & + !$acc refBottomDepthTopOfCell, & + !$acc vertCoordMovementWeights, & + !$acc meshScalingDel2, & + !$acc meshScalingDel4, & + !$acc meshDensity, & + !$acc angleEdge, & + !$acc distanceToCoast, & + !$acc boundaryEdge, & + !$acc boundaryCell, & + !$acc boundaryVertex, & + !$acc edgeSignOnCell, & + !$acc edgeSignOnVertex, & + !$acc weightsOnEdge, & + !$acc kiteAreasOnVertex, & + !$acc edgeTangentVectors,& + !$acc edgeNormalVectors, & + !$acc localVerticalUnitVectors, & + !$acc cellTangentPlane, & + !$acc coeffs_reconstruct) + + ! Reset all scalars to zero + onSphere = .false. + sphereRadius = 0.0_RKIND + nCells = 0 + nEdges = 0 + nVertices = 0 + nCellsSolve = 0 + nEdgesSolve = 0 + nVerticesSolve = 0 + maxEdges = 0 + maxEdges2 = 0 + vertexDegree = 0 + nVertLevels = 0 +! nVertLevelsP1 = 0 + + ! Now nullify all pointers to invalidate fields + ! If this becomes the only mesh structure and mesh pool is eliminated, + ! then we will want to deallocate here instead of nullify. + + nullify (nEdgesOnEdge, & + nEdgesOnCell, & + !minLevelCell, & + !minLevelEdgeTop, & + !minLevelEdgeBot, & + !minLevelVertexTop, & + !minLevelVertexBot, & + !maxLevelCell, & + !maxLevelEdgeTop, & + !maxLevelEdgeBot, & + !maxLevelVertexTop, & + !maxLevelVertexBot, & + indexToCellID, & + indexToEdgeID, & + indexToVertexID, & + edgesOnEdge, & + cellsOnEdge, & + verticesOnEdge, & + cellsOnCell, & + edgesOnCell, & + verticesOnCell, & + cellsOnVertex, & + edgesOnVertex, & + !kiteIndexOnCell, & + latCell, & + lonCell, & + xCell, & + yCell, & + zCell, & + latEdge, & + lonEdge, & + xEdge, & + yEdge, & + zEdge, & + latVertex, & + lonVertex, & + xVertex, & + yVertex, & + zVertex, & + !fEdge, & + !fVertex, & + !fCell, & + dcEdge, & + dvEdge, & + areaCell, & + areaTriangle, & + !bottomDepth, & + !refBottomDepth, & + !refBottomDepthTopOfCell, & + !vertCoordMovementWeights, & + !meshScalingDel2, & + !meshScalingDel4, & + meshDensity, & + angleEdge, & + !distanceToCoast, & + weightsOnEdge, & + kiteAreasOnVertex, & + edgeTangentVectors, & + edgeNormalVectors, & + localVerticalUnitVectors, & + cellTangentPlane, & + edgeSignOnCell, & + edgeSignOnVertex, & + coeffs_reconstruct) + +#ifdef MPAS_OPENACC + !$acc exit data delete(edgeAreaFractionOfCell) +#endif + deallocate (nEdgesHalo, & + nCellsHalo, & + nVerticesHalo, & + boundaryEdge, & + boundaryCell, & + boundaryVertex, & + edgeAreaFractionOfCell, & + invAreaCell) + +!------------------------------------------------------------------------------- + + end subroutine li_meshDestroy !}}} + +!******************************************************************************* +! +! li_meshUpdateFields +! +!> \brief Updates fields on an accelerator device +!> \author Rob Aulwes and Phil Jones +!> \date 14 Jan 2020 +!> \details +!> Many mesh fields are computed or input later in the initialization +!> phase after the meshCreate call. This routine updates these fields +!> on the device. The routine is not needed if the original mesh pool is +!> eliminated and this mesh becomes the only mesh structure and can be +!> updated directly. +! +!------------------------------------------------------------------------------- + + subroutine li_meshUpdateFields(domain) !{{{ + + ! Input arguments + + type(domain_type) :: & + domain !< [in] MPAS type to describe domain + + ! Local variables + + type(block_type), pointer :: & + block ! variables in current subblock + + type(mpas_pool_type), pointer :: & + meshPool ! mesh variables in MPAS pool structure + + ! temporary pointers for converting masks + integer k, n ! loop indices + integer, dimension(:, :), pointer :: & + edgeSignOnCellTmp, & + edgeSignOnVertexTmp + + !*** + !*** end of preamble, begin code + !*** + + ! already check during create that there should only be one block + block => domain%blocklist + + ! retrieve the mpas mesh pool + call mpas_pool_get_subpool(block%structs, 'mesh', meshPool) + + ! Mesh dimensions should not have changed so no updates + + ! Because these fields are pointers into meshPool, they do not need + ! to be updated - they capture the updates automatically. + !call mpas_pool_get_array(meshPool, 'nEdgesOnEdge', & + ! nEdgesOnEdge) + !call mpas_pool_get_array(meshPool, 'nEdgesOnCell', & + ! nEdgesOnCell) + !call mpas_pool_get_array(meshPool, 'maxLevelCell', & + ! maxLevelCell) + !call mpas_pool_get_array(meshPool, 'maxLevelEdgeTop', & + ! maxLevelEdgeTop) + !call mpas_pool_get_array(meshPool, 'maxLevelEdgeBot', & + ! maxLevelEdgeBot) + !call mpas_pool_get_array(meshPool, 'maxLevelVertexTop', & + ! maxLevelVertexTop) + !call mpas_pool_get_array(meshPool, 'maxLevelVertexBot', & + ! maxLevelVertexBot) + !call mpas_pool_get_array(meshPool, 'indexToCellID', & + ! indexToCellID) + !call mpas_pool_get_array(meshPool, 'indexToEdgeID', & + ! indexToEdgeID) + !call mpas_pool_get_array(meshPool, 'indexToVertexID', & + ! indexToVertexID) + !call mpas_pool_get_array(meshPool, 'edgesOnEdge', & + ! edgesOnEdge) + !call mpas_pool_get_array(meshPool, 'cellsOnEdge', & + ! cellsOnEdge) + !call mpas_pool_get_array(meshPool, 'verticesOnEdge', & + ! verticesOnEdge) + !call mpas_pool_get_array(meshPool, 'cellsOnCell', & + ! cellsOnCell) + !call mpas_pool_get_array(meshPool, 'edgesOnCell', & + ! edgesOnCell) + !call mpas_pool_get_array(meshPool, 'verticesOnCell', & + ! verticesOnCell) + !call mpas_pool_get_array(meshPool, 'cellsOnVertex', & + ! cellsOnVertex) + !call mpas_pool_get_array(meshPool, 'edgesOnVertex', & + ! edgesOnVertex) + !call mpas_pool_get_array(meshPool, 'kiteIndexOnCell', & + ! kiteIndexOnCell) + + ! these are also pointers that do not require updating + ! now set a number of physics and numerical properties of mesh + !call mpas_pool_get_array(meshPool, 'latCell', & + ! latCell) + !call mpas_pool_get_array(meshPool, 'lonCell', & + ! lonCell) + !call mpas_pool_get_array(meshPool, 'xCell', & + ! xCell) + !call mpas_pool_get_array(meshPool, 'yCell', & + ! yCell) + !call mpas_pool_get_array(meshPool, 'zCell', & + ! zCell) + !call mpas_pool_get_array(meshPool, 'latEdge', & + ! latEdge) + !call mpas_pool_get_array(meshPool, 'lonEdge', & + ! lonEdge) + !call mpas_pool_get_array(meshPool, 'xEdge', & + ! xEdge) + !call mpas_pool_get_array(meshPool, 'yEdge', & + ! yEdge) + !call mpas_pool_get_array(meshPool, 'zEdge', & + ! zEdge) + !call mpas_pool_get_array(meshPool, 'latVertex', & + ! latVertex) + !call mpas_pool_get_array(meshPool, 'lonVertex', & + ! lonVertex) + !call mpas_pool_get_array(meshPool, 'xVertex', & + ! xVertex) + !call mpas_pool_get_array(meshPool, 'yVertex', & + ! yVertex) + !call mpas_pool_get_array(meshPool, 'zVertex', & + ! zVertex) + !call mpas_pool_get_array(meshPool, 'fEdge', & + ! fEdge) + !call mpas_pool_get_array(meshPool, 'fVertex', & + ! fVertex) + !call mpas_pool_get_array(meshPool, 'fCell', & + ! fCell) + !call mpas_pool_get_array(meshPool, 'dcEdge', & + ! dcEdge) + !call mpas_pool_get_array(meshPool, 'dvEdge', & + ! dvEdge) + !call mpas_pool_get_array(meshPool, 'areaCell', & + ! areaCell) + !call mpas_pool_get_array(meshPool, 'areaTriangle', & + ! areaTriangle) + !call mpas_pool_get_array(meshPool, 'weightsOnEdge', & + ! weightsOnEdge) + !call mpas_pool_get_array(meshPool, 'bottomDepth', & + ! bottomDepth) + !call mpas_pool_get_array(meshPool, 'refBottomDepth', & + ! refBottomDepth) + !call mpas_pool_get_array(meshPool, 'refBottomDepthTopOfCell', & + ! refBottomDepthTopOfCell) + !call mpas_pool_get_array(meshPool, 'vertCoordMovementWeights', & + ! vertCoordMovementWeights) + !call mpas_pool_get_array(meshPool, 'meshScalingDel2', & + ! meshScalingDel2) + !call mpas_pool_get_array(meshPool, 'meshScalingDel4', & + ! meshScalingDel4) + !call mpas_pool_get_array(meshPool, 'meshDensity', & + ! meshDensity) + !call mpas_pool_get_array(meshPool, 'angleEdge', & + ! angleEdge) + !call mpas_pool_get_array(meshPool, 'distanceToCoast', & + ! distanceToCoast) + ! + !call mpas_pool_get_array(meshPool, 'weightsOnEdge', & + ! weightsOnEdge) + !call mpas_pool_get_array(meshPool, 'kiteAreasOnVertex', & + ! kiteAreasOnVertex) + !call mpas_pool_get_array(meshPool, 'edgeTangentVectors', & + ! edgeTangentVectors) + !call mpas_pool_get_array(meshPool, 'edgeNormalVectors', & + ! edgeNormalVectors) + !call mpas_pool_get_array(meshPool, 'localVerticalUnitVectors', & + ! localVerticalUnitVectors) + !call mpas_pool_get_array(meshPool, 'cellTangentPlane', & + ! cellTangentPlane) + !call mpas_pool_get_array(meshPool, 'coeffs_reconstruct', & + ! coeffs_reconstruct) + + ! For masks, we converted to real masks, so need to recompute for + ! updated values. Probably only the edgeSign masks need updating, but + ! do them all to be sure. + call mpas_pool_get_array(meshPool, 'edgeSignOnCell', & + edgeSignOnCell) + call mpas_pool_get_array(meshPool, 'edgeSignOnVertex', & + edgeSignOnVertex) + +! do n = 1, nCells+1 +! do k = 1, nVertLevels +! boundaryCell(k, n) = real(boundaryCellTmp(k, n), RKIND) +! end do +! end do +! +! ! Allocatable array was updated in ocn_init_routines_setup_sign_and_index_fields. Now the +! ! pointer needs to be udpated +! do n = 1, nCells+1 +! do k = 1, maxEdges +! edgeSignOnCellTmp(k, n) = int(edgeSignOnCell(k, n)) +! end do +! end do +! +! do n = 1, nEdges+1 +! do k = 1, nVertLevels +! boundaryEdge(k, n) = real(boundaryEdgeTmp(k, n), RKIND) +! end do +! end do +! +! do n = 1, nVertices+1 +! do k = 1, nVertLevels +! boundaryVertex(k, n) = real(boundaryVertexTmp(k, n), RKIND) +! end do +! end do +! +! ! Allocatable array was updated in ocn_init_routines_setup_sign_and_index_fields. Now the +! ! pointer needs to be udpated +! do n = 1, nVertices+1 +! do k = 1, vertexDegree +! edgeSignOnVertexTmp(k, n) = & +! int(edgeSignOnVertex(k, n)) +! end do +! end do + + ! Go ahead and update all fields on device to be safe. + ! NOTE: if we end up computing some fields on the device during + ! init, the update must go the opposite direction (update host) + + !$acc update device(onSphere, & + !$acc sphereRadius, & + !$acc nCells, & + !$acc nEdges, & + !$acc nVertices, & + !$acc nCellsSolve, & + !$acc nEdgesSolve, & + !$acc nVerticesSolve, & + !$acc maxEdges, & + !$acc maxEdges2, & + !$acc vertexDegree, & + !$acc nVertLevels, & + !$acc nVertLevelsP1, & + !$acc nEdgesHalo, & + !$acc nCellsHalo, & + !$acc nVerticesHalo, & + !$acc nEdgesOnEdge, & + !$acc nEdgesOnCell, & + !$acc minLevelCell, & + !$acc minLevelEdgeTop, & + !$acc minLevelEdgeBot, & + !$acc minLevelVertexTop, & + !$acc minLevelVertexBot, & + !$acc maxLevelCell, & + !$acc maxLevelEdgeTop, & + !$acc maxLevelEdgeBot, & + !$acc maxLevelVertexTop, & + !$acc maxLevelVertexBot, & + !$acc indexToCellID, & + !$acc indexToEdgeID, & + !$acc indexToVertexID, & + !$acc edgesOnEdge, & + !$acc cellsOnEdge, & + !$acc verticesOnEdge, & + !$acc cellsOnCell, & + !$acc edgesOnCell, & + !$acc verticesOnCell, & + !$acc cellsOnVertex, & + !$acc edgesOnVertex, & + !$acc kiteIndexOnCell, & + !$acc latCell, & + !$acc lonCell, & + !$acc xCell, & + !$acc yCell, & + !$acc zCell, & + !$acc latEdge, & + !$acc lonEdge, & + !$acc xEdge, & + !$acc yEdge, & + !$acc zEdge, & + !$acc latVertex, & + !$acc lonVertex, & + !$acc xVertex, & + !$acc yVertex, & + !$acc zVertex, & + !$acc fEdge, & + !$acc fVertex, & + !$acc fCell, & + !$acc dcEdge, & + !$acc dvEdge, & + !$acc areaCell, & + !$acc invAreaCell, & + !$acc areaTriangle, & + !$acc bottomDepth, & + !$acc refBottomDepth, & + !$acc refBottomDepthTopOfCell, & + !$acc vertCoordMovementWeights, & + !$acc meshScalingDel2, & + !$acc meshScalingDel4, & + !$acc meshDensity, & + !$acc angleEdge, & + !$acc distanceToCoast, & + !$acc boundaryEdge, & + !$acc boundaryCell, & + !$acc boundaryVertex, & + !$acc edgeSignOnCell, & + !$acc edgeSignOnVertex, & + !$acc weightsOnEdge, & + !$acc kiteAreasOnVertex, & + !$acc edgeTangentVectors, & + !$acc edgeNormalVectors, & + !$acc localVerticalUnitVectors, & + !$acc cellTangentPlane, & + !$acc coeffs_reconstruct) + +!------------------------------------------------------------------------------- + + end subroutine li_meshUpdateFields !}}} + +!*********************************************************************** +! +! routine ocn_MeshMinMaxLevel +! +!> \brief initialize max level and boundary mask variables +!> \author Doug Jacobsen, Mark Petersen, Todd Ringler +!> \date September 2011 +!> \details +!> This routine initializes max level and boundary mask variables +! +!----------------------------------------------------------------------- + + subroutine meshMinMaxLevel()!{{{ + + !----------------------------------------------------------------- + ! local variables + !----------------------------------------------------------------- + + integer :: & + iCell, iEdge, iVertex, &! loop indices for cells,edges,vertices + i, k ! loop indices for neighbors, vertical + + ! End preamble + !------------- + ! Begin code + + ! minLevelEdgeTop is the minimum (shallowest) of surrounding cells + ! if the water column is dry, set minLevelCell outside of bounds + ! to prevent dry cells from determining minimum + !do iCell = 1, nCells+1 + ! if (maxLevelCell(iCell) == 0) minLevelCell(iCell)=nVertLevels+1 + !end do + !do iEdge = 1, nEdges + ! minLevelEdgeTop(iEdge) = & + ! min( minLevelCell(cellsOnEdge(1,iEdge)), & + ! minLevelCell(cellsOnEdge(2,iEdge)) ) + !end do + !minLevelEdgeTop(nEdges+1) = 0 + + !! minLevelEdgeBot is the maximum (deepest) of surrounding cells + !! prevent dry cells from determining maximum + !do iCell = 1, nCells+1 + ! if (maxLevelCell(iCell) == 0) minLevelCell(iCell) = 1 + !end do + !do iEdge = 1, nEdges + ! minLevelEdgeBot(iEdge) = & + ! max( minLevelCell(cellsOnEdge(1,iEdge)), & + ! minLevelCell(cellsOnEdge(2,iEdge)) ) + !end do + !minLevelEdgeBot(nEdges+1) = 0 + + !! minLevelVertexBot is the maximum (deepest) of surrounding cells + !do iVertex = 1,nVertices + ! minLevelVertexBot(iVertex) = & + ! minLevelCell(cellsOnVertex(1,iVertex)) + ! do i = 2, vertexDegree + ! minLevelVertexBot(iVertex) = & + ! max( minLevelVertexBot(iVertex), & + ! minLevelCell(cellsOnVertex(i,iVertex))) + ! end do + !end do + !minLevelVertexBot(nVertices+1) = 0 + + !! minLevelVertexTop is the minimum (shallowest) of surrounding cells + !! if the water column is dry, set minLevelCell outside of bounds + !! to prevent dry cells from determining minimum + !do iCell = 1, nCells+1 + ! if (maxLevelCell(iCell) == 0) minLevelCell(iCell) = nVertLevels+1 + !end do + !do iVertex = 1,nVertices + ! minLevelVertexTop(iVertex) = minLevelCell(cellsOnVertex(1,iVertex)) + ! do i = 2, vertexDegree + ! minLevelVertexTop(iVertex) = & + ! min( minLevelVertexTop(iVertex), & + ! minLevelCell(cellsOnVertex(i,iVertex))) + ! end do + !end do + !minLevelVertexTop(nVertices+1) = 0 + + !! if the water column is dry, set minLevelCell = 1 + !do iCell = 1, nCells+1 + ! if (maxLevelCell(iCell) == 0) minLevelCell(iCell) = 1 + !end do + ! + !! maxLevelEdgeTop is the minimum (shallowest) of surrounding cells + !do iEdge = 1, nEdges + ! maxLevelEdgeTop(iEdge) = & + ! min( maxLevelCell(cellsOnEdge(1,iEdge)), & + ! maxLevelCell(cellsOnEdge(2,iEdge)) ) + !end do + !maxLevelEdgeTop(nEdges+1) = 0 + + !! maxLevelEdgeBot is the maximum (deepest) of surrounding cells + !do iEdge = 1, nEdges + ! maxLevelEdgeBot(iEdge) = & + ! max( maxLevelCell(cellsOnEdge(1,iEdge)), & + ! maxLevelCell(cellsOnEdge(2,iEdge)) ) + !end do + !maxLevelEdgeBot(nEdges+1) = 0 + + !! maxLevelVertexBot is the maximum (deepest) of surrounding cells + !do iVertex = 1,nVertices + ! maxLevelVertexBot(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex)) + ! do i = 2, vertexDegree + ! maxLevelVertexBot(iVertex) = & + ! max( maxLevelVertexBot(iVertex), & + ! maxLevelCell(cellsOnVertex(i,iVertex))) + ! end do + !end do + !maxLevelVertexBot(nVertices+1) = 0 + + !! maxLevelVertexTop is the minimum (shallowest) of surrounding cells + !do iVertex = 1,nVertices + ! maxLevelVertexTop(iVertex) = maxLevelCell(cellsOnVertex(1,iVertex)) + ! do i = 2, vertexDegree + ! maxLevelVertexTop(iVertex) = & + ! min( maxLevelVertexTop(iVertex), & + ! maxLevelCell(cellsOnVertex(i,iVertex))) + ! end do + !end do + !maxLevelVertexTop(nVertices+1) = 0 + + ! set boundary edge + boundaryEdge(1:nEdges+1)=0.0_RKIND + !edgeMask(1:nEdges+1)=0 + + ! TODO: Replace this loop over levels with MALI levels + !do iEdge = 1, nEdges + !do k= minLevelEdgeBot(iEdge),maxLevelEdgeTop(iEdge) + ! boundaryEdge(k,iEdge)=0.0_RKIND + !end do + !end do + + ! Find cells and vertices that have an edge on the boundary + boundaryCell(1:nCells+1) = 0.0_RKIND + boundaryVertex(1:nVertices+1) = 0.0_RKIND + + do iEdge = 1, nEdges + do k = 1, nVertLevels + if (boundaryEdge(iEdge) == 1) then + boundaryCell(cellsOnEdge(1,iEdge)) = 1.0_RKIND + boundaryCell(cellsOnEdge(2,iEdge)) = 1.0_RKIND + boundaryVertex(verticesOnEdge(1,iEdge)) = 1.0_RKIND + boundaryVertex(verticesOnEdge(2,iEdge)) = 1.0_RKIND + endif + end do + end do + + !do iCell = 1, nCells + !do k = 1, nVertLevels + ! if ( minLevelCell(iCell) <= k .and. & + ! maxLevelCell(iCell) >= k ) then + ! end if + !end do + !end do + + !do iVertex = 1, nVertices + !do k = 1, nVertLevels + ! if ( minLevelVertexTop(iVertex) <= k .and. & + ! maxLevelVertexBot(iVertex) >= k ) then + ! end if + !end do + !end do + + ! Note: We do not update halos on maxLevel* variables. We want + ! the outside edge of a halo to be zero on each processor. + + !-------------------------------------------------------------------- + + end subroutine meshMinMaxLevel !}}} + +!*********************************************************************** +! +! routine li_meshSignIndexFields +! +!> \brief set up sign and index fields +!> \author Doug Jacobsen, Mark Petersen, Todd Ringler +!> \date September 2011 +!> \details +!> This routine initializes edgeSignOnCell, edgeSignOnVertex, and +!> kiteIndexOnCell. +! +!----------------------------------------------------------------------- + + subroutine meshSignIndexFields()!{{{ + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + integer :: iCell, iEdge, iVertex, i, j + + ! End preamble + !------------- + ! Begin code + + ! Initialize to zero + edgeSignOnCell = 0.0_RKIND + edgeSignOnVertex = 0.0_RKIND + + do iCell = 1, nCells + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + iVertex = verticesOnCell(i, iCell) + + ! Vector points from cell 1 to cell 2 + if (iCell == cellsOnEdge(1, iEdge)) then + edgeSignOnCell(i, iCell) = -1.0_RKIND + else + edgeSignOnCell(i, iCell) = 1.0_RKIND + end if + + end do + end do + + do iVertex = 1, nVertices + do i = 1, vertexDegree + iEdge = edgesOnVertex(i, iVertex) + + ! Vector points from vertex 1 to vertex 2 + if (iVertex == verticesOnEdge(1,iEdge)) then + edgeSignOnVertex(i,iVertex) = -1.0_RKIND + else + edgeSignOnVertex(i,iVertex) = 1.0_RKIND + end if + end do + end do + + !-------------------------------------------------------------------- + + end subroutine meshSignIndexFields !}}} + +!*********************************************************************** +! +! routine li_meshVertCoordInit +! +!> \brief initialize vertical coordinate variables +!> \author Doug Jacobsen, Mark Petersen, Todd Ringler +!> \date September 2011 +!> \details +!> This routine initializes vertical coordinate variables, including +!> zlevel-type variables and adjusted initial conditions for partial +!> bottom cells. +! +!----------------------------------------------------------------------- + +! subroutine meshVertCoordInit(domain)!{{{ +! +! !----------------------------------------------------------------- +! ! Input/output variables +! !----------------------------------------------------------------- +! +! type (domain_type), intent(inout) :: domain !< [inout] ocean state +! +! !----------------------------------------------------------------- +! ! Local variables +! !----------------------------------------------------------------- +! +! type (block_type), pointer :: block ! all ocn data in a block +! +! type (mpas_pool_type), pointer :: & +! verticalMeshPool, &! vertical mesh data +! statePool, &! ocean state (layer thickness) +! forcingPool ! forcing data (land ice mask) +! +! integer :: & +! iCell, k, &! loop indices for cell, vertical loops +! kmin, kmax ! min, max active levels in column +! +! logical :: & +! consistentSSH, &! flag if SSH consistent with depth +! landIceFlag, &! flag if land ice mask exists +! checkSSH ! flag for checking SSH in cell column +! +! real (kind=RKIND) & +! thickDepth ! depth based on sum of layer thicknesses +! +! integer, dimension(:), pointer :: & +! landIceMask ! land ice mask +! +! real (kind=RKIND), dimension(:), pointer :: & +! refZMid, &! reference depth at mid-layer +! refLayerThickness ! reference profile layer thickness +! +! real (kind=RKIND), dimension(:,:), pointer :: & +! layerThickness ! current layer thickness +! +! ! End preamble +! !------------- +! ! Begin code +! +! ! Extract pools from domain +! block => domain % blocklist +! call mpas_pool_get_subpool(block%structs, 'state', statePool) +! call mpas_pool_get_subpool(block%structs, 'verticalMesh', & +! verticalMeshPool) +! call mpas_pool_get_subpool(block%structs, 'forcing', forcingPool) +! +! ! Extract needed variables from pools +! call mpas_pool_get_array(statePool, 'layerThickness', & +! layerThickness, 1) +! +! call mpas_pool_get_array(verticalMeshPool, 'refZMid', refZMid) +! call mpas_pool_get_array(verticalMeshPool, 'refLayerThickness', & +! refLayerThickness) +! +! call mpas_pool_get_array(forcingPool, 'landIceMask', landIceMask) +! if (associated(landIceMask)) then +! landIceFlag = .true. ! ice shelves +! else +! landIceFlag = .false. ! ice shelves +! endif +! +! ! Initialize reference z coordinates based on refBottomDepth +! refBottomDepthTopOfCell(1) = 0.0_RKIND +! do k = 1, nVertLevels +! refBottomDepthTopOfCell(k+1) = refBottomDepth(k) +! refLayerThickness(k) = refBottomDepth(k) - & +! refBottomDepthTopOfCell(k) +! refZMid(k) = - refBottomDepthTopOfCell(k) - & +! refLayerThickness(k)/2.0_RKIND +! end do +! +! ! Initialization of vertCoordMovementWeights. +! ! This determines how SSH perturbations are distributed throughout +! ! the column. Check for invalid choice of vert coord movement. +! +! call mpas_log_write(' Vertical coordinate movement is: ' // & +! trim(config_vert_coord_movement)) +! +! +! select case (trim(config_vert_coord_movement)) +! case ('fixed') +! +! vertCoordMovementWeights = 0.0_RKIND +! vertCoordMovementWeights(1) = 1.0_RKIND +! +! case ('uniform_stretching') +! +! vertCoordMovementWeights = 1.0_RKIND +! +! case ('tapered') +! +! ! Set weight tapering: +! ! 1.0 shallower than config_vert_taper_weight_depth_1 +! ! linear in between +! ! 0.0 deeper than config_vert_taper_weight_depth_2 +! do k = 1, nVertLevels +! vertCoordMovementWeights(k) = 1.0_RKIND + & +! (refZMid(k) + config_vert_taper_weight_depth_1) / & +! (config_vert_taper_weight_depth_2 - & +! config_vert_taper_weight_depth_1) +! ! limit range to (0,1) +! vertCoordMovementWeights(k) = max( 0.0_RKIND, & +! min( 1.0_RKIND, vertCoordMovementWeights(k) ) ) +! end do +! +! case ('impermeable_interfaces') +! +! ! weights are never used, no init needed +! +! case ('user_specified') +! +! ! weights should be input? probably should check +! +! case default +! +! call mpas_log_write( & +! 'Incorrect choice of config_vert_coord_movement.', & +! MPAS_LOG_CRIT) +! +! end select +! +! !----------------------------------------------------------------- +! ! Check consistency and validity of options +! !----------------------------------------------------------------- +! +! ! SSH consistency +! if (config_check_ssh_consistency) then +! +! ! Check if abs(ssh)>20m. If so, print warning. +! consistentSSH = .true. +! do iCell = 1,nCells +! +! ! Only check for open ocean when land ice is being used +! checkSSH = .true. +! if (landIceFlag) then +! if (landIceMask(iCell) /= 0) checkSSH = .false. +! endif +! +! ! Check whether thickness is consistent with bottom depth +! if (checkSSH) then +! kmin = minLevelCell(iCell) +! kmax = maxLevelCell(iCell) +! +! thickDepth = sum(layerThickness(kmin:kmax,iCell)) +! +! if (abs(thickDepth - bottomDepth(iCell)) > 20.0_RKIND) then +! consistentSSH = .false. +! +! call mpas_log_write( & +! ' Warning: Sea surface height outside of acceptable range (20m)', & +! MPAS_LOG_ERR) +! call mpas_log_write( & +! ' iCell: $i, minLevelCell(iCell), maxLevelCell(iCell): $i, $i, bottomDepth(iCell): $r, sum(h): $r', & +! intArgs=(/iCell, minLevelCell(iCell), maxLevelCell(iCell) /), & +! realArgs=(/ bottomDepth(iCell), thickDepth /) ) +! endif ! depth diff check +! endif ! checkSSH +! end do ! cell loop +! +! if (.not. consistentSSH) call mpas_log_write( & +! 'Warning: SSH not consistent - ' // & +! 'initial layerThickness does not match bottomDepth.') +! +! endif ! config_check_ssh_consistency +! +! ! Z-level consistency +! if (config_check_zlevel_consistency) then +! do iCell = 1,nCells +! ! Check that bottomDepth and maxLevelCell match. Some +! ! older meshs do not have the bottomDepth variable. +! kmax = maxLevelCell(iCell) +! if (bottomDepth(iCell) > refBottomDepth(kmax) .or. & +! bottomDepth(iCell) < refBottomDepthTopOfCell(kmax)) then +! call mpas_log_write( & +! ' fatal error: bottomDepth and maxLevelCell do not match:', & +! MPAS_LOG_ERR) +! call mpas_log_write( & +! ' iCell: $i, maxLevelCell(iCell): $i, bottomDepth(iCell): $r', & +! intArgs=(/iCell, kmax /), & +! realArgs=(/ bottomDepth(iCell) /) ) +! endif +! enddo ! cell loop +! endif ! check_zlevel_consistency +! +! ! pressure gradient compatibility +! if (config_vert_coord_movement /= 'impermeable_interfaces' .and. & +! config_pressure_gradient_type == 'MontgomeryPotential') then +! call mpas_log_write( 'Incompatible choice of ' // & +! 'impermeable vertical interfaces and Montgomery Potential',& +! MPAS_LOG_CRIT) +! end if +! +! ! barotropic filter mode compatibility +! if (config_filter_btr_mode .and. & +! config_vert_coord_movement /= 'fixed')then +! call mpas_log_write( & +! 'filter_btr_mode has only been tested with ' // & +! 'config_vert_coord_movement=fixed.', & +! MPAS_LOG_CRIT) +! endif +! +! !-------------------------------------------------------------------- +! +! end subroutine meshVertCoordInit !}}} + +!*********************************************************************** +! +! routine li_meshAreaWeights +! +!> \brief set up area weighting +!> \author Mark Petersen +!> \date June 2020 +!> \details +!> This routine initializes edgeAreaFractionOfCell which is defined as +!> the fractional area of cell iCell encompassed by the triangle with +!> edge iEdge connected to the cell center of cell iCell. On a perfect +!> planar hex mesh this is always 1/6. This weighting is used to +!> interpolate scalars from edges to cell centers. +!> Use 2*edgeAreaFractionOfCell to interpolate normal vectors at +!> edges to vector norms at cell centers. +! +!----------------------------------------------------------------------- + + subroutine meshAreaWeights()!{{{ + + !----------------------------------------------------------------- + ! Local variables + !----------------------------------------------------------------- + + integer :: & + iCell, iEdge, &! cell, edge addresses + i, &! edge on cell index + ierr ! internal error flag + + ! End preamble + !------------- + ! Begin code + + ! Allocate space + allocate (edgeAreaFractionOfCell(maxEdges,nCells), STAT=ierr) +#ifdef MPAS_OPENACC + !$acc enter data create(edgeAreaFractionOfCell) +#endif + if (ierr /= 0) call mpas_log_write( & + 'Error allocating edgeAreaFractionOfCell in li_mesh', & + MPAS_LOG_CRIT) + + do iCell = 1, nCells + do i = 1, nEdgesOnCell(iCell) + iEdge = edgesOnCell(i, iCell) + edgeAreaFractionOfCell(i, iCell) = 0.25_RKIND* & + dcEdge(iEdge)*dvEdge(iEdge)/areaCell(iCell) + end do + end do +#ifdef MPAS_OPENACC + !$acc update device(edgeAreaFractionOfCell) +#endif + + !-------------------------------------------------------------------- + + end subroutine meshAreaWeights !}}} + +!*********************************************************************** +! +! routine li_meshScaling +! +!> \brief Computes mesh scaling variables for mixing +!> \author Doug Jacobsen, Mark Petersen, Todd Ringler +!> \date September 2011 +!> \details +!> This routine initializes meshScalingDel2 and meshScalingDel4 +!> that are used to scale coefficients with mesh size. +! +!----------------------------------------------------------------------- + +! subroutine meshScaling() !{{{ +! +! !----------------------------------------------------------------- +! ! local variables +! !----------------------------------------------------------------- +! +! integer :: & +! iEdge, &! loop index/address for edges +! cell1, cell2 ! cell addresses for neighbors across edge +! +! real (kind=RKIND) :: & +! cellWidth, &! reference cell width for scaling +! avgDensity ! avg mesh density across edge +! +! ! End preamble +! !------------- +! ! Begin code +! +! if (config_hmix_scaleWithMesh) then +! if (config_hmix_use_ref_cell_width) then +! ! Mesh scaling is set by areaCell and and an input +! ! reference widht config_hmix_ref_cell_width +! ! See description of config_hmix_ref_cell_width in +! ! Registry.xml for more detail. +! do iEdge = 1, nEdges +! cell1 = cellsOnEdge(1,iEdge) +! cell2 = cellsOnEdge(2,iEdge) +! ! Effective cell width at edge iEdge, assuming +! ! neighboring cells are circles for this calculation. +! cellWidth = 2.0_RKIND* & +! sqrt((areaCell(cell1) + areaCell(cell2))/ & +! 2.0_RKIND/pii) +! meshScalingDel2(iEdge) = cellWidth/ & +! config_hmix_ref_cell_width +! meshScalingDel4(iEdge) = (cellWidth/ & +! config_hmix_ref_cell_width)**3 +! end do +! +! else ! not using ref cell width +! ! Mesh scaling is set by meshDensity. This is both confusing +! ! and inconvenient, as the flags like config_mom_del2 need +! ! to be reset for every resolution. It is kept for backwards +! ! compatibility, but should become defunct. +! ! del2 scales as dc**1, del4 scales as dc**3 +! do iEdge = 1, nEdges +! cell1 = cellsOnEdge(1,iEdge) +! cell2 = cellsOnEdge(2,iEdge) +! avgDensity = (meshDensity(cell1) + meshDensity(cell2))/ & +! 2.0_RKIND +! meshScalingDel2(iEdge) = 1.0_RKIND/ & +! (avgDensity/config_maxMeshDensity)**(1.0_RKIND/4.0_RKIND) +! meshScalingDel4(iEdge) = 1.0_RKIND/ & +! (avgDensity/config_maxMeshDensity)**(3.0_RKIND/4.0_RKIND) +! end do +! end if +! +! else ! no scaling with mesh +! +! ! If config_hmix_scaleWithMesh is false, hmix coefficients +! ! do not vary with cell size but remain constant across domain. +! do iEdge = 1, nEdges +! meshScalingDel2(iEdge) = 1.0_RKIND +! meshScalingDel4(iEdge) = 1.0_RKIND +! end do +! end if +! +! !-------------------------------------------------------------------- +! +! end subroutine meshScaling !}}} + +!*********************************************************************** + +end module li_mesh + +!||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| +! vim: foldmethod=marker diff --git a/components/mpas-albany-landice/src/shared/mpas_li_setup.F b/components/mpas-albany-landice/src/shared/mpas_li_setup.F index dd6b2f2431c5..e231b4733bbd 100644 --- a/components/mpas-albany-landice/src/shared/mpas_li_setup.F +++ b/components/mpas-albany-landice/src/shared/mpas_li_setup.F @@ -121,7 +121,7 @@ subroutine li_setup_config_options( domain, err ) call mpas_pool_get_config(liConfigs, 'config_tracer_advection', config_tracer_advection) if (( (trim(config_thermal_solver) == 'temperature') .or. & (trim(config_thermal_solver) == 'enthalpy') ) .and. & - (trim(config_tracer_advection) /= 'fo') )then + ( (trim(config_tracer_advection) /= 'fo') .and. (trim(config_tracer_advection) /= 'fct') ) )then call mpas_log_write("Setting config_tracer_advection='fo' because a thermal solver has been selected. " // & "(config_tracer_advection was set to: " // trim(config_tracer_advection) // ")", MPAS_LOG_WARN) config_tracer_advection = 'fo'