<p><b>xylar@lanl.gov</b> 2012-09-03 08:00:40 -0600 (Mon, 03 Sep 2012)</p><p>* BRANCH COMMIT *<br>
<br>
Added a new module ocn_ice_shelf_fluxes to compute heat, salt and fresh<br>
water fluxes under ice shelves and to add these to the temperature and<br>
salinity tendancies. 6 new variables have been added to the registry:<br>
iceShelfCellMask: a mask that is 1 where a cell is under an ice shelf, 0 elsewhere<br>
iceShelfEdgeMask: same but on edges (edges are considered to be under ice shelves only if both neighboring cells are)<br>
iceShelfFreezingTemperature: the local freezing temperature<br>
iceShelfHeatFlux: the positive upward heat flux (divided by rho Cp)<br>
iceShelfFreshWaterFlux: the positive upward fresh-water flux (m/s)<br>
iceShelfSaltFlux: the positive upward salt flux (divided by rho)<br>
<br>
The fluxes are computed at the end of ocn_diagnostic_solve; they are<br>
added to the tracer tendencies at the end of ocn_tend_scalar.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/Makefile
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/Makefile 2012-09-03 13:43:30 UTC (rev 2138)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/Makefile 2012-09-03 14:00:40 UTC (rev 2139)
@@ -36,6 +36,7 @@
mpas_ocn_vmix_coefs_rich.o \
mpas_ocn_vmix_coefs_tanh.o \
mpas_ocn_restoring.o \
+ mpas_ocn_ice_shelf_fluxes.o \
mpas_ocn_tendency.o \
mpas_ocn_tracer_advection.o \
mpas_ocn_tracer_advection_std.o \
@@ -153,6 +154,8 @@
mpas_ocn_restoring.o:
+mpas_ocn_ice_shelf_fluxes.o:
+
mpas_ocn_vmix.o: mpas_ocn_vmix_coefs_const.o mpas_ocn_vmix_coefs_rich.o mpas_ocn_vmix_coefs_tanh.o
mpas_ocn_vmix_coefs_const.o:
@@ -204,6 +207,7 @@
mpas_ocn_vmix_coefs_rich.o \
mpas_ocn_vmix_coefs_tanh.o \
mpas_ocn_restoring.o \
+ mpas_ocn_ice_shelf_fluxes.o \
mpas_ocn_tracer_advection.o \
mpas_ocn_tracer_advection_std.o \
mpas_ocn_tracer_advection_std_hadv.o \
Modified: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/Registry
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/Registry 2012-09-03 13:43:30 UTC (rev 2138)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/Registry 2012-09-03 14:00:40 UTC (rev 2139)
@@ -87,6 +87,7 @@
namelist real restore config_restoreT_timescale 90.0
namelist real restore config_restoreS_timescale 90.0
namelist logical restore config_use_monthly_forcing false
+namelist logical ice_shelf_fluxes config_iceShelfFluxesOn true
%
% dim type name_in_file name_in_code
@@ -274,7 +275,6 @@
var persistent real uSrcReconstructMeridional ( nVertLevels nCells Time ) 2 o uSrcReconstructMeridional state - -
var persistent real MontPot ( nVertLevels nCells Time ) 2 - MontPot state - -
var persistent real pressure ( nVertLevels nCells Time ) 2 - pressure state - -
-var persistent real surfacePressure ( nCells ) 1 iro surfacePressure mesh - -
var persistent real wTop ( nVertLevelsP1 nCells Time ) 2 - wTop state - -
var persistent real rhoDisplaced ( nVertLevels nCells Time ) 2 - rhoDisplaced state - -
@@ -308,3 +308,13 @@
var persistent real acc_uReconstructMeridionalVar ( nVertLevels nCells Time ) 2 o acc_uReconstructMeridionalVar state - -
var persistent real acc_u ( nVertLevels nEdges Time ) 2 o acc_u state - -
var persistent real acc_uVar ( nVertLevels nEdges Time ) 2 o acc_uVar state - -
+
+% Fields related to ice shelves
+var persistent real surfacePressure ( nCells ) 0 iro surfacePressure mesh - -
+var persistent integer iceShelfCellMask ( nCells ) 0 iro iceShelfCellMask mesh - -
+var persistent integer iceShelfEdgeMask ( nEdges ) 0 iro iceShelfEdgeMask mesh - -
+var persistent real iceShelfFreezingTemperature ( nCells Time ) 2 o iceShelfFreezingTemperature state - -
+var persistent real iceShelfHeatFlux ( nCells Time ) 2 o iceShelfHeatFlux state - -
+var persistent real iceShelfFreshWaterFlux ( nCells Time ) 2 o iceShelfFreshWaterFlux state - -
+var persistent real iceShelfSaltFlux ( nCells Time ) 2 o iceShelfSaltFlux state - -
+
Added: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_ice_shelf_fluxes.F
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_ice_shelf_fluxes.F (rev 0)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_ice_shelf_fluxes.F 2012-09-03 14:00:40 UTC (rev 2139)
@@ -0,0 +1,262 @@
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+!
+! ocn_ice_shelf_fluxes
+!
+!> \brief MPAS ocean ice-shelf fluxes
+!> \author Xylar Asay-Davis
+!> \date 26 August 2012
+!> \version SVN:$Id:$
+!> \details
+!> This module contains the main driver routine for computing
+!> tendencies from heat and salt fluxes from ice shelves.
+!
+!-----------------------------------------------------------------------
+
+module ocn_ice_shelf_fluxes
+
+ use mpas_grid_types
+ use mpas_configure
+
+ implicit none
+ private
+ save
+
+ !--------------------------------------------------------------------
+ !
+ ! Public parameters
+ !
+ !--------------------------------------------------------------------
+
+ !--------------------------------------------------------------------
+ !
+ ! Public member functions
+ !
+ !--------------------------------------------------------------------
+
+ public :: ocn_ice_shelf_fluxes_tend, &
+ ocn_ice_shelf_fluxes_diagnostic, &
+ ocn_ice_shelf_fluxes_init
+
+ !--------------------------------------------------------------------
+ !
+ ! Private module variables
+ !
+ !--------------------------------------------------------------------
+
+ logical :: iceShelvesOn !< Flag to turn on/off ice-shelf fluxes
+
+
+!***********************************************************************
+
+contains
+
+!***********************************************************************
+!
+! routine ocn_ice_shelf_fluxes_tend
+!
+!> \brief Computes tendency term for ice-shelf fluxes
+!> \author Xylar Asay-Davis
+!> \date 26 August 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the tracer tendency due to ice-shelf fluxes
+!> based on current state.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_ice_shelf_fluxes_tend(grid, h, s, tend, err)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ type (state_type), intent(in) :: s !< Input: State information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:,:), intent(inout) :: &
+ tend !< Input/Output: velocity tendency
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ integer, intent(out) :: err !< Output: Error flag
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iCell, nCellsSolve, k
+
+ real (kind=RKIND), dimension(:), pointer :: heatFlux, saltFlux
+
+ integer :: indexT, indexS
+ err = 0
+
+ if(.not.iceShelvesOn) return
+ indexT = s%index_temperature
+ indexS = s%index_salinity
+ heatFlux => s % iceShelfHeatFlux % array
+ saltFlux => s % iceShelfSaltFlux % array
+
+ nCellsSolve = grid % nCellsSolve
+
+ k = 1 ! fluxes only in top layer
+ do iCell=1,nCellsSolve
+ tend(indexT, k, iCell) = tend(indexT, k, iCell) - heatFlux(iCell)
+ tend(indexS, k, iCell) = tend(indexS, k, iCell) - saltFlux(iCell)
+ enddo
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_ice_shelf_fluxes_tend!}}}
+
+!***********************************************************************
+!
+! routine ocn_ice_shelf_fluxes_diagnostic
+!
+!> \brief Computes ice-shelf fluxes as diagnostic variables
+!> \author Xylar Asay-Davis
+!> \date 2 September 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine computes the heat, salt and fresh-water fluxes, and the
+!> freezing temperature under ice shelves based on current state.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_ice_shelf_fluxes_diagnostic(grid, h, s)!{{{
+
+ !-----------------------------------------------------------------
+ !
+ ! input variables
+ !
+ !-----------------------------------------------------------------
+
+ real (kind=RKIND), dimension(:,:), intent(in) :: &
+ h !< Input: thickness
+
+ type (mesh_type), intent(in) :: &
+ grid !< Input: grid information
+
+ !-----------------------------------------------------------------
+ !
+ ! input/output variables
+ !
+ !-----------------------------------------------------------------
+
+ type (state_type), intent(inout) :: s !< Input/Output: State information
+
+ !-----------------------------------------------------------------
+ !
+ ! output variables
+ !
+ !-----------------------------------------------------------------
+
+ !-----------------------------------------------------------------
+ !
+ ! local variables
+ !
+ !-----------------------------------------------------------------
+
+ integer :: iCell, nCellsSolve, k
+
+ real (kind=RKIND), dimension(:), pointer :: freezingTemperature, heatFlux, freshWaterFlux, saltFlux
+
+ real (kind=RKIND), parameter :: &
+ lambda_1 = -0.0589, & ! best fit linearization of potential freezing point from Jackett, McDougall,
+ lambda_2 = 0.153, & ! Feistel, Wright and Griffies (2006) over the range 30 < S < 40 psu, 0 < p < 200 bar
+ lambda_3 = -8.04e-8 ! theta_f = lambda_1*S + lambda_2 + lambda_3*p
+
+ real (kind=RKIND), parameter :: gamma = 1e-4 ! exchange velocity in m/s
+ real (kind=RKIND), parameter :: referenceSalinity = 34.4 !psu
+ real (kind=RKIND), parameter :: surfaceReferenceSalinity = 30.0 !psu
+ real (kind=RKIND), parameter :: cp_oceanWater = 3974 !J kg^−1 K^−1
+ real (kind=RKIND), parameter :: latentHeatOfFusion = 334000 !J kg^−1
+
+ integer, dimension(:), pointer :: iceShelfCellMask
+
+ real (kind=RKIND), dimension(:,:,:), pointer :: &
+ tracers
+
+ real (kind=RKIND), dimension(:,:), pointer :: &
+ pressure
+
+ integer :: indexT, indexS
+
+ if(.not.iceShelvesOn) return
+
+ pressure => s % pressure % array
+ tracers => s % tracers % array
+ indexT = s%index_temperature
+ indexS = s%index_salinity
+ freezingTemperature => s % iceShelfFreezingTemperature % array
+ heatFlux => s % iceShelfHeatFlux % array
+ freshWaterFlux => s % iceShelfFreshWaterFlux % array
+ saltFlux => s % iceShelfSaltFlux % array
+
+ nCellsSolve = grid % nCellsSolve
+ iceShelfCellMask => grid % iceShelfCellMask % array
+
+ k = 1 ! fluxes only in top layer
+ do iCell=1,nCellsSolve
+ freezingTemperature(iCell) = lambda_1*tracers(indexS, k, iCell) + lambda_2 + lambda_3*pressure(k, iCell)
+ heatFlux(iCell) = gamma*(tracers(indexT, k, iCell)-freezingTemperature(iCell))
+ freshWaterFlux(iCell) = heatFlux(iCell)*cp_oceanWater/latentHeatOfFusion
+ !saltFlux(iCell) = freshWaterFlux(iCell) * (iceShelfCellMask(iCell)*referenceSalinity &
+ ! + (1 - iceShelfCellMask(iCell))*surfaceReferenceSalinity)
+ saltFlux(iCell) = freshWaterFlux(iCell) * referenceSalinity
+ enddo
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_ice_shelf_fluxes_diagnostic!}}}
+
+!***********************************************************************
+!
+! routine ocn_ice_shelf_fluxes_init
+!
+!> \brief Initializes ocean tracer ice-shelf fluxes
+!> \author Xylar Asay-Davis
+!> \date 26 August 2012
+!> \version SVN:$Id$
+!> \details
+!> This routine in itializes a variety of quantities related to
+!> rice-shelf fluxes.
+!
+!-----------------------------------------------------------------------
+
+ subroutine ocn_ice_shelf_fluxes_init(err)!{{{
+
+ integer, intent(out) :: err !< Output: error flag
+
+ err = 0
+ iceShelvesOn = config_iceShelfFluxesOn
+
+ !--------------------------------------------------------------------
+
+ end subroutine ocn_ice_shelf_fluxes_init!}}}
+
+!***********************************************************************
+
+end module ocn_ice_shelf_fluxes
+
+!|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
+! vim: foldmethod=marker
Modified: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_mpas_core.F
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_mpas_core.F 2012-09-03 13:43:30 UTC (rev 2138)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_mpas_core.F 2012-09-03 14:00:40 UTC (rev 2139)
@@ -23,6 +23,7 @@
use ocn_tracer_hmix
use ocn_gm
use ocn_restoring
+ use ocn_ice_shelf_fluxes
use ocn_equation_of_state
@@ -83,6 +84,9 @@
call ocn_restoring_init(err_tmp)
err = ior(err, err_tmp)
+ call ocn_ice_shelf_fluxes_init(err_tmp)
+ err = ior(err, err_tmp)
+
call ocn_vmix_init(err_tmp)
err = ior(err, err_tmp)
Modified: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_tendency.F
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_tendency.F 2012-09-03 13:43:30 UTC (rev 2138)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_tendency.F 2012-09-03 14:00:40 UTC (rev 2139)
@@ -37,6 +37,7 @@
use ocn_tracer_vadv
use ocn_tracer_hmix
use ocn_restoring
+ use ocn_ice_shelf_fluxes
use ocn_equation_of_state
use ocn_vmix
@@ -51,6 +52,7 @@
type (timer_node), pointer :: thickHadvTimer, thickVadvTimer
type (timer_node), pointer :: velCorTimer, velVadvTimer, velPgradTimer, velHmixTimer, velForceTimer, velExpVmixTimer
type (timer_node), pointer :: tracerHadvTimer, tracerVadvTimer, tracerHmixTimer, tracerExpVmixTimer, tracerRestoringTimer
+ type (timer_node), pointer :: iceShelfFluxesTimer
!--------------------------------------------------------------------
!
@@ -284,7 +286,7 @@
real (kind=RKIND), intent(in) :: dt !< Input: Time step
real (kind=RKIND), dimension(:,:), pointer :: &
- uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh
+ uTransport, h,wTop, h_edge, vertDiffTopOfCell, tend_h, uh, pressure
real (kind=RKIND), dimension(:,:,:), pointer :: &
tracers, tend_tr
@@ -298,6 +300,7 @@
tracers => s % tracers % array
h_edge => s % h_edge % array
vertDiffTopOfCell => d % vertDiffTopOfCell % array
+ pressure => s % pressure % array
tend_tr => tend % tracers % array
tend_h => tend % h % array
@@ -369,6 +372,12 @@
call mpas_timer_stop("restoring", tracerRestoringTimer)
+ call mpas_timer_start("ice_shelf_fluxes", .false., iceShelfFluxesTimer)
+
+ call ocn_ice_shelf_fluxes_tend(grid, h, s, tend_tr, err)
+
+ call mpas_timer_stop("ice_shelf_fluxes", iceShelfFluxesTimer)
+
10 format(2i8,10e20.10)
call mpas_timer_stop("ocn_tend_scalar")
@@ -799,10 +808,8 @@
do iCell=1,nCells
! pressure for generalized coordinates
- ! assume atmospheric pressure at the surface is zero for now.
- ! mrp 120724 use same rho as top cell right now.
- ! Later this would be a rho that represents the land ice.
+ ! atmospheric or ice pressure given by surfacePressure (currently time independent)
pressure(1,iCell) = surfacePressure(iCell) + rho(1,iCell)*gravity*0.5*h(1,iCell)
do k=2,maxLevelCell(iCell)
@@ -853,6 +860,8 @@
uBolusGM = 0.0
end if
+ call ocn_ice_shelf_fluxes_diagnostic(grid, h, s)
+
end subroutine ocn_diagnostic_solve!}}}
!***********************************************************************
Modified: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2012-09-03 13:43:30 UTC (rev 2138)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_vel_forcing_bottomdrag.F 2012-09-03 14:00:40 UTC (rev 2139)
@@ -110,6 +110,7 @@
integer :: iEdge, nEdgesSolve, k
integer, dimension(:), pointer :: maxLevelEdgeTop
integer, dimension(:,:), pointer :: edgeMask
+ integer, dimension(:), pointer :: iceShelfEdgeMask
!-----------------------------------------------------------------
!
@@ -126,6 +127,7 @@
nEdgesSolve = grid % nEdgesSolve
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
edgeMask => grid % edgeMask % array
+ iceShelfEdgeMask => grid % iceShelfEdgeMask % array
do iEdge=1,grid % nEdgesSolve
@@ -137,6 +139,10 @@
tend(k,iEdge) = tend(k,iEdge)-edgeMask(k,iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
+ ! top drag under ice shelves
+ k = 1
+ tend(k,iEdge) = tend(k,iEdge)-iceShelfEdgeMask(iEdge)*(bottomDragCoef*u(k,iEdge)*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge))
+
enddo
Modified: branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_vmix.F
===================================================================
--- branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_vmix.F 2012-09-03 13:43:30 UTC (rev 2138)
+++ branches/ocean_projects/surface_pressure_ice_shelves/src/core_ocean/mpas_ocn_vmix.F 2012-09-03 14:00:40 UTC (rev 2139)
@@ -295,6 +295,8 @@
integer, dimension(:,:), pointer :: cellsOnEdge
+ integer, dimension(:), pointer :: iceShelfEdgeMask
+
real (kind=RKIND), dimension(:), allocatable :: A, B, C, uTemp
err = 0
@@ -305,6 +307,7 @@
nVertLevels = grid % nVertLevels
maxLevelEdgeTop => grid % maxLevelEdgeTop % array
cellsOnEdge => grid % cellsOnEdge % array
+ iceShelfEdgeMask => grid % iceShelfEdgeMask % array
allocate(A(nVertLevels),B(nVertLevels),C(nVertLevels),uTemp(nVertLevels))
A(1)=0
@@ -343,9 +346,15 @@
enddo
! Apply bottom drag boundary condition on the viscous term
- B(N) = 1 - A(N) + dt*config_bottom_drag_coeff &
+ k = N
+ B(k) = 1 - A(k) + dt*config_bottom_drag_coeff &
*sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+ ! Apply top drag boundary condition on the viscous term
+ k = 1
+ B(k) = B(k) + iceShelfEdgeMask(iEdge)*dt*config_bottom_drag_coeff &
+ *sqrt(2.0*ke_edge(k,iEdge))/h_edge(k,iEdge)
+
call tridiagonal_solve(A(2:N),B,C(1:N-1),u(:,iEdge),uTemp,N)
u(1:N,iEdge) = uTemp(1:N)
</font>
</pre>