<p><b>akt@lanl.gov</b> 2013-04-19 12:44:55 -0600 (Fri, 19 Apr 2013)</p><p>Initial versions of velocity solver. Divergence of stress calculated with both weak integral formulation and variational discretization on both triangles and arbitary polygons<br>
</p><hr noshade><pre><font color="gray">Modified: branches/cice_projects/initial_cice_core/src/core_cice/Makefile
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/Makefile 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/Makefile 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,35 +1,29 @@
.SUFFIXES: .F .o
-OBJS = mpas_cice_kinds_mod.o \
- mpas_cice_constants.o \
- mpas_cice_domain_size.o \
- mpas_cice_state.o \
- mpas_cice_test.o \
- mpas_cice_init.o \
- mpas_cice_timestep.o \
- mpas_cice_mpas_core.o
+OBJS = mpas_cice_mpas_core.o \
+ mpas_cice_time_integration.o \
+ mpas_cice_dynamics_evp_tri.o \
+ mpas_cice_dynamics_evp_hex.o \
+ mpas_cice_dynamics_shared.o \
+ mpas_cice_testing.o
all: core_cice
core_cice: $(OBJS)
ar -ru libdycore.a $(OBJS)
-mpas_cice_kinds_mod.o:
+mpas_cice_testing.o:
-mpas_cice_constants.o: mpas_cice_kinds_mod.o
+mpas_cice_dynamics_shared.o: mpas_cice_testing.o
-mpas_cice_domain_size.o: mpas_cice_kinds_mod.o
+mpas_cice_dynamics_evp_tri.o: mpas_cice_testing.o mpas_cice_dynamics_shared.o
-mpas_cice_state.o: mpas_cice_kinds_mod.o mpas_cice_domain_size.o
+mpas_cice_dynamics_evp_hex.o: mpas_cice_testing.o mpas_cice_dynamics_shared.o
-mpas_cice_init.o: mpas_cice_state.o
+mpas_cice_time_integration.o:
-mpas_cice_test.o: mpas_cice_kinds_mod.o mpas_cice_constants.o
+mpas_cice_mpas_core.o: mpas_cice_dynamics_evp_hex.o mpas_cice_dynamics_evp_tri.o mpas_cice_time_integration.o
-mpas_cice_timestep.o: mpas_cice_kinds_mod.o mpas_cice_test.o
-
-mpas_cice_mpas_core.o: mpas_cice_test.o mpas_cice_init.o mpas_cice_timestep.o
-
clean:
$(RM) *.o *.mod *.f90 libdycore.a
Modified: branches/cice_projects/initial_cice_core/src/core_cice/Registry
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/Registry 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/Registry 2013-04-19 18:44:55 UTC (rev 2780)
@@ -36,14 +36,12 @@
dim vertexDegree vertexDegree
dim nVertLevels nVertLevels
dim nTracers nTracers
-dim nCategories nCategories
-dim nIceLayers nIceLayers
-dim nSnowLayers nSnowLayers
%
% var persistence type name_in_file ( dims ) time_levs iro- name_in_code struct super-array array_class
%
var persistent text xtime ( Time ) 2 ro xtime state - -
+var persistent integer somevar ( TWO ) 2 - somevar state - -
var persistent real latCell ( nCells ) 0 iro latCell mesh - -
var persistent real lonCell ( nCells ) 0 iro lonCell mesh - -
@@ -96,6 +94,7 @@
var persistent real fEdge ( nEdges ) 0 iro fEdge mesh - -
var persistent real fVertex ( nVertices ) 0 iro fVertex mesh - -
var persistent real fCell ( nCells ) 0 iro fCell mesh - -
+var persistent real h_s ( nCells ) 0 iro h_s mesh - -
% Space needed for advection
var persistent real deriv_two ( FIFTEEN TWO nEdges ) 0 o deriv_two mesh - -
@@ -104,37 +103,148 @@
% !! NOTE: the following arrays are needed to allow the use
% !! of the module_advection.F w/o alteration
% Space needed for deformation calculation weights
-%var persistent real defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
-%var persistent real defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
-%var persistent real kdiff ( nVertLevels nCells Time ) 0 - kdiff mesh - -
+var persistent real defc_a ( maxEdges nCells ) 0 - defc_a mesh - -
+var persistent real defc_b ( maxEdges nCells ) 0 - defc_b mesh - -
+var persistent real kdiff ( nVertLevels nCells Time ) 0 - kdiff mesh - -
% Arrays required for reconstruction of velocity field
var persistent real coeffs_reconstruct ( R3 maxEdges nCells ) 0 - coeffs_reconstruct mesh - -
% Boundary conditions: read from input, saved in restart and written to output
-%var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
-%var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
-%var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
-%var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
+var persistent integer boundaryEdge ( nVertLevels nEdges ) 0 iro boundaryEdge mesh - -
+var persistent integer boundaryVertex ( nVertLevels nVertices ) 0 iro boundaryVertex mesh - -
+var persistent integer boundaryCell ( nVertLevels nCells ) 0 iro boundaryCell mesh - -
+var persistent real u_src ( nVertLevels nEdges ) 0 iro u_src mesh - -
-%--------------------------------------------------------------------------------------
+%-------------------------------------------------------------------
+% ice state
+var persistent real iceAreaCell ( nCells ) 0 - iceAreaCell icestate - -
+var persistent real iceAreaVertex ( nVertices ) 0 - iceAreaVertex icestate - -
+var persistent real iceVolumeCell ( nCells ) 0 - iceVolumeCell icestate - -
+var persistent real totalMassCell ( nCells ) 0 - totalMassCell icestate - -
+var persistent real totalMassVertex ( nVertices ) 0 - totalMassVertex icestate - -
-% state
-var persistent real iceAreaCell ( nCells ) 0 - iceAreaCell states - -
-var persistent real iceVolumeCell ( nCells ) 0 - iceVolumeCell states - -
-var persistent real snowVolumeCell ( nCells ) 0 - snowVolumeCell states - -
-var persistent real iceAreaCategory ( nCells nCategories ) 0 - iceAreaCategory states - -
-var persistent real iceVolumeCategory ( nCells nCategories ) 0 - iceVolumeCategory states - -
-var persistent real snowVolumeCategory ( nCells nCategories ) 0 - snowVolumeCategory states - -
+% ice dynamics
+var persistent real uVelocity ( nCells ) 0 - uVelocity dynamics - -
+var persistent real vVelocity ( nCells ) 0 - vVelocity dynamics - -
+var persistent real icePressure ( nCells ) 0 - icePressure dynamics - -
+var persistent real stress1 ( vertexDegree nVertices ) 0 - stress1 dynamics - -
+var persistent real stress2 ( vertexDegree nVertices ) 0 - stress2 dynamics - -
+var persistent real stress12 ( vertexDegree nVertices ) 0 - stress12 dynamics - -
+var persistent real stressDivergenceU ( nCells ) 0 - stressDivergenceU dynamics - -
+var persistent real stressDivergenceV ( nCells ) 0 - stressDivergenceV dynamics - -
+var persistent real airStressU ( nCells ) 0 - airStressU dynamics - -
+var persistent real airStressV ( nCells ) 0 - airStressV dynamics - -
+var persistent real oceanStressU ( nCells ) 0 - oceanStressU dynamics - -
+var persistent real oceanStressV ( nCells ) 0 - oceanStressV dynamics - -
+var persistent real oceanStressCoeff ( nCells ) 0 - oceanStressCoeff dynamics - -
+var persistent real surfaceTiltForceU ( nCells ) 0 - surfaceTiltForceU dynamics - -
+var persistent real surfaceTiltForceV ( nCells ) 0 - surfaceTiltForceV dynamics - -
-% tracers
-var persistent real iceEnthalpy ( nIceLayers nCategories nCells ) 0 - iceEnthalpy tracer iceVolumeLayerTracers tracers
-var persistent real iceSalinity ( nIceLayers nCategories nCells ) 0 - iceSalinity tracer iceVolumeLayerTracers tracers
-var persistent real snowEnthalpy ( nSnowLayers nCategories nCells ) 0 - snowEnthalpy tracer snowVolumeLayerTracers wibble
-var persistent real snowSalinity ( nSnowLayers nCategories nCells ) 0 - snowSalinity tracer snowVolumeLayerTracers wibble
+% basis
+var persistent real basisGradientU ( vertexDegree nVertices ) 0 - basisGradientU basis - -
+var persistent real basisGradientV ( vertexDegree nVertices ) 0 - basisGradientV basis - -
+var persistent real basisIntegrals ( vertexDegree nVertices ) 0 - basisIntegrals basis - -
+var persistent integer triangleCornerAtCellCenter ( maxEdges nCells ) 0 - triangleCornerAtCellCenter basis - -
+var persistent integer edgeOppositeTriangleCorner ( vertexDegree nVertices ) 0 - edgeOppositeTriangleCorner basis - -
-% flux
-var persistent real seaSurfaceTemperature ( nCells ) 0 - seaSurfaceTemperature flux - -
-var persistent real seaFreezingTemperature ( nCells ) 0 - seaFreezingTemperature flux - -
-var persistent real iceLateralMeltFraction ( nCells ) 0 - iceLateralMeltFraction flux - -
-var persistent real seaFreezeMeltPotential ( nCells ) 0 - seaFreezeMeltPotential flux - -
+% forcing
+var persistent real uAirVelocity ( nCells ) 0 - uAirVelocity forcing - -
+var persistent real vAirVelocity ( nCells ) 0 - vAirVelocity forcing - -
+var persistent real uOceanVelocity ( nCells ) 0 - uOceanVelocity forcing - -
+var persistent real vOceanVelocity ( nCells ) 0 - vOceanVelocity forcing - -
+
+% boundary
+var persistent integer boundaryCell ( nCells ) 0 - boundaryCell boundary - -
+var persistent integer boundaryCell2 ( nCells ) 0 - boundaryCell2 boundary - -
+var persistent integer boundaryEdge ( nEdges ) 0 - boundaryEdge boundary - -
+var persistent integer boundaryVertex ( nVertices ) 0 - boundaryVertex boundary - -
+var persistent integer boundaryVertex2 ( nVertices ) 0 - boundaryVertex2 boundary - -
+var persistent integer interiorVertex ( nVertices ) 0 - interiorVertex boundary - -
+var persistent integer interiorVertex2 ( nVertices ) 0 - interiorVertex2 boundary - -
+
+%-------------------------------------------------------------------
+% hex basis
+%var persistent real xLocal ( maxEdges nCells ) 0 - xLocal hexbas - -
+%var persistent real yLocal ( maxEdges nCells ) 0 - yLocal hexbas - -
+%var persistent real wachspressA ( maxEdges nCells ) 0 - wachspressA hexbas - -
+%var persistent real wachspressB ( maxEdges nCells ) 0 - wachspressB hexbas - -
+%var persistent real wachspressKappa ( maxEdges maxEdges nCells ) 0 - wachspressKappa hexbas - -
+%var persistent real basisIntegralsU ( maxEdges maxEdges nCells ) 0 - basisIntegralsU hexbas - -
+%var persistent real basisIntegralsV ( maxEdges maxEdges nCells ) 0 - basisIntegralsV hexbas - -
+%var persistent real basisGradientU ( maxEdges maxEdges nCells ) 0 - basisGradientU hexbas - -
+%var persistent real basisGradientV ( maxEdges maxEdges nCells ) 0 - basisGradientV hexbas - -
+%var persistent integer cellVerticesAtVertex ( vertexDegree nVertices ) 0 - cellVerticesAtVertex hexbas - -
+
+var persistent real xLocal ( vertexDegree nVertices ) 0 - xLocal hexbas - -
+var persistent real yLocal ( vertexDegree nVertices ) 0 - yLocal hexbas - -
+var persistent real wachspressA ( vertexDegree nVertices ) 0 - wachspressA hexbas - -
+var persistent real wachspressB ( vertexDegree nVertices ) 0 - wachspressB hexbas - -
+var persistent real wachspressKappa ( vertexDegree vertexDegree nVertices ) 0 - wachspressKappa hexbas - -
+var persistent real basisIntegralsU ( vertexDegree vertexDegree nVertices ) 0 - basisIntegralsU hexbas - -
+var persistent real basisIntegralsV ( vertexDegree vertexDegree nVertices ) 0 - basisIntegralsV hexbas - -
+var persistent real basisGradientU ( vertexDegree vertexDegree nVertices ) 0 - basisGradientU hexbas - -
+var persistent real basisGradientV ( vertexDegree vertexDegree nVertices ) 0 - basisGradientV hexbas - -
+var persistent integer cellVerticesAtVertex ( maxEdges nCells ) 0 - cellVerticesAtVertex hexbas - -
+
+
+% hex dynamics
+var persistent real uVelocity ( nVertices ) 0 - uVelocity hexdyn - -
+var persistent real vVelocity ( nVertices ) 0 - vVelocity hexdyn - -
+var persistent real icePressure ( nCells ) 0 - icePressure hexdyn - -
+%var persistent real stress1 ( maxEdges nCells ) 0 - stress1 hexdyn - -
+%var persistent real stress2 ( maxEdges nCells ) 0 - stress2 hexdyn - -
+%var persistent real stress12 ( maxEdges nCells ) 0 - stress12 hexdyn - -
+%var persistent real stressDivergenceU ( nVertices ) 0 - stressDivergenceU hexdyn - -
+%var persistent real stressDivergenceV ( nVertices ) 0 - stressDivergenceV hexdyn - -
+var persistent real airStressU ( nVertices ) 0 - airStressU hexdyn - -
+var persistent real airStressV ( nVertices ) 0 - airStressV hexdyn - -
+var persistent real oceanStressU ( nVertices ) 0 - oceanStressU hexdyn - -
+var persistent real oceanStressV ( nVertices ) 0 - oceanStressV hexdyn - -
+var persistent real oceanStressCoeff ( nVertices ) 0 - oceanStressCoeff hexdyn - -
+var persistent real surfaceTiltForceU ( nVertices ) 0 - surfaceTiltForceU hexdyn - -
+var persistent real surfaceTiltForceV ( nVertices ) 0 - surfaceTiltForceV hexdyn - -
+
+var persistent real stress1 ( vertexDegree nVertices ) 0 - stress1 hexdyn - -
+var persistent real stress2 ( vertexDegree nVertices ) 0 - stress2 hexdyn - -
+var persistent real stress12 ( vertexDegree nVertices ) 0 - stress12 hexdyn - -
+var persistent real stressDivergenceU ( nCells ) 0 - stressDivergenceU hexdyn - -
+var persistent real stressDivergenceV ( nCells ) 0 - stressDivergenceV hexdyn - -
+
+% hex forcing
+var persistent real uAirVelocity ( nVertices ) 0 - uAirVelocity hexfor - -
+var persistent real vAirVelocity ( nVertices ) 0 - vAirVelocity hexfor - -
+var persistent real uOceanVelocity ( nVertices ) 0 - uOceanVelocity hexfor - -
+var persistent real vOceanVelocity ( nVertices ) 0 - vOceanVelocity hexfor - -
+
+%-------------------------------------------------------------------
+% weak integral formulation
+var persistent real uVelocity ( nVertices ) 0 - uVelocity weak - -
+var persistent real vVelocity ( nVertices ) 0 - vVelocity weak - -
+var persistent real icePressure ( nCells ) 0 - icePressure weak - -
+var persistent real strain11 ( nCells ) 0 - strain11 weak - -
+var persistent real strain22 ( nCells ) 0 - strain22 weak - -
+var persistent real strain12 ( nCells ) 0 - strain12 weak - -
+var persistent real stress11 ( nCells ) 0 - stress11 weak - -
+var persistent real stress22 ( nCells ) 0 - stress22 weak - -
+var persistent real stress12 ( nCells ) 0 - stress12 weak - -
+var persistent real principalStress1 ( nCells ) 0 - principalStress1 weak - -
+var persistent real principalStress2 ( nCells ) 0 - principalStress2 weak - -
+var persistent real stressDivergenceU ( nVertices ) 0 - stressDivergenceU weak - -
+var persistent real stressDivergenceV ( nVertices ) 0 - stressDivergenceV weak - -
+var persistent real normalVectorPolygon ( TWO maxEdges nCells ) 0 - normalVectorPolygon weak - -
+var persistent real normalVectorTriangle ( TWO vertexDegree nVertices ) 0 - normalVectorTriangle weak - -
+
+
+var persistent real airStressU ( nVertices ) 0 - airStressU weak - -
+var persistent real airStressV ( nVertices ) 0 - airStressV weak - -
+var persistent real oceanStressU ( nVertices ) 0 - oceanStressU weak - -
+var persistent real oceanStressV ( nVertices ) 0 - oceanStressV weak - -
+var persistent real oceanStressCoeff ( nVertices ) 0 - oceanStressCoeff weak - -
+var persistent real surfaceTiltForceU ( nVertices ) 0 - surfaceTiltForceU weak - -
+var persistent real surfaceTiltForceV ( nVertices ) 0 - surfaceTiltForceV weak - -
+
+var persistent real uAirVelocity ( nVertices ) 0 - uAirVelocity weak - -
+var persistent real vAirVelocity ( nVertices ) 0 - vAirVelocity weak - -
+var persistent real uOceanVelocity ( nVertices ) 0 - uOceanVelocity weak - -
+var persistent real vOceanVelocity ( nVertices ) 0 - vOceanVelocity weak - -
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_constants.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_constants.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_constants.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,242 +0,0 @@
-!=======================================================================
-!BOP
-!
-! !MODULE: ice_constants - sets physical constants
-!
-! !DESCRIPTION:
-!
-! This module defines a variety of physical and numerical constants
-! used throughout the ice model \\
-!
-! Code originally based on constants.F in POP
-!
-! !REVISION HISTORY:
-! SVN:$Id: ice_constants.F90 567 2013-01-07 02:57:36Z eclare $
-!
-! author Elizabeth C. Hunke, LANL
-!
-! !INTERFACE:
-
- module cice_constants
-!
-! !USES:
-!
- use cice_kinds_mod
-!
-!EOP
-!
- implicit none
- !save
-
- !-----------------------------------------------------------------
- ! physical constants
- !-----------------------------------------------------------------
-
-#ifdef AOMIP
- real (kind=dbl_kind), parameter :: &
- rhos = 300.0_dbl_kind ,&! density of snow (kg/m^3)
- rhoi = 900.0_dbl_kind ,&! density of ice (kg/m^3)
- rhow = 1025.0_dbl_kind ,&! density of seawater (kg/m^3)
- cp_air = 1004.0_dbl_kind ,&! specific heat of air (J/kg/K)
- emissivity= 0.98_dbl_kind ,&! emissivity of snow and ice
- cp_ice = 2090._dbl_kind ,&! specific heat of fresh ice (J/kg/K)
- cp_brine = 4650._dbl_kind ,&! specific heat of brine (J/kg/K)
- ! from Taylor and Feltham (2004)
- cp_ocn = 4190._dbl_kind ,&! specific heat of ocn (J/kg/K)
- depressT = 0.0575_dbl_kind ,&! Tf:brine salinity ratio (C/ppt)
- dragio = 0.0055_dbl_kind ,&! ice-ocn drag coefficient
- albocn = 0.10_dbl_kind ! ocean albedo
-#else
-! CICE default parameters
- real (kind=dbl_kind), parameter :: &
- rhos = 330.0_dbl_kind ,&! density of snow (kg/m^3)
- rhoi = 917.0_dbl_kind ,&! density of ice (kg/m^3)
- rhow = 1026.0_dbl_kind ,&! density of seawater (kg/m^3)
- cp_air = 1005.0_dbl_kind ,&! specific heat of air (J/kg/K)
- ! (Briegleb JGR 97 11475-11485 July 1992)
- emissivity= 0.95_dbl_kind ,&! emissivity of snow and ice
- cp_ice = 2106._dbl_kind ,&! specific heat of fresh ice (J/kg/K)
- cp_ocn = 4218._dbl_kind ,&! specific heat of ocn (J/kg/K)
- ! freshwater value needed for enthalpy
- depressT = 0.054_dbl_kind ,&! Tf:brine salinity ratio (C/ppt)
- dragio = 0.00536_dbl_kind ,&! ice-ocn drag coefficient
- albocn = 0.06_dbl_kind ! ocean albedo
-#endif
-
-! UNESCO melting temperature parameters
- real (kind=dbl_kind), parameter :: &
- mlt_a = 0.0575_dbl_kind ,& !nonlinear melting T coefficient (oC/ppt)
- mlt_b = 1.710523e-3_dbl_kind ,& !(oC/ppt^(3/2))
- mlt_c = 2.154996e-4_dbl_kind !(oC/ppt^2)
-
-
- real (kind=dbl_kind), parameter :: &
- gravit = 9.80616_dbl_kind ,&! gravitational acceleration (m/s^2)
- omega = 7.292e-5_dbl_kind ,&! angular velocity of earth (rad/sec)
- radius = 6.37e6_dbl_kind ! earth radius (m)
-
- real (kind=dbl_kind), parameter :: &
- pi = 3.14159265358979323846_dbl_kind,&! pi
- secday = 86400.0_dbl_kind ,&! seconds in calendar day
- Tocnfrz = -1.8_dbl_kind ,&! freezing temp of seawater (C),
- ! used as Tsfcn for open water
- rhofresh = 1000.0_dbl_kind ,&! density of fresh water (kg/m^3)
- zvir = 0.606_dbl_kind ,&! rh2o/rair - 1.0
- vonkar = 0.4_dbl_kind ,&! von Karman constant
- cp_wv = 1.81e3_dbl_kind ,&! specific heat of water vapor (J/kg/K)
- stefan_boltzmann = 567.0e-10_dbl_kind,&! W/m^2/K^4
- Tffresh = 273.15_dbl_kind ,&! freezing temp of fresh ice (K)
- Lsub = 2.835e6_dbl_kind ,&! latent heat, sublimation freshwater (J/kg)
- Lvap = 2.501e6_dbl_kind ,&! latent heat, vaporization freshwater (J/kg)
- Lfresh = Lsub-Lvap ,&! latent heat of melting of fresh ice (J/kg)
- Timelt = 0.0_dbl_kind ,&! melting temperature, ice top surface (C)
- Tsmelt = 0.0_dbl_kind ,&! melting temperature, snow top surface (C)
- ice_ref_salinity = 4._dbl_kind ,&! (ppt)
-! ocn_ref_salinity = 34.7_dbl_kind,&! (ppt)
-! rho_air = 1.2_dbl_kind ,&! ambient air density (kg/m^3)
- spval_dbl = 1.0e30_dbl_kind ! special value (double precision)
-
- real (kind=real_kind), parameter :: &
- spval = 1.0e30_real_kind ! special value for netCDF output
-
- real (kind=dbl_kind), parameter :: &
- iceruf = 0.0005_dbl_kind ,&! ice surface roughness (m)
-
- ! (Ebert, Schramm and Curry JGR 100 15965-15975 Aug 1995)
- kappav = 1.4_dbl_kind ,&! vis extnctn coef in ice, wvlngth<700nm (1/m)
- kappan = 17.6_dbl_kind,&! vis extnctn coef in ice, wvlngth<700nm (1/m)
-
- kice = 2.03_dbl_kind ,&! thermal conductivity of fresh ice(W/m/deg)
- kseaice= 2.00_dbl_kind ,&! thermal conductivity of sea ice (W/m/deg)
- ! (used in zero layer thermodynamics option)
- ksno = 0.30_dbl_kind ,&! thermal conductivity of snow (W/m/deg)
- zref = 10._dbl_kind ,&! reference height for stability (m)
- hsmin = 0.0001_dbl_kind,&! minimum allowed snow depth (m) for DE
- snowpatch = 0.02_dbl_kind,&! parameter for fractional snow area (m)
- salt_loss = 0.65_dbl_kind ! zbgc fraction of available salt retained (1.0 = no loss)
-
- ! weights for albedos
-! 4 Jan 2007 BPB Following are appropriate for complete cloud
-! in a summer polar atmosphere with 1.5m bare sea ice surface:
-! .636/.364 vis/nir with only 0.5% direct for each band.
- real (kind=dbl_kind), parameter :: & ! currently used only
- awtvdr = 0.00318_dbl_kind, &! visible, direct ! for history and
- awtidr = 0.00182_dbl_kind, &! near IR, direct ! diagnostics
- awtvdf = 0.63282_dbl_kind, &! visible, diffuse
- awtidf = 0.36218_dbl_kind ! near IR, diffuse
-
- real (kind=dbl_kind), parameter :: &
- qqqice = 11637800._dbl_kind ,&! for qsat over ice
- TTTice = 5897.8_dbl_kind ,&! for qsat over ice
- qqqocn = 627572.4_dbl_kind ,&! for qsat over ocn
- TTTocn = 5107.4_dbl_kind ! for qsat over ocn
-
- ! these are currently set so as to have no effect on the decomposition
- real (kind=dbl_kind), parameter :: &
- shlat = 30.0_dbl_kind ,&! artificial masking edge (deg)
- nhlat = -30.0_dbl_kind ! artificial masking edge (deg)
-
- !-----------------------------------------------------------------
- ! numbers
- !-----------------------------------------------------------------
-
- real (kind=dbl_kind), parameter :: &
- c0 = 0.0_dbl_kind, &
- c1 = 1.0_dbl_kind, &
- c1p5 = 1.5_dbl_kind, &
- c2 = 2.0_dbl_kind, &
- c3 = 3.0_dbl_kind, &
- c4 = 4.0_dbl_kind, &
- c5 = 5.0_dbl_kind, &
- c6 = 6.0_dbl_kind, &
- c8 = 8.0_dbl_kind, &
- c9 = 9.0_dbl_kind, &
- c10 = 10.0_dbl_kind, &
- c12 = 12.0_dbl_kind, &
- c15 = 15.0_dbl_kind, &
- c16 = 16.0_dbl_kind, &
- c20 = 20.0_dbl_kind, &
- c25 = 25.0_dbl_kind, &
- c30 = 30.0_dbl_kind, &
- c100 = 100.0_dbl_kind, &
- c180 = 180.0_dbl_kind, &
- c360 = 360.0_dbl_kind, &
- c365 = 365.0_dbl_kind, &
- c400 = 400.0_dbl_kind, &
- c3600= 3600.0_dbl_kind, &
- c1000= 1000.0_dbl_kind, &
- p001 = 0.001_dbl_kind, &
- p01 = 0.01_dbl_kind, &
- p025 = 0.025_dbl_kind, &
- p1 = 0.1_dbl_kind, &
- p2 = 0.2_dbl_kind, &
- p4 = 0.4_dbl_kind, &
- p5 = 0.5_dbl_kind, &
- p6 = 0.6_dbl_kind, &
- p05 = 0.05_dbl_kind, &
- p15 = 0.15_dbl_kind, &
- p25 = 0.25_dbl_kind, &
- p75 = 0.75_dbl_kind, &
- p166 = c1/c6, &
- p333 = c1/c3, &
- p666 = c2/c3, &
- p111 = c1/c9, &
- p055 = p111*p5, &
- p027 = p055*p5, &
- p222 = c2/c9, &
- eps11 = 1.0e-11_dbl_kind, &
- eps13 = 1.0e-13_dbl_kind, &
- eps16 = 1.0e-16_dbl_kind, &
- puny = eps11, &
- bignum = 1.0e+30_dbl_kind, &
- pih = p5*pi, &
- piq = p5*pih, &
- pi2 = c2*pi, &
- days_per_4c = 146097.0_dbl_kind, &
- days_per_c = 36524.0_dbl_kind, &
- days_per_4y = 1461.0_dbl_kind, &
- days_per_y = 365.0_dbl_kind
-
- !-----------------------------------------------------------------
- ! location of fields for staggered grids
- !-----------------------------------------------------------------
-
- integer (int_kind), parameter :: &
- field_loc_unknown = 0, &
- field_loc_noupdate = -1, &
- field_loc_center = 1, &
- field_loc_NEcorner = 2, &
- field_loc_Nface = 3, &
- field_loc_Eface = 4, &
- field_loc_Wface = 5
-
-
- !-----------------------------------------------------------------
- ! field type attribute - necessary for handling
- ! changes of direction across tripole boundary
- !-----------------------------------------------------------------
-
- integer (int_kind), parameter :: &
- field_type_unknown = 0, &
- field_type_noupdate = -1, &
- field_type_scalar = 1, &
- field_type_vector = 2, &
- field_type_angle = 3
-
- !-----------------------------------------------------------------
- ! conversion factors
- !-----------------------------------------------------------------
-
- real (kind=dbl_kind), parameter :: &
- cm_to_m = 0.01_dbl_kind ,&! cm to meters
- m_to_cm = 100._dbl_kind ,&! meters to cm
- m2_to_km2 = 1.e-6_dbl_kind ,&! m^2 to km^2
- kg_to_g = 1000._dbl_kind ,&! kilograms to grams
- mps_to_cmpdy = 8.64e6_dbl_kind ,&! m per s to cm per day
- rad_to_deg = 180._dbl_kind/pi ! degree-radian conversion
-
-!=======================================================================
-
- end module cice_constants
-
-!=======================================================================
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_domain_size.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_domain_size.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_domain_size.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,102 +0,0 @@
-!=======================================================================
-!BOP
-!
-! !MODULE: ice_domain_size
-!
-! !DESCRIPTION:
-!
-! Defines the global domain size and number of categories and layers.
-! Code originally based on domain_size.F in POP
-!
-! !REVISION HISTORY:
-! SVN:$Id: ice_domain_size.F90 574 2013-01-13 16:20:47Z eclare $
-!
-! author Elizabeth C. Hunke, LANL
-! 2004: Block structure and snow parameters added by William Lipscomb
-! Renamed (used to be ice_model_size)
-! 2006: Converted to free source form (F90) by Elizabeth Hunke
-! Removed hardwired sizes (NX...can now be set in compile scripts)
-!
-! !INTERFACE:
-!
- module cice_domain_size
-!
-! !USES:
-!
- use cice_kinds_mod
-!
-!EOP
-!=======================================================================
-
- implicit none
- save
-
- ! added instead of preprocessor
- integer (kind=int_kind), parameter :: &
- NXGLOB = 1 , &
- NYGLOB = 1 , &
- NICECAT = 1 , &
- NICELYR = 4 , &
- NSNWLYR = 1 , &
- NTRAERO = 0 , &
- NBGCLYR = 1 , &
- TRAGE = 0 , &
- TRFY = 0 , &
- TRLVL = 0 , &
- TRPND = 0 , &
- TRBRI = 0 , &
- BLCKX = 1 , &
- BLCKY = 1 , &
- MXBLCKS = 1
-
- integer (kind=int_kind), parameter :: &
- nx_global = NXGLOB , & ! i-axis size
- ny_global = NYGLOB , & ! j-axis size
- ncat = NICECAT , & ! number of categories
- nilyr = NICELYR , & ! number of ice layers per category
- nslyr = NSNWLYR , & ! number of snow layers per category
- max_aero = 6 , & ! maximum number of aerosols
- n_aero = NTRAERO , & ! number of aerosols in use
-
- nblyr = NBGCLYR , & ! number of biology layers per category
- nblyr_hist = nblyr+2 , & ! number of ice layer plus boundary points
- nltrcr = 6 , & ! number of layer bgc tracers
- ! number of biology layer tracers
- nbltrcr = 1 , & ! 4 basic bio; 1 nitrate only
-
- ! number of tracers (defined in ice_init)
- max_ntrcr = 1 & ! 1 = surface temperature
- + nilyr & ! ice salinity
- + nilyr & ! ice enthalpy
- + nslyr & ! snow enthalpy
- !!!!! optional tracers:
- + TRAGE & ! age
- + TRFY & ! first-year area
- + TRLVL*2 & ! level/deformed ice
- + TRPND*3 & ! ponds
- + n_aero*4 & ! number of aerosols * 4 aero layers
- + TRBRI & ! brine height
- + TRBRI*nltrcr*nblyr,&! zbgc (off if TRBRI=0)
- max_nstrm = 5 ! max number of history output streams
-
- integer (kind=int_kind), parameter :: &
- block_size_x = BLCKX , & ! size of block in first horiz dimension
- block_size_y = BLCKY ! size of block in second horiz dimension
-
- !*** The model will inform the user of the correct
- !*** values for the parameter below. A value higher than
- !*** necessary will not cause the code to fail, but will
- !*** allocate more memory than is necessary. A value that
- !*** is too low will cause the code to exit.
- !*** A good initial guess is found using
- !*** max_blocks = (nx_global/block_size_x)*(ny_global/block_size_y)/
- !*** num_procs
-
- integer (kind=int_kind), parameter :: &
- max_blocks = MXBLCKS ! max number of blocks per processor
-
-!=======================================================================
-
- end module cice_domain_size
-
-!=======================================================================
Added: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_evp_hex.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_evp_hex.F (rev 0)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_evp_hex.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -0,0 +1,1430 @@
+module cice_dynamics_evp_hex
+
+ use mpas_grid_types
+
+ implicit none
+
+ private
+ public :: init_hexagonal_basis, &
+ test_divergence_stress_tensor
+
+contains
+
+ !-------------------------------------------------------------
+
+ subroutine init_hexagonal_basis(mesh, hexbas, boundary)
+
+ use cice_testing, only: writeout_edgecell, writeout_edgeedgecell
+
+ use cice_dynamics_shared, only: boundary_cells, boundary_vertices
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(hexbas_type), pointer :: hexbas
+ type(boundary_type), pointer :: boundary
+
+ integer :: iCell, iVertex, jVertex
+
+ real(kind=RKIND) :: &
+ integration
+
+ call calc_local_coords(mesh, &
+ hexbas % xLocal % array, &
+ hexbas % yLocal % array)
+
+ call calc_wachspress_coefficients(mesh, &
+ hexbas % wachspressKappa % array, &
+ hexbas % wachspressA % array, &
+ hexbas % wachspressB % array, &
+ hexbas % xLocal % array, &
+ hexbas % yLocal % array)
+
+ call calculate_wachspress_derivatives(mesh, &
+ hexbas % basisGradientU % array, &
+ hexbas % basisGradientV % array, &
+ hexbas % wachspressA % array, &
+ hexbas % wachspressB % array, &
+ hexbas % wachspressKappa % array)
+
+ !call writeout_edgeedgecell(mesh, hexbas % wachspressKappa % array, "wachspressKappa.txt")
+ !call writeout_edgecell(mesh, hexbas % wachspressA % array, "wachspressA.txt")
+ !call writeout_edgecell(mesh, hexbas % wachspressB % array, "wachspressB.txt")
+
+ call integrate_wachspress(mesh, hexbas)
+
+ !call writeout_edgeedgecell(mesh, hexbas % basisIntegralsU % array, "integralsu.txt")
+ !call writeout_edgeedgecell(mesh, hexbas % basisIntegralsV % array, "integralsv.txt")
+
+ call calc_cell_vertices_at_vertex(mesh, &
+ hexbas % cellVerticesAtVertex % array)
+
+ !call boundary_cells(mesh, &
+ ! boundary % boundaryCell % array, &
+ ! boundary % boundaryCell2 % array)
+
+ !call boundary_vertices(mesh, &
+ ! boundary % boundaryVertex % array, &
+ ! boundary % interiorVertex % array, &
+ ! boundary % interiorVertex2 % array, &
+ ! boundary % boundaryCell % array)
+
+ end subroutine init_hexagonal_basis
+
+ !-------------------------------------------------------------
+
+ subroutine calc_local_coords(mesh, xLocal, yLocal)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ xLocal, &
+ yLocal
+
+ integer :: &
+ iCell, &
+ iVertex, &
+ iVertexOnCell
+
+ do iCell = 1, mesh % nCells
+
+ do iVertex = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iVertexOnCell = mesh % verticesOnCell % array(iVertex, iCell)
+
+ xLocal(iVertex,iCell) = mesh % xVertex % array(iVertexOnCell) - mesh % xCell % array(iCell)
+ yLocal(iVertex,iCell) = mesh % yVertex % array(iVertexOnCell) - mesh % yCell % array(iCell)
+
+ enddo ! iVertex
+
+ enddo ! iCell
+
+ end subroutine calc_local_coords
+
+ !-------------------------------------------------------------
+
+ subroutine calc_wachspress_coefficients(mesh, wachspressKappa, wachspressA, wachspressB, xLocal, yLocal)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:,:), intent(out) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ xLocal, &
+ yLocal
+
+ integer :: &
+ iCell, &
+ iVertex, &
+ i0, &
+ i1, &
+ i2, &
+ jVertex, &
+ nEdgesOnCell
+
+ ! loop over cells
+ do iCell = 1, mesh % nCells
+
+ ! number of edges
+ nEdgesOnCell = mesh % nEdgesOnCell % array(iCell)
+
+ ! loop over vertices
+ do iVertex = 1, nEdgesOnCell
+
+ ! end points of line segment
+ i1 = iVertex - 1
+ i2 = iVertex
+ if (i1 < 1) i1 = i1 + nEdgesOnCell
+
+ ! solve for the line segment equation
+ wachspressA(iVertex, iCell) = (yLocal(i2,iCell) - yLocal(i1,iCell)) / (xLocal(i1,iCell) * yLocal(i2,iCell) - xLocal(i2,iCell) * yLocal(i1,iCell))
+ wachspressB(iVertex, iCell) = (xLocal(i1,iCell) - xLocal(i2,iCell)) / (xLocal(i1,iCell) * yLocal(i2,iCell) - xLocal(i2,iCell) * yLocal(i1,iCell))
+
+ write(*,*) iCell, iVertex, i1, i2, wachspressA(iVertex, iCell), wachspressB(iVertex, iCell)
+
+ enddo ! iVertex
+
+ ! loop over vertices
+ do iVertex = 1, nEdgesOnCell
+
+ ! determine kappa
+ wachspressKappa(1,iVertex,iCell) = 1.0_RKIND
+
+ do jVertex = 2, nEdgesOnCell
+
+ ! previous, this and next vertex
+ i0 = jVertex - 1
+ i1 = jVertex
+ i2 = jVertex + 1
+ if (i2 > nEdgesOnCell) i2 = i2 - nEdgesOnCell
+
+ wachspressKappa(jVertex,iVertex,iCell) = wachspressKappa(jVertex-1,iVertex,iCell) * &
+ (wachspressA(i2,iCell) * (xLocal(i0,iCell) - xLocal(i1,iCell)) + wachspressB(i2,iCell) * (yLocal(i0,iCell) - yLocal(i1,iCell))) / &
+ (wachspressA(i0,iCell) * (xLocal(i1,iCell) - xLocal(i0,iCell)) + wachspressB(i0,iCell) * (yLocal(i1,iCell) - yLocal(i0,iCell)))
+
+ enddo ! jVertex
+
+ enddo ! iVertex
+
+ enddo ! iCell
+
+ end subroutine calc_wachspress_coefficients
+
+ !-------------------------------------------------------------
+
+ function wachspress_basis_function(nEdgesOnCell, iVertex, x, y, wachspressKappa, wachspressA, wachspressB) result(wachpress)
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ iVertex
+
+ real(kind=RKIND), intent(in) :: &
+ x, &
+ y
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: &
+ wachpress
+
+ real(kind=RKIND) :: &
+ numerator, &
+ denominator, &
+ numerator_ivertex
+
+ integer :: &
+ jVertex
+
+ ! sum over numerators to get denominator
+ denominator = 0.0_RKIND
+
+ do jVertex = 1, nEdgesOnCell
+
+ numerator = wachspress_numerator(nEdgesOnCell, jVertex, iVertex, x, y, wachspressKappa, wachspressA, wachspressB)
+
+ denominator = denominator + numerator
+
+ if (jvertex == iVertex) then
+ numerator_ivertex = numerator
+ endif
+
+ enddo ! jVertex
+
+ wachpress = numerator_ivertex / denominator
+
+ end function wachspress_basis_function
+
+ !-------------------------------------------------------------
+
+ function wachspress_basis_derivative(nEdgesOnCell, iVertex, x, y, wachspressKappa, wachspressA, wachspressB, iDerivativeType) result(wachspress)
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ iVertex, &
+ iDerivativeType
+
+ real(kind=RKIND), intent(in) :: &
+ x, &
+ y
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: &
+ wachspress
+
+ real(kind=RKIND) :: &
+ numerator, &
+ derivative, &
+ denominator, &
+ sum_of_derivatives, &
+ numerator_ivertex, &
+ derivative_ivertex
+
+ integer :: &
+ jVertex
+
+ ! sum over numerators to get denominator
+ denominator = 0.0_RKIND
+ sum_of_derivatives = 0.0_RKIND
+
+ do jVertex = 1, nEdgesOnCell
+
+ numerator = wachspress_numerator(nEdgesOnCell, jVertex, iVertex, x, y, wachspressKappa, wachspressA, wachspressB)
+ derivative = wachspress_numerator_derivative(nEdgesOnCell, jVertex, iVertex, x, y, wachspressKappa, wachspressA, wachspressB, iDerivativeType)
+
+ denominator = denominator + numerator
+ sum_of_derivatives = sum_of_derivatives + derivative
+
+ if (jvertex == iVertex) then
+ numerator_ivertex = numerator
+ derivative_ivertex = derivative
+ endif
+
+ enddo ! jVertex
+
+ wachspress = derivative_ivertex / denominator - (numerator_ivertex / denominator**2) * sum_of_derivatives
+
+ end function wachspress_basis_derivative
+
+ !-------------------------------------------------------------
+
+ function wachspress_numerator(nEdgesOnCell, jVertex, iVertex, x, y, wachspressKappa, wachspressA, wachspressB) result(numerator)
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ jVertex, &
+ iVertex
+
+ real(kind=RKIND), intent(in) :: &
+ x, &
+ y
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: numerator
+
+ integer :: &
+ kVertex, &
+ i1, &
+ i2
+
+ i1 = jVertex
+ i2 = jVertex + 1
+ if (i2 > nEdgesOnCell) i2 = i2 - nEdgesOnCell
+
+ numerator = 1.0_RKIND
+
+ do kVertex = 1, nEdgesOnCell
+
+ if (kVertex /= i1 .and. kVertex /= i2) then
+ numerator = numerator * wachspress_edge_equation(kVertex, x, y, wachspressA, wachspressB)
+ endif
+
+ enddo ! jVertex
+
+ numerator = numerator * wachspressKappa(jVertex,iVertex)
+
+ end function wachspress_numerator
+
+ !-------------------------------------------------------------
+
+ function wachspress_numerator_derivative(nEdgesOnCell, jVertex, iVertex, x, y, wachspressKappa, wachspressA, wachspressB, iDerivativeType) result(derivative)
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ jVertex, &
+ iVertex, &
+ iDerivativeType
+
+ real(kind=RKIND), intent(in) :: &
+ x, &
+ y
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: &
+ derivative
+
+ real(kind=RKIND) :: &
+ sum_of_products, &
+ product
+
+ integer :: &
+ kVertex, &
+ lVertex, &
+ i1, &
+ i2
+
+ i1 = jVertex
+ i2 = jVertex + 1
+ if (i2 > nEdgesOnCell) i2 = i2 - nEdgesOnCell
+
+ sum_of_products = 0.0_RKIND
+
+ do kVertex = 1, nEdgesOnCell
+
+ if (kVertex /= i1 .and. kVertex /= i2) then
+
+ product = 1.0_RKIND
+
+ do lVertex = 1, nEdgesOnCell
+
+ if (lVertex /= i1 .and. lVertex /= i2) then
+
+ if (lVertex == kVertex) then
+
+ product = product * wachspress_edge_equation_derivative(lVertex, wachspressA, wachspressB, iDerivativeType)
+
+ else
+
+ product = product * wachspress_edge_equation(lVertex, x, y, wachspressA, wachspressB)
+
+ endif
+
+ endif
+
+ enddo ! jVertex
+
+ sum_of_products = sum_of_products + product
+
+ endif
+
+ enddo ! jVertex
+
+ derivative = sum_of_products * wachspressKappa(jVertex,iVertex)
+
+ end function wachspress_numerator_derivative
+
+ !-------------------------------------------------------------
+
+ function wachspress_edge_equation(iVertex, x, y, wachspressA, wachspressB) result(edge_equation)
+
+ integer, intent(in) :: &
+ iVertex
+
+ real(kind=RKIND), intent(in) :: &
+ x, &
+ y
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: edge_equation
+
+ edge_equation = 1.0_RKIND - wachspressA(iVertex) * x - wachspressB(iVertex) * y
+
+ end function wachspress_edge_equation
+
+ !-------------------------------------------------------------
+
+ function wachspress_edge_equation_derivative(iVertex, wachspressA, wachspressB, iDerivativeType) result(derivative)
+
+ integer, intent(in) :: &
+ iVertex, &
+ iDerivativeType
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: derivative
+
+ if (iDerivativeType == 1) then
+
+ derivative = -wachspressA(iVertex)
+
+ else if (iDerivativeType == 2) then
+
+ derivative = -wachspressB(iVertex)
+
+ endif
+
+ end function wachspress_edge_equation_derivative
+
+ !-------------------------------------------------------------
+
+ subroutine calculate_wachspress_derivatives(mesh, &
+ basisGradientU, basisGradientV, &
+ wachspressA, wachspressB, &
+ wachspressKappa)
+
+ ! basisGradientUV(jVertexOnCell,iVertexOnCell,iCell)
+ ! iCell : The cell the gradients are based in
+ ! iVertexOnCell : The vertex basis function the gradient is calculated from
+ ! jVertexOnCell : The vertex location the gradients are calculated at
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:,:), intent(out) :: &
+ basisGradientU, &
+ basisGradientV
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND), dimension(:,:,:), intent(in) :: &
+ wachspressKappa
+
+ integer :: &
+ nEdgesOnCell
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ jVertexOnCell, &
+ jVertex
+
+ ! loop over cells
+ do iCell = 1, mesh % nCells
+
+ nEdgesOnCell = mesh % nEdgesOnCell % array(iCell)
+
+ ! loop over vertices - basis function
+ do iVertexOnCell = 1, nEdgesOnCell
+
+ ! loop over vertices again - derivative position
+ do jVertexOnCell = 1, nEdgesOnCell
+
+ jVertex = mesh % verticesOnCell % array(jVertexOnCell,iCell)
+
+ basisGradientU(jVertexOnCell,iVertexOnCell,iCell) = &
+ wachspress_basis_derivative(nEdgesOnCell, &
+ iVertexOnCell, &
+ mesh % xVertex % array(jVertex), &
+ mesh % yVertex % array(jVertex), &
+ wachspressKappa(:,:,iCell), &
+ wachspressA(:,iCell), &
+ wachspressB(:,iCell), &
+ 1)
+
+ basisGradientV(jVertexOnCell,iVertexOnCell,iCell) = &
+ wachspress_basis_derivative(nEdgesOnCell, &
+ iVertexOnCell, &
+ mesh % xVertex % array(jVertex), &
+ mesh % yVertex % array(jVertex), &
+ wachspressKappa(:,:,iCell), &
+ wachspressA(:,iCell), &
+ wachspressB(:,iCell), &
+ 2)
+
+ enddo ! jVertexOnCell
+
+ enddo ! iVertexOnCell
+
+ enddo ! iCell
+
+ end subroutine calculate_wachspress_derivatives
+
+ !-------------------------------------------------------------
+
+ subroutine get_triangle_mapping(mapping, x1, y1, x2, y2, u1, v1, u2, v2)
+
+ real(kind=RKIND), dimension(2,2), intent(out) :: &
+ mapping
+
+ real(kind=RKIND), intent(in) :: &
+ x1, & ! input
+ y1, &
+ x2, &
+ y2, &
+ u1, & ! output
+ v1, &
+ u2, &
+ v2
+
+ mapping(1,1) = (u2*y1 - u1*y2) / (x2*y1 - x1*y2)
+ mapping(1,2) = (u1*x2 - u2*x1) / (y1*x2 - y2*x1)
+
+ mapping(2,1) = (v2*y1 - v1*y2) / (x2*y1 - x1*y2)
+ mapping(2,2) = (v1*x2 - v2*x1) / (y1*x2 - y2*x1)
+
+ end subroutine get_triangle_mapping
+
+ !-------------------------------------------------------------
+
+ subroutine use_triangle_mapping(u, v, x, y, mapping)
+
+ real(kind=RKIND), intent(out) :: &
+ u, v
+
+ real(kind=RKIND), intent(in) :: &
+ x, y
+
+ real(kind=RKIND), dimension(2,2), intent(in) :: &
+ mapping
+
+ u = mapping(1,1) * x + mapping(1,2) * y
+ v = mapping(2,1) * x + mapping(2,2) * y
+
+ end subroutine use_triangle_mapping
+
+ !-------------------------------------------------------------
+
+ subroutine integrate_wachspress(mesh, hexbas)
+
+ ! basisIntegralsUV (iVertexOnCell,jVertexOnCell,iCell)
+ ! iCell : cell integrals are performed on
+ ! iVertexOnCell : vertex number of Wachspress function
+ ! jVertexOnCell : vertex number of Wachspress derivative function
+ ! Sij
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(hexbas_type), pointer :: hexbas
+
+ real(kind=RKIND) :: &
+ integration
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ jVertexOnCell, &
+ iDerivativeType, &
+ nEdgesOnCell
+
+ do iCell = 1, mesh % nCells
+
+ nEdgesOnCell = mesh % nEdgesOnCell % array(iCell)
+
+ do iVertexOnCell = 1, nEdgesOnCell
+
+ do jVertexOnCell = 1, nEdgesOnCell
+
+ do iDerivativeType = 1, 2
+
+ call integrate_wachspress_polygon(&
+ integration, nEdgesOnCell, iVertexOnCell, jvertexOnCell, &
+ hexbas % xLocal % array(:,iCell), &
+ hexbas % yLocal % array(:,iCell), &
+ hexbas % wachspressKappa % array(:,:,iCell), &
+ hexbas % wachspressA % array(:,iCell), &
+ hexbas % wachspressB % array(:,iCell), &
+ iDerivativeType)
+
+ if (iDerivativeType == 1) then
+
+ hexbas % basisIntegralsU % array(iVertexOnCell,jVertexOnCell,iCell) = &
+ integration
+
+ else if (iDerivativeType == 2) then
+
+ hexbas % basisIntegralsV % array(iVertexOnCell,jVertexOnCell,iCell) = &
+ integration
+
+ endif
+
+ enddo ! iDerivativeType
+
+ enddo ! jVertex
+
+ enddo ! iVertex
+
+ enddo ! iCell
+
+ end subroutine integrate_wachspress
+
+ !-------------------------------------------------------------
+
+ subroutine integrate_wachspress_polygon(integration, nEdgesOnCell, iVertexOnCell, jvertexOnCell, xLocal, yLocal, &
+ wachspressKappa, wachspressA, wachspressB, iDerivativeType)
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ iVertexOnCell, &
+ jvertexOnCell, &
+ iDerivativeType
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ xLocal, &
+ yLocal
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND) :: &
+ integration, &
+ integration_subtriangle
+
+ real(kind=RKIND), dimension(2,2) :: &
+ mapping
+
+ integer :: &
+ iSubTriangle, &
+ i1, &
+ i2
+
+ integration = 0.0_RKIND
+
+ do iSubTriangle = 1, nEdgesOnCell
+
+ i1 = iSubTriangle
+ i2 = iSubTriangle + 1
+ if (i2 > nEdgesOnCell) i2 = i2 - nEdgesOnCell
+
+ call get_triangle_mapping(mapping, &
+ 1.0_RKIND, 0.0_RKIND, &
+ 0.0_RKIND, 1.0_RKIND, &
+ xLocal(i1), yLocal(i1), &
+ xLocal(i2), yLocal(i2))
+
+ call integrate_wachspress_subtriangle(integration_subtriangle, nEdgesOnCell, iVertexOnCell, jvertexOnCell, &
+ wachspressKappa(:,:), wachspressA(:), wachspressB(:), iDerivativeType, mapping)
+
+ integration = integration + integration_subtriangle
+
+ enddo ! iSubTriangle
+
+ end subroutine integrate_wachspress_polygon
+
+ !-------------------------------------------------------------
+
+ subroutine integrate_wachspress_subtriangle(integration,nEdgesOnCell, iVertexOnCell, jvertexOnCell, &
+ wachspressKappa, wachspressA, wachspressB, iDerivativeType, mapping)
+
+ real(kind=RKIND), intent(out) :: &
+ integration
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ iVertexOnCell, &
+ jvertexOnCell, &
+ iDerivativeType
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND), dimension(2,2), intent(in) :: &
+ mapping
+
+ real(kind=RKIND) :: &
+ jacobian
+
+ real(kind=RKIND) :: &
+ scaling, &
+ x, &
+ y, &
+ u, &
+ v
+
+
+ integer :: &
+ i, j, k
+
+ integer, parameter :: n = 10
+
+ integration = 0.0_RKIND
+
+ do i = 0, n
+ do j = 0, n
+
+ if (i<=n-j) then
+
+ if (i==n .or. j==n .or. (i==0 .and. j==0)) then
+
+ scaling = 1.0_RKIND
+
+ else if ((j==0 .and. i/=0 .and. i/=n) .or. (i==0 .and. j/=0 .and. j/=n) .or. (i==n-j .and. i/=0 .and. j/=0)) then
+
+ scaling = 3.0_RKIND
+
+ else
+
+ scaling = 6.0_RKIND
+
+ endif
+
+ u = real(i,RKIND) / real(n,RKIND)
+ v = real(j,RKIND) / real(n,RKIND)
+
+ call use_triangle_mapping(x, y, u, v, mapping)
+
+ jacobian = mapping(1,1) * mapping(2,2) - mapping(1,2) * mapping(2,1)
+
+ ! area test
+ !integration = integration + scaling * jacobian
+
+ ! actual integration
+ integration = integration + scaling * jacobian * &
+ wachspress_basis_function(nEdgesOnCell, iVertexOnCell, x, y, wachspressKappa, wachspressA, wachspressB) * &
+ wachspress_basis_derivative(nEdgesOnCell, jVertexOnCell, x, y, wachspressKappa, wachspressA, wachspressB, iDerivativeType)
+
+ endif
+
+ enddo ! i
+ enddo ! j
+
+ integration = integration / (6.0_RKIND * real(n,RKIND)**2)
+
+ end subroutine integrate_wachspress_subtriangle
+
+ !-------------------------------------------------------------
+
+ subroutine calc_cell_vertices_at_vertex(mesh, cellVerticesAtVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:,:), intent(out) :: &
+ cellVerticesAtVertex
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell, &
+ iVertexOnCell, &
+ jVertex
+
+ do iVertex = 1, mesh % nVertices
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ cellVerticesAtVertex(iVertexDegree,iVertex) = 0
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ jVertex = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ if (iVertex == jVertex) then
+
+ cellVerticesAtVertex(iVertexDegree,iVertex) = iVertexOnCell
+
+ endif
+
+ enddo ! iVertexOnCell
+
+ enddo ! iVertexDegree
+
+ enddo ! iVertex
+
+ end subroutine calc_cell_vertices_at_vertex
+
+ !-------------------------------------------------------------
+ ! time step
+ !-------------------------------------------------------------
+
+ subroutine run_dynamics_evp(mesh, icestate, hexbas, hexdyn, hexfor, boundary)
+
+ use cice_dynamics_shared, only: ice_strength, air_stress, ocean_stress, surface_tilt, solve_velocity
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(hexdyn_type), pointer :: hexdyn
+ type(hexbas_type), pointer :: hexbas
+ type(hexfor_type), pointer :: hexfor
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND) :: timeStep, elasticTimeStep
+
+ call ice_strength(mesh, &
+ hexdyn % icePressure % array, &
+ icestate % iceAreaCell % array, &
+ icestate % iceVolumeCell % array)
+
+ call stress_tensor(mesh, &
+ hexdyn % stress1 % array, &
+ hexdyn % stress2 % array, &
+ hexdyn % stress12 % array, &
+ hexdyn % icePressure % array, &
+ hexdyn % uVelocity % array, &
+ hexdyn % vVelocity % array, &
+ hexbas % basisGradientU % array, &
+ hexbas % basisGradientV % array, &
+ timeStep, elasticTimeStep)
+
+ call divergence_stress_tensor(mesh, &
+ hexdyn % stressDivergenceU % array, &
+ hexdyn % stressDivergenceV % array, &
+ hexdyn % stress1 % array, &
+ hexdyn % stress2 % array, &
+ hexdyn % stress12 % array, &
+ hexbas % basisIntegralsU % array, &
+ hexbas % basisIntegralsV % array, &
+ hexbas % cellVerticesAtVertex % array, &
+ boundary % interiorVertex % array)
+
+ call air_stress(&
+ hexdyn % airStressU % array, &
+ hexdyn % airStressV % array, &
+ hexfor % uAirVelocity % array, &
+ hexfor % vAirVelocity % array, &
+ icestate % iceAreaVertex % array)
+
+ !call ocean_stress(&
+ ! hexdyn % oceanStressU % array, &
+ ! hexdyn % oceanStressV % array, &
+ ! hexdyn % oceanStressCoeff % array, &
+ ! hexfor % uOceanVelocity % array, &
+ ! hexfor % vOceanVelocity % array, &
+ ! hexdyn % uVelocity % array, &
+ ! hexdyn % vVelocity % array, &
+ ! icestate % iceAreaCell % array)
+
+ call surface_tilt(&
+ hexdyn % surfaceTiltForceU % array, &
+ hexdyn % surfaceTiltForceV % array)
+
+ call solve_velocity(mesh % nCells, &
+ boundary % interiorVertex % array, &
+ hexdyn % uVelocity % array, &
+ hexdyn % vVelocity % array, &
+ iceState % totalMassCell % array, &
+ mesh % fCell % array, &
+ hexdyn % stressDivergenceU % array, &
+ hexdyn % stressDivergenceV % array, &
+ hexdyn % airStressU % array, &
+ hexdyn % airStressV % array, &
+ hexdyn % surfaceTiltForceU % array, &
+ hexdyn % surfaceTiltForceV % array, &
+ hexdyn % oceanStressU % array, &
+ hexdyn % oceanStressV % array, &
+ hexdyn % oceanStressCoeff % array, &
+ timeStep)
+
+ end subroutine run_dynamics_evp
+
+ !-------------------------------------------------------------
+
+ subroutine stress_tensor(mesh, &
+ stress1, stress2, &
+ stress12, icePressure, &
+ uVelocity, vVelocity, &
+ basisGradientU, basisGradientV, &
+ timeStep, elasticTimeStep)
+
+ use cice_dynamics_shared, only: evp_constitutive_relation, linear_constitutive_relation
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ icePressure, &
+ uVelocity, &
+ vVelocity
+
+ real(kind=RKIND), dimension(:,:,:), intent(in) :: &
+ basisGradientU, &
+ basisGradientV
+
+ real(kind=RKIND), intent(in) :: &
+ timeStep, &
+ elasticTimeStep
+
+ real(kind=RKIND) :: &
+ strainDivergence, &
+ strainTension, &
+ strainShearing
+
+ integer :: &
+ iCell, &
+ jVertexOnCell, &
+ iVertexOnCell, &
+ iVertex
+
+ ! loop over cells
+ do iCell = 1, mesh % nCells
+
+ ! loop over velocity points surrounding cell - location of stress and derivative
+ do jVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ strainDivergence = 0.0_RKIND
+ strainTension = 0.0_RKIND
+ strainShearing = 0.0_RKIND
+
+ ! loop over basis functions
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ strainDivergence = strainDivergence + &
+ uVelocity(iVertex) * basisGradientU(jVertexOnCell,iVertexOnCell,iCell) + &
+ vVelocity(iVertex) * basisGradientV(jVertexOnCell,iVertexOnCell,iCell)
+
+ strainTension = strainTension + &
+ uVelocity(iVertex) * basisGradientU(jVertexOnCell,iVertexOnCell,iCell) - &
+ vVelocity(iVertex) * basisGradientV(jVertexOnCell,iVertexOnCell,iCell)
+
+ strainShearing = strainShearing + &
+ uVelocity(iVertex) * basisGradientV(jVertexOnCell,iVertexOnCell,iCell) + &
+ vVelocity(iVertex) * basisGradientU(jVertexOnCell,iVertexOnCell,iCell)
+
+ enddo ! iVertexOnCell
+
+ !call evp_constitutive_relation(stress1(jVertexOnCell,iCell), &
+ ! stress2(jVertexOnCell,iCell), &
+ ! stress12(jVertexOnCell,iCell), &
+ ! strainDivergence, strainTension, strainShearing, &
+ ! icePressure(iCell), &
+ ! timeStep, elasticTimeStep)
+
+ call linear_constitutive_relation(stress1(jVertexOnCell,iCell), &
+ stress2(jVertexOnCell,iCell), &
+ stress12(jVertexOnCell,iCell), &
+ strainDivergence, strainTension, strainShearing)
+
+
+ enddo ! jVertexOnCell
+
+ enddo ! iCell
+
+ end subroutine stress_tensor
+
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_tensor(mesh, &
+ stressDivergenceU, stressDivergenceV, &
+ stress1, stress2, &
+ stress12, &
+ basisIntegralsU, basisIntegralsV, &
+ cellVerticesAtVertex, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ stressDivergenceU, &
+ stressDivergenceV
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND), dimension(:,:,:), intent(in) :: &
+ basisIntegralsU, &
+ basisIntegralsV
+
+ integer, dimension(:,:), intent(in) :: &
+ cellVerticesAtVertex
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ real(kind=RKIND) :: &
+ stressDivergenceUCell, &
+ stressDivergenceVCell
+
+ real(kind=RKIND) :: &
+ stress1contribU, &
+ stress2contribU, &
+ stress12contribU, &
+ stress1contribV, &
+ stress2contribV, &
+ stress12contribV
+
+ integer :: &
+ iVertex, &
+ iCellOnVertex, &
+ iCell, &
+ iVertexOnCell, &
+ jVertexOnCell
+
+ ! loop over velocity positions
+ do iVertex = 1, mesh % nVertices
+
+ stressDivergenceU(iVertex) = 0.0_RKIND
+ stressDivergenceV(iVertex) = 0.0_RKIND
+
+ !if (interiorVertex(iVertex) == 1) then
+
+ stress1contribU = 0.0_RKIND
+ stress2contribU = 0.0_RKIND
+ stress12contribU = 0.0_RKIND
+
+ stress1contribV = 0.0_RKIND
+ stress2contribV = 0.0_RKIND
+ stress12contribV = 0.0_RKIND
+
+ ! loop over surrounding cells
+ do iCellOnVertex = 1, mesh % vertexDegree
+
+ ! get the cell number of this cell
+ iCell = mesh % cellsOnVertex % array(iCellOnVertex, iVertex)
+
+ ! get the vertexOnCell number of the iVertex velocity point from cell iCell
+ jVertexOnCell = cellVerticesAtVertex(iCellOnVertex,iVertex)
+
+ stressDivergenceUCell = 0.0_RKIND
+ stressDivergenceVCell = 0.0_RKIND
+
+ ! loop over the vertices of the surrounding cell
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+
+ stress1contribU = stress1contribU + stress1(iVertexOnCell,iCell) * basisIntegralsU(iVertexOnCell,jVertexOnCell,iCell)
+ stress2contribU = 0.0_RKIND
+ stress12contribU = stress12contribU + stress12(iVertexOnCell,iCell) * basisIntegralsV(iVertexOnCell,jVertexOnCell,iCell)
+
+ stress1contribV = 0.0_RKIND
+ stress2contribV = stress2contribV + stress2(iVertexOnCell,iCell) * basisIntegralsV(iVertexOnCell,jVertexOnCell,iCell)
+ stress12contribV = stress12contribV + stress12(iVertexOnCell,iCell) * basisIntegralsU(iVertexOnCell,jVertexOnCell,iCell)
+
+
+ stressDivergenceUCell = stressDivergenceUCell + &
+ stress1(iVertexOnCell,iCell) * basisIntegralsU(iVertexOnCell,jVertexOnCell,iCell) * 0.5_RKIND + &
+ stress2(iVertexOnCell,iCell) * basisIntegralsU(iVertexOnCell,jVertexOnCell,iCell) * 0.5_RKIND + &
+ stress12(iVertexOnCell,iCell) * basisIntegralsV(iVertexOnCell,jVertexOnCell,iCell)
+
+ stressDivergenceVCell = stressDivergenceVCell + &
+ stress1(iVertexOnCell,iCell) * basisIntegralsV(iVertexOnCell,jVertexOnCell,iCell) * 0.5_RKIND - &
+ stress2(iVertexOnCell,iCell) * basisIntegralsV(iVertexOnCell,jVertexOnCell,iCell) * 0.5_RKIND + &
+ stress12(iVertexOnCell,iCell) * basisIntegralsU(iVertexOnCell,jVertexOnCell,iCell)
+
+ !write(*,*) iVertex, iCellOnVertex, iVertexOnCell, stressDivergenceU(iVertex), stress1(iVertexOnCell,iCell), basisIntegralsU(iVertexOnCell,jVertexOnCell,iCell)
+
+ enddo ! jVertex
+
+
+
+ stressDivergenceU(iVertex) = stressDivergenceU(iVertex) + stressDivergenceUCell / mesh % areaCell % array(iCell)
+ stressDivergenceV(iVertex) = stressDivergenceV(iVertex) + stressDivergenceVCell / mesh % areaCell % array(iCell)
+
+ !write(*,*) iVertex, iVertexOnCell, stressDivergenceUCell / mesh % areaCell % array(iCell), stressDivergenceVCell / mesh % areaCell % array(iCell)
+
+
+ enddo ! iCellOnVertex
+
+ !write(*,*) iVertex, stressDivergenceU(iVertex), stressDivergenceV(iVertex)
+
+
+ stress1contribU = stress1contribU / mesh % areaCell % array(iCell)
+ stress2contribU = stress2contribU / mesh % areaCell % array(iCell)
+ stress12contribU = stress12contribU / mesh % areaCell % array(iCell)
+
+ stress1contribV = stress1contribV / mesh % areaCell % array(iCell)
+ stress2contribV = stress2contribV / mesh % areaCell % array(iCell)
+ stress12contribV = stress12contribV / mesh % areaCell % array(iCell)
+
+ write(*,*) iVertex, iCellOnVertex, stress1contribU, stress12contribU, &
+ stress2contribV, stress12contribV, &
+ stress1contribU+stress2contribU+stress12contribU,stress1contribV+stress2contribV+stress12contribV
+
+
+ !endif ! interiorVertex
+
+ enddo ! iVertex
+
+ write(*,*) "divergence"
+
+ end subroutine divergence_stress_tensor
+
+ !-------------------------------------------------------------
+
+ subroutine test_divergence_stress_tensor(mesh, hexbas, hexdyn, boundary)
+
+ use cice_testing, only: divergence_stress_test_velocity_set, writeout_on_vertex_real, &
+ writeout_on_vertex_on_cell_real, gnuplot_triangle, divergence_stress_test_stress_set_hex, gnuplot_cell
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(hexbas_type), pointer :: hexbas
+ type(hexdyn_type), pointer :: hexdyn
+ type(boundary_type), pointer :: boundary
+
+ !call gnuplot_cell(mesh, mesh % xCell % array, "cells.txt", .false., .false.)
+
+ ! set the velocity field
+ !call divergence_stress_test_velocity_set(hexdyn % uVelocity % array, &
+ ! hexdyn % vVelocity % array, &
+ ! mesh % xVertex % array, &
+ ! mesh % yVertex % array, &
+ ! "div1")
+
+ !call writeout_on_vertex_real(mesh, hexdyn % uVelocity % array, "u.txt")
+ !call writeout_on_vertex_real(mesh, hexdyn % vVelocity % array, "v.txt")
+
+ ! set the stress field from the velocity field with a linear constitutive relation
+ !call stress_tensor(mesh, &
+ ! hexdyn % stress1 % array, hexdyn % stress2 % array, &
+ ! hexdyn % stress12 % array, hexdyn % icePressure % array, &
+ ! hexdyn % uVelocity % array, hexdyn % vVelocity % array, &
+ ! hexbas % basisGradientU % array, hexbas % basisGradientV % array, &
+ ! 0.0_RKIND, 0.0_RKIND)
+
+ call divergence_stress_test_stress_set_hex(mesh, &
+ hexdyn % stress1 % array, &
+ hexdyn % stress2 % array, &
+ hexdyn % stress12 % array)
+
+ !call writeout_on_vertex_on_cell_real(mesh, hexdyn % stress1 % array, "stress1.txt")
+ !call writeout_on_vertex_on_cell_real(mesh, hexdyn % stress2 % array, "stress2.txt")
+ !call writeout_on_vertex_on_cell_real(mesh, hexdyn % stress12 % array, "stress12.txt")
+
+ ! set the divergence of the stress
+ write(*,*) "divergence test"
+ call divergence_stress_tensor(mesh, &
+ hexdyn % stressDivergenceU % array, hexdyn % stressDivergenceV % array, &
+ hexdyn % stress1 % array, hexdyn % stress2 % array, &
+ hexdyn % stress12 % array, &
+ hexbas % basisIntegralsU % array, hexbas % basisIntegralsV % array, &
+ hexbas % cellVerticesAtVertex % array, boundary % interiorVertex % array)
+
+ ! writeout the divergence
+ !call writeout_on_vertex_real(mesh, hexdyn % stressDivergenceU % array, "divu.txt")
+ !call writeout_on_vertex_real(mesh, hexdyn % stressDivergenceV % array, "divv.txt")
+
+ !call gnuplot_triangle(mesh, hexdyn % stressDivergenceU % array, "divugnu.txt")
+
+ stop
+
+ end subroutine test_divergence_stress_tensor
+
+ !-------------------------------------------------------------
+ ! Hex specific plotting stuff
+ !-------------------------------------------------------------
+
+ subroutine plot_wachpress(mesh, iCell, nEdgesOnCell, iVertex, wachspressKappa, wachspressA, wachspressB, xLocal, yLocal)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, intent(in) :: &
+ iCell, &
+ nEdgesOnCell, &
+ iVertex
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ xLocal, &
+ yLocal
+
+ real(kind=RKIND) :: &
+ xmin, xmax, ymin, ymax
+
+ real(kind=RKIND), dimension(2,2) :: &
+ mapping
+
+ integer :: &
+ iSubTriangle, &
+ i1, &
+ i2, &
+ iVertexOnCell, &
+ iVertex2
+
+ xmin = 1e30
+ xmax = -1e30
+ ymin = 1e30
+ ymax = -1e30
+
+ do iVertexOnCell = 1, nEdgesOnCell
+
+ iVertex2 = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ xmin = min(xmin,mesh % xVertex % array(iVertex2)) - 1000.0_RKIND
+ xmax = max(xmax,mesh % xVertex % array(iVertex2)) + 1000.0_RKIND
+ ymin = min(ymin,mesh % yVertex % array(iVertex2)) - 1000.0_RKIND
+ ymax = max(ymax,mesh % yVertex % array(iVertex2)) + 1000.0_RKIND
+
+ open(55,file='test.txt')
+ if (iVertexOnCell == iVertex) then
+ write(55,*) mesh % xVertex % array(iVertex2), mesh % yVertex % array(iVertex2)
+ endif
+ close(55)
+
+ enddo ! iVertex
+
+ open(55,file="wachspress.txt")
+
+ !write(55,fmt='(a,f10.2,a,f10.2,a)') "set xrange [",xmin,":",xmax,"]"
+ !write(55,fmt='(a,f10.2,a,f10.2,a)') "set yrange [",ymin,":",ymax,"]"
+
+ ! loop over subtriangles
+ do iSubTriangle = 1, nEdgesOnCell
+
+ i1 = iSubTriangle
+ i2 = iSubTriangle + 1
+ if (i2 > nEdgesOnCell) i2 = i2 - nEdgesOnCell
+
+ ! get the triangle mapping
+ call get_triangle_mapping(mapping, &
+ 1.0_RKIND,0.0_RKIND,&
+ 0.0_RKIND,1.0_RKIND,&
+ xLocal(i1),yLocal(i1),&
+ xLocal(i2),yLocal(i2))
+
+ ! plot from subtriangle
+ call plot_subtriangle(nEdgesOnCell, iSubTriangle, iVertex, &
+ wachspressKappa, wachspressA, wachspressB, &
+ mapping, mesh % xCell % array(iCell), mesh % yCell % array(iCell))
+
+ enddo ! iSubTriangle
+
+ close(55)
+stop
+ end subroutine plot_wachpress
+
+ !-------------------------------------------------------------
+
+ subroutine plot_subtriangle(nEdgesOnCell, iSubTriangle, iVertex, wachspressKappa, wachspressA, wachspressB, mapping, x0, y0)
+
+ use cice_testing, only: matlab_jet
+
+ integer, intent(in) :: &
+ nEdgesOnCell, &
+ iSubTriangle, &
+ iVertex
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ wachspressKappa
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ wachspressA, &
+ wachspressB
+
+ real(kind=RKIND), dimension(2,2), intent(in) :: &
+ mapping
+
+ real(kind=RKIND), intent(in) :: &
+ x0, y0
+
+
+ real(kind=RKIND) :: &
+ u, v, &
+ x, y, &
+ wachspress, &
+ wachspress1
+
+ integer :: &
+ i, j, iObject
+
+ real(kind=RKIND) :: &
+ x1, x2, x3, x4, x5, &
+ y1, y2, y3, y4, y5, &
+ d, minf, maxf, fValue, &
+ dx, dy
+
+ real(kind=RKIND), parameter :: &
+ fMin = -1.7759450831950581E-004_RKIND, &
+ fMax = 5.5313696223098112E-005_RKIND
+
+ integer, parameter :: n = 50
+
+ character(len=7) :: &
+ color
+
+ logical, parameter :: &
+ lfinitedifference = .true.
+
+ iObject = ((iSubTriangle - 1) * (n+1) * (n+2)) / 2
+
+ d = 1.0_RKIND / real(n, RKIND)
+
+ minf = 1e30
+ maxf = -1e30
+
+ do i = 0, n
+ do j = 0, n
+
+ if (i<=n-j-1) then
+
+ iObject = iObject + 1
+
+ u = real(i,RKIND) / real(n,RKIND)
+ v = real(j,RKIND) / real(n,RKIND)
+
+ call use_triangle_mapping(x, y, u, v, mapping)
+
+ if (lfinitedifference) then
+ wachspress = &
+ wachspress_basis_function(nEdgesOnCell, iVertex, x, y, wachspressKappa, wachspressA, wachspressB)
+
+ dx = 0.0_RKIND
+ dy = 100.0_RKIND
+
+ x1 = x + dx
+ y1 = y + dy
+
+ wachspress1 = &
+ wachspress_basis_function(nEdgesOnCell, iVertex, x1, y1, wachspressKappa, wachspressA, wachspressB)
+
+ wachspress = (wachspress1 - wachspress) / dy
+
+ else
+
+ wachspress = &
+ wachspress_basis_derivative(nEdgesOnCell, iVertex, x, y, wachspressKappa, wachspressA, wachspressB, 2)
+
+ endif
+
+ fValue = (wachspress - fMin) / (fMax - fMin)
+
+ minf = min(minf,wachspress)
+ maxf = max(maxf,wachspress)
+
+ if (i<n-j-1) then
+
+ call use_triangle_mapping(x1, y1, u, v, mapping)
+ call use_triangle_mapping(x2, y2, u, v+d, mapping)
+ call use_triangle_mapping(x3, y3, u+d, v+d, mapping)
+ call use_triangle_mapping(x4, y4, u+d, v, mapping)
+ x5 = x1
+ y5 = y1
+
+ x1 = x1 + x0 ; y1 = y1 + y0
+ x2 = x2 + x0 ; y2 = y2 + y0
+ x3 = x3 + x0 ; y3 = y3 + y0
+ x4 = x4 + x0 ; y4 = y4 + y0
+ x5 = x5 + x0 ; y5 = y5 + y0
+
+ !write(55,fmt='(a,i5,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2)') "set object ",iObject," polygon from ",&
+ ! x1,",",y1," to ",x2,",",y2," to ",x3,",",y3," to ",x4,",",y4," to ",x5,",",y5
+
+ color = matlab_jet(max(min(fValue,1.0_RKIND),0.0_RKIND))
+ !write(55,fmt='(a,i5,a,a,a)') "set object ",iObject,' fc rgb "', color, '" fillstyle solid'
+
+ else
+
+ call use_triangle_mapping(x1, y1, u, v, mapping)
+ call use_triangle_mapping(x2, y2, u, v+d, mapping)
+ call use_triangle_mapping(x4, y4, u+d, v, mapping)
+ x5 = x1
+ y5 = y1
+
+ x1 = x1 + x0 ; y1 = y1 + y0
+ x2 = x2 + x0 ; y2 = y2 + y0
+ x4 = x4 + x0 ; y4 = y4 + y0
+ x5 = x5 + x0 ; y5 = y5 + y0
+
+ !write(55,fmt='(a,i5,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2,a,f10.2)') "set object ",iObject," polygon from ",&
+ ! x1,",",y1," to ",x2,",",y2," to ",x4,",",y4," to ",x5,",",y5
+
+ color = matlab_jet(max(min(fValue,1.0_RKIND),0.0_RKIND))
+ !write(55,fmt='(a,i5,a,a,a)') "set object ",iObject,' fc rgb "', color, '" fillstyle solid'
+
+ endif
+
+ endif
+
+ enddo ! j
+ enddo ! i
+
+ !write(0,*) minf, maxf
+
+ end subroutine plot_subtriangle
+
+ !-------------------------------------------------------------
+
+end module cice_dynamics_evp_hex
Added: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_evp_tri.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_evp_tri.F (rev 0)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_evp_tri.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -0,0 +1,1064 @@
+module cice_dynamics_evp_tri
+
+ use mpas_grid_types
+
+ implicit none
+
+ private
+ public :: init_dynamics_evp, &
+ test_tri_divergence_stress_tensor
+
+contains
+
+
+
+ !-------------------------------------------------------------
+ ! Initialization
+ !-------------------------------------------------------------
+
+ subroutine init_dynamics_evp(mesh, basis, boundary)
+
+ use cice_testing, only: gnuplot_cell, gnuplot_triangle, plot_boundary_triangles, &
+ boundary_locations, labels_verticesOnCell, labels_verticesDegreeOnVertex, &
+ writeout_minmax
+
+ use cice_dynamics_shared, only: boundary_cells, boundary_vertices, boundary_edges
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(basis_type), pointer :: basis
+ type(boundary_type), pointer :: boundary
+
+ integer :: iCell, iVertex
+
+ call boundary_cells(mesh, &
+ boundary % boundaryCell % array, &
+ boundary % boundaryCell2 % array)
+
+ call boundary_vertices(mesh, &
+ boundary % boundaryVertex % array, &
+ boundary % interiorVertex % array, &
+ boundary % interiorVertex2 % array, &
+ boundary % boundaryCell % array)
+
+ call boundary_edges(mesh, &
+ boundary % boundaryEdge % array)
+
+ !call gnuplot_cell(mesh, mesh % xCell % array, "cells.txt", .false., .true.)
+ call gnuplot_triangle(mesh, mesh % xVertex % array+mesh % yVertex % array, "triangles.txt", .false., .true.)
+
+ call writeout_minmax()
+
+ !call labels_verticesOnCell(mesh, mesh % verticesOnCell % array, "cell_vertices.txt")
+ !call labels_verticesDegreeOnVertex(mesh, mesh % cellsOnVertex % array, boundary % interiorVertex % array, "triangle_corners.txt")
+
+ call calc_triangle_corner_at_cell_center(mesh, basis % triangleCornerAtCellCenter % array)
+ call calc_edge_opposite_triangle_corner(mesh, basis % edgeOppositeTriangleCorner % array, boundary % interiorVertex % array)
+
+ !call labels_verticesOnCell(mesh, basis % triangleCornerAtCellCenter % array, "triangleCornerAtCellCenter.txt")
+ !call labels_verticesDegreeOnVertex(mesh, basis % edgeOppositeTriangleCorner % array, boundary % interiorVertex % array, "edgeOppositeTriangleCorner.txt")
+
+
+ !stop
+
+
+ call basis_gradients_planar(mesh, &
+ basis % basisGradientU % array, &
+ basis % basisGradientV % array, &
+ boundary % interiorVertex % array)
+
+ call basis_integrals(mesh, &
+ basis % basisIntegrals % array, &
+ basis % edgeOppositeTriangleCorner % array)
+
+ ! plotting stuff
+ !call gnuplot_cell(mesh, mesh % xCell % array, "xCell.txt")
+ !call gnuplot_triangle(mesh, mesh % xVertex % array, "xVertex.txt")
+
+ !call boundary_locations(mesh, &
+ ! boundary % boundaryCell % array, &
+ ! boundary % boundaryVertex % array, &
+ ! boundary % interiorVertex % array, &
+ ! boundary % boundaryEdge % array)
+
+ !call plot_boundary_triangles(mesh, &
+ ! boundary % boundaryVertex % array, &
+ ! boundary % interiorVertex % array, &
+ ! boundary % boundaryEdge % array)
+
+ end subroutine init_dynamics_evp
+
+ !-------------------------------------------------------------
+
+ subroutine calc_triangle_corner_at_cell_center(mesh, triangleCornerAtCellCenter)
+
+ ! this finds the triangle corner number of a cell centre for an adjacent triangle
+
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:,:), intent(out) :: &
+ triangleCornerAtCellCenter
+
+ integer :: &
+ iCell, & !
+ iVertexOnCell, & !
+ iVertex, & !
+ iVertexDegree, & !
+ iCellVertexDegree !
+
+ ! loop over velocity positions
+ do iCell = 1, mesh % nCells
+
+ ! loop over triangles surrounding velocity position
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ ! get the vertex number associated with this cell vertex
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell, iCell)
+
+ ! loop over triangle corners of this triangle
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ ! get the cell number associated with this triangle corner
+ iCellVertexDegree = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ ! test to see if it is the cell number of interest
+ if (iCellVertexDegree == iCell) then
+
+ triangleCornerAtCellCenter(iVertexOnCell,iCell) = iVertexDegree
+ exit
+
+ endif
+
+ enddo ! iVertexDegree
+
+ enddo ! iVertexOnCell
+
+ enddo ! iCell
+
+ end subroutine calc_triangle_corner_at_cell_center
+
+ !-------------------------------------------------------------
+
+ subroutine calc_edge_opposite_triangle_corner(mesh, edgeOppositeTriangleCorner, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:,:), intent(out) :: &
+ edgeOppositeTriangleCorner
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell, &
+ jVertexDegree, &
+ jEdge, &
+ jEdgeOnCell
+
+ integer, dimension(3) :: &
+ edgeMatched, &
+ jEdgeVertexDegree
+
+ ! loop over all triangles
+ do iVertex = 1, mesh % nVertices
+
+ edgeOppositeTriangleCorner(:,iVertex) = 0
+
+ ! check are at an interior triangle
+ if (interiorVertex(iVertex) == 1) then
+
+ ! find opposite edge for each corner of triangle
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ ! get cell number of this triangle corner
+ iCell = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ ! initialize the matched array to none
+ edgeMatched = 0
+
+ ! check each triangle corner to see if NOT on iCell
+ do jVertexDegree = 1, mesh % vertexDegree
+
+ ! edge of this triangle corner
+ jEdgeVertexDegree(jVertexDegree) = mesh % edgesOnVertex % array(jVertexDegree, iVertex)
+
+ ! loop over edges on this cell see whether match with
+ do jEdge = 1, mesh % nEdgesOnCell % array(iCell)
+
+ ! edge number of this edge on the cell
+ jEdgeOnCell = mesh % edgesOnCell % array(jEdge,iCell)
+
+ ! see if edges are the same
+ if (jEdgeVertexDegree(jVertexDegree) == jEdgeOnCell) then
+
+ ! matched the vertex edge so this ISNT the right edge
+ edgeMatched(jVertexDegree) = 1
+
+ endif
+
+ enddo ! iEdgeOnCell
+
+ enddo ! jVertexDegree
+
+ ! check have matched right number
+ if (sum(edgeMatched) /= 2) then
+ write(*,*) "calc_edge_opposite_triangle_corner: matched failed"
+ stop
+ endif
+
+ ! find the right nonmatched edge
+ do jVertexDegree = 1, mesh % vertexDegree
+
+ if (edgeMatched(jVertexDegree) == 0) then
+ edgeOppositeTriangleCorner(iVertexDegree,iVertex) = jEdgeVertexDegree(jVertexDegree)
+ endif
+
+ enddo ! jVertexDegree
+
+ enddo ! iVertexDegree
+
+ endif ! interiorVertex
+
+ enddo ! iVertex
+
+ end subroutine calc_edge_opposite_triangle_corner
+
+ !-------------------------------------------------------------
+
+ subroutine basis_gradients_planar(mesh, basisGradientU, basisGradientV, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ basisGradientU, &
+ basisGradientV
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ integer :: &
+ iVertex, &
+ iVertexDegree
+
+ integer, dimension(3) :: &
+ iCells
+
+ write(*,*) "basis_gradients_planar"
+
+ ! loop over the triangles
+ do iVertex = 1, mesh % nVertices
+
+ if (interiorVertex(iVertex) == 1) then
+
+ ! loop over corners of triangle
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCells(iVertexDegree) = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ enddo ! iVertexDegree
+
+ call basis_gradients_planar_single(&
+ basisGradientU(:,iVertex), &
+ basisGradientV(:,iVertex), &
+ mesh % xVertex % array(iVertex), mesh % yVertex % array(iVertex), &
+ mesh % xCell % array(iCells(1)), mesh % yCell % array(iCells(1)), &
+ mesh % xCell % array(iCells(2)), mesh % yCell % array(iCells(2)), &
+ mesh % xCell % array(iCells(3)), mesh % yCell % array(iCells(3)))
+
+ endif ! iVertex
+
+ enddo
+
+ end subroutine basis_gradients_planar
+
+ !-------------------------------------------------------------
+
+ subroutine basis_gradients_planar_single(basisGradientU, basisGradientV, &
+ x0, y0, x1, y1, x2, y2, x3, y3)
+
+ use cice_dynamics_shared, only: solve_linear_basis_system
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ basisGradientU, &
+ basisGradientV
+
+ real(kind=RKIND), intent(in) :: &
+ x0, y0, &
+ x1, y1, &
+ x2, y2, &
+ x3, y3
+
+ real(kind=RKIND), dimension(3,3) :: &
+ leftMatrix
+
+ real(kind=RKIND), dimension(3) :: &
+ rightHandSide, &
+ solutionVector
+
+ integer :: &
+ iTriangleCorner
+
+ leftMatrix(1,1) = x1 - x0
+ leftMatrix(1,2) = y1 - y0
+ leftMatrix(1,3) = 1.0_RKIND
+
+ leftMatrix(2,1) = x2 - x0
+ leftMatrix(2,2) = y2 - y0
+ leftMatrix(2,3) = 1.0_RKIND
+
+ leftMatrix(3,1) = x3 - x0
+ leftMatrix(3,2) = y3 - y0
+ leftMatrix(3,3) = 1.0_RKIND
+
+ ! loop over corners of triangle
+ do iTriangleCorner = 1, 3
+
+ rightHandSide(:) = 0.0_RKIND
+ rightHandSide(iTriangleCorner) = 1.0_RKIND
+
+ call solve_linear_basis_system(leftMatrix, rightHandSide, solutionVector)
+
+ basisGradientU(iTriangleCorner) = solutionVector(1)
+ basisGradientV(iTriangleCorner) = solutionVector(2)
+
+ enddo ! iVertexDegree
+
+ end subroutine basis_gradients_planar_single
+
+ !-------------------------------------------------------------
+
+ subroutine basis_integrals(mesh, basisIntegral, edgeOppositeTriangleCorner)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ basisIntegral
+
+ integer, dimension(:,:), intent(in) :: &
+ edgeOppositeTriangleCorner
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iEdge
+
+ real(kind=RKIND) :: &
+ a, b, c
+
+ write(*,*) "basis_integrals"
+
+ ! loop over triangles
+ do iVertex = 1, mesh % nVertices
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iEdge = edgeOppositeTriangleCorner(modulo(iVertexDegree-1,3)+1, iVertex)
+ c = mesh % dcEdge % array(iEdge)
+
+ iEdge = edgeOppositeTriangleCorner(modulo(iVertexDegree,3)+1, iVertex)
+ a = mesh % dcEdge % array(iEdge)
+
+ iEdge = edgeOppositeTriangleCorner(modulo(iVertexDegree-2,3)+1, iVertex)
+ b = mesh % dcEdge % array(iEdge)
+
+ basisIntegral(iVertexDegree,iVertex) = &
+ basis_integral(a, b, c, mesh % areaTriangle % array(iVertex))
+
+ !write(*,*) iVertex, iVertexDegree, basisIntegral(iVertexDegree,iVertex)
+
+ enddo ! iVertexDegree
+
+ enddo ! iVertex
+
+ end subroutine basis_integrals
+
+ !-------------------------------------------------------------
+
+ function basis_integral(a, b, c, areaTriangle) result(integral)
+
+ real(kind=RKIND), intent(in) :: &
+ a, &
+ b, &
+ c, &
+ areaTriangle
+
+ real(kind=RKIND) :: &
+ integral
+
+ real(kind=RKIND) :: &
+ d, &
+ cosGamma
+
+ cosGamma = (a**2 + b**2 - c**2) / (2.0_RKIND * a * b)
+
+ d = ((a*b) / c) * sqrt(1.0_RKIND - cosGamma**2)
+
+ !integral = (c*d) / (6.0_RKIND * areaTriangle)
+ integral = (c*d) / 6.0_RKIND
+
+ !write(*,*) a, b, c, d, areaTriangle, integral
+
+ end function basis_integral
+
+ !-------------------------------------------------------------
+
+ function greatCircleLength(x1,y1,z1,x2,y2,z2) result(length)
+
+ real(kind=RKIND), intent(in) :: &
+ x1, y1, z1, &
+ x2, y2, z2
+
+ real(kind=RKIND) :: length
+
+ real(kind=RKIND), dimension(3) :: &
+ crossProduct
+
+ real(kind=RKIND) :: &
+ dotProduct, &
+ absCrossProduct
+
+ !http://en.wikipedia.org/wiki/Great-circle_distance#Vector_version
+
+ dotProduct = x1 * x2 + y1 * y2 + z1 * z2
+
+ crossProduct(1) = y1 * z2 - z1 * y2
+ crossProduct(2) = z1 * x2 - x1 * z2
+ crossProduct(3) = x1 * y2 - y1 * x2
+
+ absCrossProduct = sqrt(crossProduct(1)**2 + crossProduct(2)**2 + crossProduct(3)**2)
+
+ length = atan(absCrossProduct / dotProduct) ! atan or atan2?
+
+ end function greatCircleLength
+
+ !-------------------------------------------------------------
+
+ function planarLength(x1,y1,z1,x2,y2,z2) result(length)
+
+ real(kind=RKIND), intent(in) :: &
+ x1, y1, z1, &
+ x2, y2, z2
+
+ real(kind=RKIND) :: length
+
+ length = sqrt((x2-x1)**2 + (y2-y1)**2 + (z2-z1)**2)
+
+ end function planarLength
+
+ !-------------------------------------------------------------
+
+ function sphericalTriangleArea(radius,a,b,c) result(area)
+
+ use mpas_constants, only: pii
+
+ real(kind=RKIND), intent(in) :: &
+ radius, a, b, c
+
+ real(kind=RKIND) :: area
+
+ real(kind=RKIND) :: &
+ sphericalExcess, &
+ s
+
+ s = 0.5_RKIND * (a + b + c)
+
+ sphericalExcess = 4.0_RKIND * atan(sqrt(tan(s/2.0_RKIND)*tan((s-a)/2.0_RKIND)*tan((s-b)/2.0_RKIND)*tan((s-c)/2.0_RKIND)))
+
+ area = (pii * radius**2 * sphericalExcess) / 180.0_RKIND
+
+ end function sphericalTriangleArea
+
+ !-------------------------------------------------------------
+
+ function planarTriangleArea(radius,a,b,c) result(area)
+
+ real(kind=RKIND), intent(in) :: &
+ radius, a, b, c
+
+ real(kind=RKIND) :: area
+
+ !http://en.wikipedia.org/wiki/Triangle#Using_Heron.27s_formula
+
+ real(kind=RKIND) :: &
+ s
+
+ s = 0.5_RKIND * (a + b + c)
+
+ area = sqrt(s * (s-a) * (s-b) * (s-c))
+
+ end function planarTriangleArea
+
+ !-------------------------------------------------------------
+ ! time step
+ !-------------------------------------------------------------
+
+ subroutine run_dynamics_evp(mesh, icestate, dynamics, basis, forcing, boundary)
+
+ use cice_dynamics_shared, only: ice_strength, air_stress, ocean_stress, surface_tilt, solve_velocity
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(dynamics_type), pointer :: dynamics
+ type(basis_type), pointer :: basis
+ type(forcing_type), pointer :: forcing
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND) :: timeStep, elasticTimeStep
+
+ call ice_strength(mesh, &
+ dynamics % icePressure % array, &
+ icestate % iceAreaCell % array, &
+ icestate % iceVolumeCell % array)
+
+ call stress_tensor(mesh, &
+ dynamics % stress1 % array, &
+ dynamics % stress2 % array, &
+ dynamics % stress12 % array, &
+ dynamics % icePressure % array, &
+ dynamics % uVelocity % array, &
+ dynamics % vVelocity % array, &
+ basis % basisGradientU % array, &
+ basis % basisGradientV % array, &
+ boundary % interiorVertex2 % array, &
+ timeStep, elasticTimeStep)
+
+ call divergence_stress_tensor(mesh, &
+ dynamics % stressDivergenceU % array, &
+ dynamics % stressDivergenceV % array, &
+ dynamics % stress1 % array, &
+ dynamics % stress2 % array, &
+ dynamics % stress12 % array, &
+ basis % basisGradientU % array, &
+ basis % basisGradientV % array, &
+ basis % basisIntegrals % array, &
+ basis % triangleCornerAtCellCenter % array, &
+ boundary % boundaryCell % array, &
+ boundary % interiorVertex2 % array)
+
+ call air_stress(&
+ dynamics % airStressU % array, &
+ dynamics % airStressV % array, &
+ forcing % uAirVelocity % array, &
+ forcing % vAirVelocity % array, &
+ icestate % iceAreaVertex % array)
+
+ !call ocean_stress(&
+ ! dynamics % oceanStressU % array, &
+ ! dynamics % oceanStressV % array, &
+ ! dynamics % oceanStressCoeff % array, &
+ ! forcing % uOceanVelocity % array, &
+ ! forcing % vOceanVelocity % array, &
+ ! dynamics % uVelocity % array, &
+ ! dynamics % vVelocity % array, &
+ ! icestate % iceAreaCell % array)
+
+ call surface_tilt(&
+ dynamics % surfaceTiltForceU % array, &
+ dynamics % surfaceTiltForceV % array)
+
+ call solve_velocity(mesh % nCells, &
+ boundary % interiorVertex % array, &
+ dynamics % uVelocity % array, &
+ dynamics % vVelocity % array, &
+ iceState % totalMassCell % array, &
+ mesh % fCell % array, &
+ dynamics % stressDivergenceU % array, &
+ dynamics % stressDivergenceV % array, &
+ dynamics % airStressU % array, &
+ dynamics % airStressV % array, &
+ dynamics % surfaceTiltForceU % array, &
+ dynamics % surfaceTiltForceV % array, &
+ dynamics % oceanStressU % array, &
+ dynamics % oceanStressV % array, &
+ dynamics % oceanStressCoeff % array, &
+ timeStep)
+
+ end subroutine run_dynamics_evp
+
+ !-------------------------------------------------------------
+
+ subroutine stress_tensor(mesh, &
+ stress1, stress2, &
+ stress12, &
+ icePressure, &
+ uVelocity, vVelocity, &
+ basisGradientU, basisGradientV, &
+ interiorVertex, &
+ timeStep, elasticTimeStep)
+
+ use cice_dynamics_shared, only: evp_constitutive_relation
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uVelocity, &
+ vVelocity
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ basisGradientU, &
+ basisGradientV
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ icePressure
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ real(kind=RKIND), intent(in) :: &
+ timeStep, &
+ elasticTimeStep
+
+ real(kind=RKIND) :: &
+ strainDivergence, & ! divergence of the strain rate tensor
+ strainTension, & ! horizontal tension of the strain rate tensor
+ strainShearing ! horizontal shearing of the strain rate tensor
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell
+
+ ! loop over the triangles
+ do iVertex = 1, mesh % nVertices
+
+ if (interiorVertex(iVertex) == 1) then
+
+ strainDivergence = 0.0_RKIND
+ strainTension = 0.0_RKIND
+ strainShearing = 0.0_RKIND
+
+ ! loop over corners of triangle
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ ! divergence of the strain rate tensor
+ strainDivergence = strainDivergence + &
+ uVelocity(iCell) * basisGradientU(iVertexDegree,iVertex) + &
+ vVelocity(iCell) * basisGradientV(iVertexDegree,iVertex)
+
+ ! horizontal tension of the strain rate tensor
+ strainTension = strainTension + &
+ uVelocity(iCell) * basisGradientU(iVertexDegree,iVertex) - &
+ vVelocity(iCell) * basisGradientV(iVertexDegree,iVertex)
+
+ ! horizontal shearing of the strain rate tensor
+ strainShearing = strainShearing + &
+ uVelocity(iCell) * basisGradientV(iVertexDegree,iVertex) + &
+ vVelocity(iCell) * basisGradientU(iVertexDegree,iVertex)
+
+ enddo ! iVertexDegree
+
+ ! loop over corners of triangle
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ !call evp_constitutive_relation(stress1(iVertexDegree,iVertex), &
+ ! stress2(iVertexDegree,iVertex), &
+ ! stress12(iVertexDegree,iVertex), &
+ ! strainDivergence, strainTension, strainShearing, &
+ ! icePressure(iCell), mesh % areaCell % array(iCell), &
+ ! timeStep, elasticTimeStep)
+
+ enddo ! iVertexDegree
+
+ endif ! interiorVertex
+
+ enddo ! iVertex
+
+ end subroutine stress_tensor
+
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_tensor(mesh, &
+ divergenceStressU, divergenceStressV, &
+ stress1, stress2, &
+ stress12, &
+ basisGradientU, basisGradientV, &
+ basisIntegrals, &
+ triangleCornerAtCellCenter, &
+ boundaryCell, interiorVertex2)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ divergenceStressU, &
+ divergenceStressV
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ stress1, &
+ stress2, &
+ stress12, &
+ basisGradientU, &
+ basisGradientV, &
+ basisIntegrals
+
+ integer, dimension(:,:), intent(in) :: &
+ triangleCornerAtCellCenter
+
+ integer, dimension(:), intent(in) :: &
+ boundaryCell, &
+ interiorVertex2
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ iVertex, &
+ iVertexDegree, &
+ jVertexDegree
+
+
+ real(kind=RKIND) :: &
+ stress1contribU, &
+ stress2contribU, &
+ stress12contribU, &
+ stress1contribV, &
+ stress2contribV, &
+ stress12contribV
+
+ write(*,*) "divergence start", mesh % nCells
+
+ ! loop over velocity positions
+ do iCell = 1, mesh % nCells
+
+ ! only cells interior to our smaller domain
+ if (boundaryCell(iCell) == 0) then
+
+ divergenceStressU(iCell) = 0.0_RKIND
+ divergenceStressV(iCell) = 0.0_RKIND
+
+ ! loop over triangles surrounding velocity position
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ ! vertex index of the current triangle
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell, iCell)
+
+ stress1contribU = 0.0_RKIND
+ stress2contribU = 0.0_RKIND
+ stress12contribU = 0.0_RKIND
+
+ stress1contribV = 0.0_RKIND
+ stress2contribV = 0.0_RKIND
+ stress12contribV = 0.0_RKIND
+
+ ! only internal triangles at the mo
+ !if (interiorVertex2(iVertex) == 1) then
+
+ ! vertex number of vertex triangle from this cell
+ jVertexDegree = triangleCornerAtCellCenter(iVertexOnCell, iCell)
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ stress1contribU = stress1contribU + stress1(iVertexDegree,iVertex) * basisGradientU(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+ stress2contribU = 0.0_RKIND
+ stress12contribU = stress12contribU + stress12(iVertexDegree,iVertex) * basisGradientV(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+
+ stress1contribV = 0.0_RKIND
+ stress2contribV = stress2contribV + stress2(iVertexDegree,iVertex) * basisGradientV(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+ stress12contribV = stress12contribV + stress12(iVertexDegree,iVertex) * basisGradientU(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+
+ ! divergence of the stress tensor in u direction
+ !divergenceStressU(iCell) = divergenceStressU(iCell) + &
+ ! stress1(iVertexDegree,iVertex) * basisGradientU(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ ! stress2(iVertexDegree,iVertex) * basisGradientU(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ ! stress12(iVertexDegree,iVertex) * basisGradientV(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+
+ ! divergence of the stress tensor in v direction
+ !divergenceStressV(iCell) = divergenceStressV(iCell) + &
+ ! stress1(iVertexDegree,iVertex) * basisGradientV(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ ! stress2(iVertexDegree,iVertex) * basisGradientV(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ ! stress12(iVertexDegree,iVertex) * basisGradientU(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+
+ write(*,*) iCell, iVertexOnCell, iVertexDegree, stress1(iVertexDegree,iVertex), basisGradientU(iVertexDegree,iVertex), basisIntegrals(iVertexDegree,iVertex)
+
+ enddo ! iVertexDegree
+
+ !endif ! interiorVertex2
+
+ divergenceStressU(iCell) = divergenceStressU(iCell) + &
+ (stress1contribU + stress12contribU) / mesh % areaTriangle % array(iVertex)
+
+ divergenceStressV(iCell) = divergenceStressV(iCell) + &
+ (stress2contribU + stress12contribU) / mesh % areaTriangle % array(iVertex)
+
+ enddo ! iVertexOnCell
+
+ write(*,*) iCell, divergenceStressU(iCell), divergenceStressV(iCell)
+
+ endif ! boundaryCell
+
+ enddo ! iCell
+
+ write(*,*) "divergence"
+
+ end subroutine divergence_stress_tensor
+
+ !-------------------------------------------------------------
+
+ subroutine test_tri_divergence_stress_tensor(mesh, dynamics, basis, boundary)
+
+ use cice_testing, only: divergence_stress_test_stress_set_tri
+
+ type(mesh_type), intent(in) :: mesh
+ type(dynamics_type), pointer :: dynamics
+ type(basis_type), pointer :: basis
+ type(boundary_type), pointer :: boundary
+
+ write(*,*) "divergence test"
+
+ call divergence_stress_test_stress_set_tri(mesh, &
+ dynamics % stress1 % array, &
+ dynamics % stress2 % array, &
+ dynamics % stress12 % array)
+
+ call divergence_stress_tensor(mesh, &
+ dynamics % stressDivergenceU % array, &
+ dynamics % stressDivergenceV % array, &
+ dynamics % stress1 % array, &
+ dynamics % stress2 % array, &
+ dynamics % stress12 % array, &
+ basis % basisGradientU % array, &
+ basis % basisGradientV % array, &
+ basis % basisIntegrals % array, &
+ basis % triangleCornerAtCellCenter % array, &
+ boundary % boundaryCell % array, &
+ boundary % interiorVertex % array)
+
+ stop
+
+ end subroutine test_tri_divergence_stress_tensor
+
+ !-------------------------------------------------------------
+ ! boundary
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_tensor_boundary(mesh, &
+ divergenceStressU, divergenceStressV, &
+ uVelocity, vVelocity, &
+ icePressure, basisIntegrals, &
+ basisGradientU, basisGradientV, &
+ boundaryCell, boundaryCell2, &
+ boundaryVertex2, &
+ timeStep, elasticTimeStep)
+
+ use cice_dynamics_shared, only: evp_constitutive_relation
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(inout) :: &
+ divergenceStressU, &
+ divergenceStressV
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uVelocity, &
+ vVelocity, &
+ icePressure
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ basisIntegrals, &
+ basisGradientU, &
+ basisGradientV
+
+ integer, dimension(:), intent(in) :: &
+ boundaryCell, &
+ boundaryCell2, &
+ boundaryVertex2
+
+ real(kind=RKIND), intent(in) :: &
+ timeStep, &
+ elasticTimeStep
+
+ real(kind=RKIND) :: &
+ strainDivergence, &
+ strainTension, &
+ strainShearing, &
+ stress1, &
+ stress2, &
+ stress12, &
+ divergenceStressUBoundary, &
+ divergenceStressVBoundary, &
+ areaCorrection
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ iVertex, &
+ iVertexDegree, &
+ iCellTriangleCorner, &
+ jVertexDegree
+
+ ! loop over all velocity points
+ do iCell = 1, mesh % nCells
+
+ ! interior boundary velocity points
+ if (boundaryCell2(iCell) == 1) then
+
+ ! find the triangles we still need to sum
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ if (boundaryVertex2(iVertex) == 1) then
+
+ strainDivergence = 0.0_RKIND
+ strainTension = 0.0_RKIND
+ strainShearing = 0.0_RKIND
+
+ areaCorrection = 0.0_RKIND
+
+ ! loop over triangle vertices
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCellTriangleCorner = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ if (boundaryCell(iCellTriangleCorner) == 0) then
+ ! interior triangle corner
+ areaCorrection = areaCorrection + mesh % kiteAreasOnVertex % array(iVertexDegree,iVertex)
+
+ ! divergence of the strain rate tensor
+ strainDivergence = strainDivergence + &
+ uVelocity(iCellTriangleCorner) * basisGradientU(iVertexDegree,iVertex) + &
+ vVelocity(iCellTriangleCorner) * basisGradientV(iVertexDegree,iVertex)
+
+ ! horizontal tension of the strain rate tensor
+ strainTension = strainTension + &
+ uVelocity(iCellTriangleCorner) * basisGradientU(iVertexDegree,iVertex) - &
+ vVelocity(iCellTriangleCorner) * basisGradientV(iVertexDegree,iVertex)
+
+ ! horizontal shearing of the strain rate tensor
+ strainShearing = strainShearing + &
+ uVelocity(iCellTriangleCorner) * basisGradientV(iVertexDegree,iVertex) + &
+ vVelocity(iCellTriangleCorner) * basisGradientU(iVertexDegree,iVertex)
+
+ else
+ ! exterior triangle corner
+
+ ! divergence of the strain rate tensor
+ strainDivergence = strainDivergence - &
+ uVelocity(iCell) * basisGradientU(iVertexDegree,iVertex) - &
+ vVelocity(iCell) * basisGradientV(iVertexDegree,iVertex)
+
+ ! horizontal tension of the strain rate tensor
+ strainTension = strainTension - &
+ uVelocity(iCell) * basisGradientU(iVertexDegree,iVertex) + &
+ vVelocity(iCell) * basisGradientV(iVertexDegree,iVertex)
+
+ ! horizontal shearing of the strain rate tensor
+ strainShearing = strainShearing - &
+ uVelocity(iCell) * basisGradientV(iVertexDegree,iVertex) - &
+ vVelocity(iCell) * basisGradientU(iVertexDegree,iVertex)
+
+ endif ! boundaryCell
+
+ enddo ! iVertexDegree
+
+ areaCorrection = areaCorrection / mesh % areaCell % array(iCell)
+
+ ! calculate stress assuming all ice strength equal that at iCell
+ !call evp_constitutive_relation(stress1, stress2, stress12, &
+ ! strainDivergence, strainTension, strainShearing, &
+ ! icePressure(iCell), mesh % areaCell % array(iCell), &
+ ! timeStep, elasticTimeStep)
+
+ ! vertex number of vertex triangle from this cell
+ !jVertexDegree = triangleCornerAtCellCenter(iVertexOnCell, iCell) ????
+
+ divergenceStressUBoundary = 0.0_RKIND
+ divergenceStressVBoundary = 0.0_RKIND
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ ! divergence of the stress tensor in u direction
+ divergenceStressUBoundary = divergenceStressUBoundary + &
+ stress1 * basisGradientU(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ stress2 * basisGradientU(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ stress12 * basisGradientV(iVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+
+ ! divergence of the stress tensor in v direction
+ divergenceStressVBoundary = divergenceStressVBoundary + &
+ stress1 * basisGradientV(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ stress2 * basisGradientV(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex) * 0.5_RKIND + &
+ stress12 * basisGradientU(jVertexDegree,iVertex) * basisIntegrals(iVertexDegree,iVertex)
+
+ enddo ! iVertexDegree
+
+ divergenceStressU(iCell) = divergenceStressU(iCell) + divergenceStressUBoundary * areaCorrection
+ divergenceStressV(iCell) = divergenceStressV(iCell) + divergenceStressVBoundary * areaCorrection
+
+ endif ! boundaryVertex2
+
+ enddo ! iVertexOnCell
+
+ endif ! boundaryCell2
+
+ enddo ! iCell
+
+ end subroutine divergence_stress_tensor_boundary
+
+ !-------------------------------------------------------------
+
+ subroutine set_boundary_velocity(mesh, &
+ uVelocityBoundary, vVelocityBoundary, &
+ uVelocity, vVelocity, &
+ boundaryCell)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(out) :: &
+ uVelocityBoundary, &
+ vVelocityBoundary
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uVelocity, &
+ vVelocity
+
+ integer, dimension(:), intent(in) :: &
+ boundaryCell
+
+ integer :: &
+ iCell, &
+ iEdgeOnCell, &
+ iCellNeighbour
+
+ ! loop over all possible velocity points
+ do iCell = 1, mesh % nCells
+
+ ! only boundary velocity points
+ if (boundaryCell(iCell) == 1) then
+
+ ! loop over surrounding cells
+ do iEdgeOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ ! cell number of neighbouring cell
+ iCellNeighbour = mesh % cellsOnCell % array(iEdgeOnCell,iCell)
+
+ ! check if this is interior cell
+ if (boundaryCell(iCellNeighbour) == 0) then
+
+ uVelocityBoundary(iEdgeOnCell,iCell) = -uVelocity(iCellNeighbour)
+ vVelocityBoundary(iEdgeOnCell,iCell) = -vVelocity(iCellNeighbour)
+
+ endif ! boundaryCell
+
+ enddo ! iEdgeOnCell
+
+ endif ! boundaryCell
+
+ enddo ! iCell
+
+ end subroutine set_boundary_velocity
+
+ !-------------------------------------------------------------
+
+end module cice_dynamics_evp_tri
+
Added: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_shared.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_shared.F (rev 0)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_dynamics_shared.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -0,0 +1,1448 @@
+module cice_dynamics_shared
+
+ use mpas_grid_types
+
+ implicit none
+
+
+ private
+
+ ! Weak formulation
+ public :: init_weak_dynamics, &
+ run_dynamics, &
+ test_stress_divergence
+
+ ! Legacy
+ public :: boundary_cells, &
+ boundary_vertices, &
+ boundary_edges, &
+ ice_strength, &
+ evp_constitutive_relation, &
+ linear_constitutive_relation, &
+ air_stress, &
+ ocean_stress, &
+ solve_velocity, &
+ surface_tilt, &
+ solve_linear_basis_system
+
+ real(kind=RKIND), parameter, public :: &
+ sinOceanTurningAngle = 0.0_RKIND, &
+ cosOceanTurningAngle = 1.0_RKIND
+
+
+ real(kind=RKIND), parameter :: &
+ Pstar = 2.75e4_RKIND, & ! constant in Hibler strength formula
+ Cstar = 20.0_RKIND , & ! constant in Hibler strength formula
+ eccentricity = 2.0_RKIND, &
+ dampingTimescaleParameter = 0.36_RKIND, &
+ puny = 1.0e-11_RKIND
+
+ real(kind=RKIND) :: &
+ dtDynamics, &
+ dtElastic, &
+ dampingTimescale, &
+ elapsedTime
+
+ integer, parameter :: &
+ nElasticSubcycle = 120 ! will be from namelist
+
+contains
+
+ !-------------------------------------------------------------
+ ! Initialization
+ !-------------------------------------------------------------
+
+ subroutine init_weak_dynamics(mesh, weak, boundary, dt)
+
+ type(mesh_type), intent(inout) :: mesh
+
+ type(weak_type), pointer :: weak
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND), intent(in) :: dt
+
+ call boundary_vertices_weak(mesh, &
+ boundary % interiorVertex % array)
+
+ call normal_vectors_planar(mesh, &
+ weak % normalVectorPolygon % array, &
+ weak % normalVectorTriangle % array, &
+ boundary % interiorVertex % array)
+
+ call init_evp(dt)
+
+ mesh % fVertex % array = 0.0_RKIND
+
+ end subroutine init_weak_dynamics
+
+ !-------------------------------------------------------------
+
+ subroutine boundary_vertices_weak(mesh, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:), intent(out) :: &
+ interiorVertex
+
+ integer :: &
+ nInteriorAdjacentCells
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell
+
+ ! boundary vertices
+ do iVertex = 1, mesh % nVertices
+
+ interiorVertex(iVertex) = 0
+
+ nInteriorAdjacentCells = 0
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ if (iCell >= 1 .and. iCell <= mesh % nCells) then
+
+ nInteriorAdjacentCells = nInteriorAdjacentCells + 1
+
+ endif
+
+ enddo ! iVertexDegree
+
+ ! vertex points we directly calculate on
+ if (nInteriorAdjacentCells == 3) then
+
+ interiorVertex(iVertex) = 1
+
+ endif
+
+ enddo ! iVertex
+
+ end subroutine boundary_vertices_weak
+
+ !-------------------------------------------------------------
+
+ subroutine normal_vectors_planar(mesh, normalVectorPolygon, normalVectorTriangle, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:,:), intent(out) :: &
+ normalVectorPolygon, &
+ normalVectorTriangle
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ call normal_vectors_planar_polygon(mesh, normalVectorPolygon)
+
+ call normal_vectors_planar_triangle(mesh, normalVectorTriangle, interiorVertex)
+
+ end subroutine normal_vectors_planar
+
+ !-------------------------------------------------------------
+
+ subroutine normal_vectors_planar_polygon(mesh, normalVectorPolygon)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:,:), intent(out) :: &
+ normalVectorPolygon
+
+ integer :: &
+ iCell, &
+ iEdgeOnCell, &
+ iEdge
+
+ real(kind=RKIND) :: &
+ dx, dy
+
+ do iCell = 1, mesh % nCells
+
+ do iEdgeOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iEdge = mesh % edgesOnCell % array(iEdgeOnCell,iCell)
+
+ dx = mesh % xEdge % array(iEdge) - mesh % xCell % array(iCell)
+ dy = mesh % yEdge % array(iEdge) - mesh % yCell % array(iCell)
+
+ normalVectorPolygon(1,iEdgeOnCell,iCell) = dx / sqrt(dx**2 + dy**2)
+ normalVectorPolygon(2,iEdgeOnCell,iCell) = dy / sqrt(dx**2 + dy**2)
+
+ enddo ! iEdgeOnCell
+
+ enddo ! iCell
+
+ end subroutine normal_vectors_planar_polygon
+
+ !-------------------------------------------------------------
+
+ subroutine normal_vectors_planar_triangle(mesh, normalVectorTriangle, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:,:), intent(out) :: &
+ normalVectorTriangle
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iEdge
+
+ real(kind=RKIND) :: &
+ dx, dy
+
+ do iVertex = 1, mesh % nVertices
+
+ normalVectorTriangle(:,:,iVertex) = 0.0_RKIND
+
+ if (interiorVertex(iVertex) == 1) then
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iEdge = mesh % edgesOnVertex % array(iVertexDegree,iVertex)
+
+ dx = mesh % xEdge % array(iEdge) - mesh % xVertex % array(iVertex)
+ dy = mesh % yEdge % array(iEdge) - mesh % yVertex % array(iVertex)
+
+ normalVectorTriangle(1,iVertexDegree,iVertex) = dx / sqrt(dx**2 + dy**2)
+ normalVectorTriangle(2,iVertexDegree,iVertex) = dy / sqrt(dx**2 + dy**2)
+
+ enddo ! iVertexDegree
+
+ endif ! interiorVertex
+
+ enddo ! iVertex
+
+ end subroutine normal_vectors_planar_triangle
+
+ !-------------------------------------------------------------
+
+ subroutine init_evp(dt)
+
+ real(kind=RKIND), intent(in) :: dt
+
+ dtDynamics = dt
+
+ dtElastic = dtDynamics / real(nElasticSubcycle,RKIND)
+
+ dampingTimescale = dampingTimescaleParameter * dtDynamics
+
+ elapsedTime = 0.0_RKIND
+
+ end subroutine init_evp
+
+ !-------------------------------------------------------------
+ ! Time stepping
+ !-------------------------------------------------------------
+
+ subroutine run_dynamics(mesh, icestate, weak, boundary, dt)
+
+ use cice_testing, only: gnuplot_triangle, gnuplot_cell, writeout_minmax
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(weak_type), pointer :: weak
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND), intent(in) :: &
+ dt
+
+ call ice_strength(mesh, &
+ weak % icePressure % array, &
+ icestate % iceAreaCell % array, &
+ icestate % iceVolumeCell % array)
+
+ call interpolate_cell_to_vertex(mesh, &
+ icestate % iceAreaVertex % array, &
+ icestate % iceAreaCell % array)
+
+ call air_stress(&
+ weak % airStressU % array, &
+ weak % airStressV % array, &
+ weak % uAirVelocity % array, &
+ weak % vAirVelocity % array, &
+ icestate % iceAreaVertex % array)
+
+ call ocean_stress(&
+ weak % oceanStressU % array, &
+ weak % oceanStressV % array, &
+ weak % uOceanVelocity % array, &
+ weak % vOceanVelocity % array)
+
+ !call surface_tilt(&
+ ! weak % surfaceTiltForceU % array, &
+ ! weak % surfaceTiltForceV % array)
+
+ call subcycle_dynamics(mesh, icestate, weak, boundary, &
+ dtDynamics, dtElastic)
+
+ call principal_stresses(mesh, &
+ weak % principalStress1 % array, &
+ weak % principalStress2 % array, &
+ weak % stress11 % array, &
+ weak % stress22 % array, &
+ weak % stress12 % array, &
+ weak % icePressure % array)
+
+
+ call gnuplot_triangle(mesh, weak % stressDivergenceU % array, "test.txt")
+ !call gnuplot_cell(mesh, weak % uVelocity % array, "test.txt")
+ write(*,*) minval(weak % uVelocity % array), maxval(weak % uVelocity % array)
+ call writeout_minmax()
+ stop
+
+ elapsedTime = elapsedTime + dt
+
+ end subroutine run_dynamics
+
+ !-------------------------------------------------------------
+
+ subroutine subcycle_dynamics(mesh, icestate, weak, boundary, &
+ dtDynamics, dtElastic)
+
+ use cice_testing, only: gnuplot_triangle, gnuplot_cell, writeout_minmax, gnuplot_vertexvector
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(weak_type), pointer :: weak
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND), intent(in) :: &
+ dtDynamics, &
+ dtElastic
+
+ integer :: &
+ iElasticSubcycle
+
+ do iElasticSubcycle = 1, nElasticSubcycle
+
+ write(*,*) "subcycle: ", iElasticSubcycle
+ call single_subcycle_dynamics(mesh, icestate, weak, boundary, &
+ dtDynamics, dtElastic)
+
+ enddo
+
+ call gnuplot_triangle(mesh, weak % stressDivergenceU % array, "test.txt")
+ !call gnuplot_triangle(mesh, weak % uVelocity % array, "test.txt")
+ !call gnuplot_cell(mesh, weak % stress11 % array, "test.txt")
+ stop
+
+ end subroutine subcycle_dynamics
+
+ !-------------------------------------------------------------
+
+ subroutine single_subcycle_dynamics(mesh, icestate, weak, boundary, &
+ dtDynamics, dtElastic)
+
+ use cice_testing, only: gnuplot_triangle, gnuplot_cell, writeout_minmax, gnuplot_vertexvector
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(weak_type), pointer :: weak
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND), intent(in) :: &
+ dtDynamics, &
+ dtElastic
+
+ !weak % uVelocity % array = mesh % xVertex % array
+ !weak % vVelocity % array = 0.0_RKIND
+
+ call strain_tensor(mesh, &
+ weak % strain11 % array, &
+ weak % strain22 % array, &
+ weak % strain12 % array, &
+ weak % uVelocity % array, &
+ weak % vVelocity % array, &
+ weak % normalVectorPolygon % array)
+
+ call stress_tensor(mesh, &
+ weak % stress11 % array, &
+ weak % stress22 % array, &
+ weak % stress12 % array, &
+ weak % strain11 % array, &
+ weak % strain22 % array, &
+ weak % strain12 % array, &
+ weak % icePressure % array, &
+ dtDynamics, &
+ dtElastic)
+
+ call stress_divergence(mesh, &
+ weak % stressDivergenceU % array, &
+ weak % stressDivergenceV % array, &
+ weak % stress11 % array, &
+ weak % stress22 % array, &
+ weak % stress12 % array, &
+ weak % normalVectorTriangle % array, &
+ boundary % interiorVertex % array)
+
+ call interpolate_cell_to_vertex(mesh, &
+ icestate % totalMassVertex % array, &
+ icestate % totalMassCell % array)
+
+ call ocean_stress_coefficient(&
+ weak % oceanStressCoeff % array, &
+ weak % uOceanVelocity % array, &
+ weak % vOceanVelocity % array, &
+ weak % uVelocity % array, &
+ weak % vVelocity % array, &
+ icestate % iceAreaVertex % array)
+
+ call solve_velocity(mesh % nVertices, &
+ boundary % interiorVertex % array, &
+ weak % uVelocity % array, &
+ weak % vVelocity % array, &
+ icestate % totalMassVertex % array, &
+ mesh % fVertex % array, &
+ weak % stressDivergenceU % array, &
+ weak % stressDivergenceV % array, &
+ weak % airStressU % array, &
+ weak % airStressV % array, &
+ weak % surfaceTiltForceU % array, &
+ weak % surfaceTiltForceV % array, &
+ weak % oceanStressU % array, &
+ weak % oceanStressV % array, &
+ weak % oceanStressCoeff % array, &
+ dtElastic)
+
+ end subroutine single_subcycle_dynamics
+
+ !-------------------------------------------------------------
+
+ subroutine ice_strength(mesh, icePressure, &
+ iceAreaCell, iceVolumeCell)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ icePressure
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ iceAreaCell, &
+ iceVolumeCell
+
+ integer :: &
+ iCell
+
+ do iCell = 1, mesh % nCells
+
+ icePressure(iCell) = Pstar * iceVolumeCell(iCell) * exp(-Cstar*(1.0_RKIND-iceAreaCell(iCell)))
+
+ enddo ! iCell
+
+ end subroutine ice_strength
+
+ !-------------------------------------------------------------
+
+ subroutine interpolate_cell_to_vertex(mesh, variableVertex, variableCell)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ variableVertex
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ variableCell
+
+ real(kind=RKIND) :: &
+ totalArea
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell
+
+ do iVertex = 1, mesh % nVertices
+
+ variableVertex(iVertex) = 0.0_RKIND
+ totalArea = 0.0_RKIND
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ variableVertex(iVertex) = variableVertex(iVertex) + mesh % areaCell % array(iCell) * variableCell(iCell)
+
+ totalArea = totalArea + mesh % areaCell % array(iCell)
+
+ enddo ! iVertexDegree
+
+ variableVertex(iVertex) = variableVertex(iVertex) / totalArea
+
+ enddo ! iVertex
+
+ end subroutine interpolate_cell_to_vertex
+
+ !-------------------------------------------------------------
+
+ subroutine evp_constitutive_relation(stress1, stress2, stress12, &
+ strainDivergence, strainTension, strainShearing, &
+ icePressure, areaCell, &
+ dtDynamics, dtElastic,iCell)
+
+ real(kind=RKIND), intent(inout) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND), intent(inout) :: &
+ strainDivergence, &
+ strainTension, &
+ strainShearing
+
+ real(kind=RKIND), intent(in) :: &
+ icePressure, &
+ areaCell, &
+ dtDynamics, &
+ dtElastic
+
+ real(kind=RKIND) :: &
+ Delta, &
+ pressureCoefficient, stress10
+
+ integer :: iCell
+
+ !write(*,*) stress1
+
+ !strainDivergence = 1.0_RKIND
+ !strainTension = 1.0_RKIND
+ !strainShearing = 0.0_RKIND
+
+ !stress1 = 0.0_RKIND
+ !stress2 = 0.0_RKIND
+ !stress12 = 0.0_RKIND
+
+ Delta = sqrt(strainDivergence**2 + &
+ (strainTension**2 + strainShearing**2) / eccentricity**2)
+
+ pressureCoefficient = (icePressure * dtElastic) / (2.0_RKIND * dampingTimescale * max(Delta,puny))
+!stress10 = stress1
+ stress1 = (stress1 + pressureCoefficient * (strainDivergence - Delta)) / &
+ (1.0_RKIND + (0.5_RKIND * dtElastic) / dampingTimescale)
+
+ stress2 = (stress2 + pressureCoefficient * strainTension) / &
+ (1.0_RKIND + (0.5_RKIND * dtElastic * eccentricity**2) / dampingTimescale)
+
+ stress12 = (stress12 + pressureCoefficient * strainShearing * 0.5_RKIND) / &
+ (1.0_RKIND + (0.5_RKIND * dtElastic * eccentricity**2) / dampingTimescale)
+
+ !write(*,*) stress1, icePressure / max(Delta,puny), icePressure
+ !write(*,*) (1.0_RKIND + (0.5_RKIND * dtElastic) / dampingTimescale), dampingTimescale, (0.5_RKIND * dtElastic) / dampingTimescale
+ !stop
+
+ !write(*,*) stress10, stress1, strainDivergence, Delta, icePressure / max(Delta,areaCell * puny), (1.0_RKIND + (0.5_RKIND * dtElastic) / dampingTimescale), (0.5_RKIND * dtElastic) / dampingTimescale
+
+ end subroutine evp_constitutive_relation
+
+ !-------------------------------------------------------------
+
+ subroutine strain_tensor(mesh, strain11, strain22, strain12, uVelocity, vVelocity, normalVectorPolygon)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ strain11, &
+ strain22, &
+ strain12
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uVelocity, &
+ vVelocity
+
+ real(kind=RKIND), dimension(:,:,:), intent(in) :: &
+ normalVectorPolygon
+
+
+ integer :: &
+ iCell, &
+ iEdgeOnCell, &
+ iEdge, &
+ iVertexOnEdge, &
+ iVertex
+
+ real(kind=RKIND) :: &
+ uVelocityEdge, &
+ vVelocityEdge
+
+ do iCell = 1, mesh % nCells
+
+ strain11(iCell) = 0.0_RKIND
+ strain22(iCell) = 0.0_RKIND
+ strain12(iCell) = 0.0_RKIND
+
+ do iEdgeOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iEdge = mesh % edgesOnCell % array(iEdgeOnCell,iCell)
+
+ uVelocityEdge = 0.0_RKIND
+ vVelocityEdge = 0.0_RKIND
+
+ do iVertexOnEdge = 1, 2
+
+ iVertex = mesh % verticesOnEdge % array(iVertexOnEdge,iEdge)
+
+ uVelocityEdge = uVelocityEdge + uVelocity(iVertex)
+ vVelocityEdge = vVelocityEdge + vVelocity(iVertex)
+
+ enddo ! iVertexOnEdge
+
+ uVelocityEdge = uVelocityEdge / 2.0_RKIND
+ vVelocityEdge = vVelocityEdge / 2.0_RKIND
+
+ strain11(iCell) = strain11(iCell) + uVelocityEdge * normalVectorPolygon(1,iEdgeOnCell,iCell) * mesh % dvEdge % array(iEdge)
+ strain22(iCell) = strain22(iCell) + vVelocityEdge * normalVectorPolygon(2,iEdgeOnCell,iCell) * mesh % dvEdge % array(iEdge)
+ strain12(iCell) = strain12(iCell) + 0.5_RKIND * ( &
+ uVelocityEdge * normalVectorPolygon(2,iEdgeOnCell,iCell) + &
+ vVelocityEdge * normalVectorPolygon(1,iEdgeOnCell,iCell) ) * mesh % dvEdge % array(iEdge)
+
+ !write(*,*) iCell, uVelocityEdge, normalVectorPolygon(1,iEdgeOnCell,iCell), mesh % dvEdge % array(iEdge), &
+ ! uVelocityEdge * normalVectorPolygon(1,iEdgeOnCell,iCell) * mesh % dvEdge % array(iEdge)
+
+
+ enddo ! iEdgeOnCell
+
+ !write(*,*) iCell, strain11(iCell), mesh % areaCell % array(iCell), strain11(iCell) / mesh % areaCell % array(iCell)
+
+ strain11(iCell) = strain11(iCell) / mesh % areaCell % array(iCell)
+ strain22(iCell) = strain22(iCell) / mesh % areaCell % array(iCell)
+ strain12(iCell) = strain12(iCell) / mesh % areaCell % array(iCell)
+
+
+
+ enddo ! iCell
+
+ end subroutine strain_tensor
+
+ !-------------------------------------------------------------
+
+ subroutine stress_tensor(mesh, stress11, stress22, stress12, strain11, strain22, strain12, icePressure, dtDynamics, dtElastic)
+
+ use cice_testing, only: find_nearest_cell
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ stress11, &
+ stress22, &
+ stress12
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ strain11, &
+ strain22, &
+ strain12, &
+ icePressure
+
+ real(kind=RKIND), intent(in) :: &
+ dtDynamics, &
+ dtElastic
+
+ real(kind=RKIND) :: &
+ strainDivergence, &
+ strainTension, &
+ strainShearing, &
+ stress1, &
+ stress2
+
+ integer :: iCell
+
+ ! start testing!
+ integer :: iCellEx
+ real(kind=RKIND), parameter :: x0 = 1.2e6_RKIND
+ real(kind=RKIND), parameter :: y0 = 6.0e5_RKIND
+
+ ! stop testing!
+
+ do iCell = 1, mesh % nCells
+
+ strainDivergence = strain11(iCell) + strain22(iCell)
+ strainTension = strain11(iCell) - strain22(iCell)
+ strainShearing = strain12(iCell) * 2.0_RKIND
+
+ stress1 = stress11(iCell) + stress22(iCell)
+ stress2 = stress11(iCell) - stress22(iCell)
+
+ call evp_constitutive_relation(stress1, stress2, stress12(iCell), &
+ strainDivergence, strainTension, strainShearing, &
+ icePressure(iCell), mesh % areaCell % array(iCell), &
+ dtDynamics, dtElastic,iCell)
+
+ !call linear_constitutive_relation(stress1, stress2, stress12(iCell), &
+ ! strainDivergence, strainTension, strainShearing)
+
+ stress11(iCell) = 0.5_RKIND * (stress1 + stress2)
+ stress22(iCell) = 0.5_RKIND * (stress1 - stress2)
+
+ ! start testing!
+ !write(*,*) iCell, strain11(iCell), strain22(iCell), strain12(iCell), stress11(iCell), stress22(iCell), stress12(iCell)
+
+ !stress11(iCell) = strain11(iCell)
+ !stress22(iCell) = strain22(iCell)
+ !stress12(iCell) = strain12(iCell)
+ ! stop testing!
+
+ end do ! iCell
+
+ ! start testing!
+ !iCellEx = find_nearest_cell(mesh,x0,y0)
+ !write(53,*) "stress:", iCellEx, stress11(iCellEx), stress22(iCellEx), stress12(iCellEx), strain11(iCellEx), strain22(iCellEx), strain12(iCellEx)
+ ! stop testing!
+
+ end subroutine stress_tensor
+
+ !-------------------------------------------------------------
+
+ subroutine stress_divergence(mesh, stressDivergenceU, stressDivergenceV, stress11, stress22, stress12, normalVectorTriangle, interiorVertex)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ stressDivergenceU, &
+ stressDivergenceV
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ stress11, &
+ stress22, &
+ stress12
+
+ real(kind=RKIND), dimension(:,:,:), intent(in) :: &
+ normalVectorTriangle
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ real(kind=RKIND) :: &
+ stress11Edge, &
+ stress22Edge, &
+ stress12Edge
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iEdge, &
+ iCellOnEdge, &
+ iCell
+
+ do iVertex = 1, mesh % nVertices
+
+ if (interiorVertex(iVertex) == 1) then
+
+ stressDivergenceU(iVertex) = 0.0_RKIND
+ stressDivergenceV(iVertex) = 0.0_RKIND
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iEdge = mesh % edgesOnVertex % array(iVertexDegree,iVertex)
+
+ stress11Edge = 0.0_RKIND
+ stress22Edge = 0.0_RKIND
+ stress12Edge = 0.0_RKIND
+
+ do iCellOnEdge = 1, 2
+
+ iCell = mesh % cellsOnEdge % array(iCellOnEdge,iEdge)
+
+ stress11Edge = stress11Edge + stress11(iCell)
+ stress22Edge = stress22Edge + stress22(iCell)
+ stress12Edge = stress12Edge + stress12(iCell)
+
+ !write(*,*) stress11(iCell), stress22(iCell), stress12(iCell)
+
+ enddo ! iCellOnEdge
+
+ stress11Edge = stress11Edge / 2.0_RKIND
+ stress22Edge = stress22Edge / 2.0_RKIND
+ stress12Edge = stress12Edge / 2.0_RKIND
+
+ stressDivergenceU(iVertex) = stressDivergenceU(iVertex) + &
+ (stress11Edge * normalVectorTriangle(1,iVertexDegree,iVertex) + &
+ stress12Edge * normalVectorTriangle(2,iVertexDegree,iVertex)) * mesh % dcEdge % array(iEdge)
+
+ stressDivergenceV(iVertex) = stressDivergenceV(iVertex) + &
+ (stress22Edge * normalVectorTriangle(2,iVertexDegree,iVertex) + &
+ stress12Edge * normalVectorTriangle(1,iVertexDegree,iVertex)) * mesh % dcEdge % array(iEdge)
+
+ enddo ! iVertexDegree
+
+ stressDivergenceU(iVertex) = stressDivergenceU(iVertex) / mesh % areaTriangle % array(iVertex)
+ stressDivergenceV(iVertex) = stressDivergenceV(iVertex) / mesh % areaTriangle % array(iVertex)
+
+ !write(*,*) iVertex, stressDivergenceU(iVertex), stressDivergenceV(iVertex)
+
+ endif ! interiorVertex
+
+ enddo ! iVertex
+
+ end subroutine stress_divergence
+
+ !-------------------------------------------------------------
+
+ subroutine solve_velocity(nPoints, solvePoints, &
+ uVelocity, vVelocity, &
+ mass, coriolisParameter, &
+ stressDivergenceU, stressDivergenceV, &
+ airStressU, airStressV, &
+ surfaceTiltForceU, surfaceTiltForceV, &
+ oceanStressU, oceanStressV, &
+ oceanStressCoeff, dtElastic)
+
+ integer, intent(in) :: &
+ nPoints
+
+ integer, dimension(:), intent(in) :: &
+ solvePoints
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ uVelocity, &
+ vVelocity
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ mass, &
+ coriolisParameter, &
+ stressDivergenceU, &
+ stressDivergenceV, &
+ airStressU, &
+ airStressV, &
+ surfaceTiltForceU, &
+ surfaceTiltForceV, &
+ oceanStressU, &
+ oceanStressV, &
+ oceanStressCoeff
+
+ real(kind=RKIND), intent(in) :: &
+ dtElastic
+
+ real(kind=RKIND), dimension(2) :: &
+ rightHandSide
+
+ real(kind=RKIND), dimension(2,2) :: &
+ leftMatrix
+
+ real(kind=RKIND) :: &
+ solutionDenominator
+
+ integer :: &
+ iPoint
+
+ do iPoint = 1, nPoints
+
+ if (solvePoints(iPoint) == 1) then
+
+ ! U equation
+ leftMatrix(1,1) = mass(iPoint) / dtElastic + oceanStressCoeff(iPoint) * cosOceanTurningAngle
+ leftMatrix(1,2) = -mass(iPoint) * coriolisParameter(iPoint) - oceanStressCoeff(iPoint) * sinOceanTurningAngle
+
+ ! V equation
+ leftMatrix(2,1) = mass(iPoint) * coriolisParameter(iPoint) + oceanStressCoeff(iPoint) * sinOceanTurningAngle
+ leftMatrix(2,2) = mass(iPoint) / dtElastic + oceanStressCoeff(iPoint) * cosOceanTurningAngle
+
+ ! right hand side of matrix solve
+ rightHandSide(1) = stressDivergenceU(iPoint) + airStressU(iPoint) + surfaceTiltForceU(iPoint) + &
+ oceanStressCoeff(iPoint) * oceanStressU(iPoint) + (mass(iPoint) * uVelocity(iPoint)) / dtElastic
+
+ rightHandSide(2) = stressDivergenceV(iPoint) + airStressV(iPoint) + surfaceTiltForceV(iPoint) + &
+ oceanStressCoeff(iPoint) * oceanStressV(iPoint) + (mass(iPoint) * vVelocity(iPoint)) / dtElastic
+
+ ! solve the equation
+ solutionDenominator = leftMatrix(1,1) * leftMatrix(2,2) - leftMatrix(1,2) * leftMatrix(2,1)
+
+ uVelocity(iPoint) = (leftMatrix(2,2) * rightHandSide(1) - leftMatrix(1,2) * rightHandSide(2)) / solutionDenominator
+ vVelocity(iPoint) = (leftMatrix(1,1) * rightHandSide(2) - leftMatrix(2,1) * rightHandSide(1)) / solutionDenominator
+
+ !write(*,*) iPoint, leftMatrix(1,1), leftMatrix(1,2), leftMatrix(2,1), leftMatrix(2,2), rightHandSide(1), rightHandSide(2), uVelocity(iPoint), vVelocity(iPoint)
+ !write(*,*) stressDivergenceU(iPoint), airStressU(iPoint), surfaceTiltForceU(iPoint), &
+ ! oceanStressCoeff(iPoint) * oceanStressU(iPoint), mass(iPoint), uVelocity(iPoint), dtElastic
+ !write(*,*) iPoint, uVelocity(iPoint), stressDivergenceU(iPoint), airStressU(iPoint), mass(iPoint) / dtElastic
+
+ !write(*,*) iPoint, mass(iPoint), dtElastic
+
+
+ endif ! solvePoints
+
+ enddo ! iCell
+
+ !stop
+
+ end subroutine solve_velocity
+
+ !-------------------------------------------------------------
+
+ subroutine air_stress(airStressU, airStressV, uAirVelocity, vAirVelocity, iceAreaVertex)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ airStressU, &
+ airStressV
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uAirVelocity, &
+ vAirVelocity, &
+ iceAreaVertex
+
+ real(kind=RKIND) :: &
+ windSpeed
+
+ integer :: &
+ nPoints, &
+ iPoint
+
+ real(kind=RKIND), parameter :: &
+ airStressCoeff = 0.0012_RKIND, &
+ densityAir = 1.3_RKIND
+
+ nPoints = size(airStressU,1)
+
+ do iPoint = 1, nPoints
+
+ windSpeed = sqrt(uAirVelocity(iPoint)**2 + vAirVelocity(iPoint)**2)
+
+ airStressU(iPoint) = densityAir * windSpeed * airStressCoeff * uAirVelocity(iPoint) * iceAreaVertex(iPoint)
+ airStressV(iPoint) = densityAir * windSpeed * airStressCoeff * vAirVelocity(iPoint) * iceAreaVertex(iPoint)
+
+ enddo ! iPoint
+
+ end subroutine air_stress
+
+ !-------------------------------------------------------------
+
+ subroutine ocean_stress(oceanStressU, oceanStressV, &
+ uOceanVelocity, vOceanVelocity)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ oceanStressU, &
+ oceanStressV
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uOceanVelocity, &
+ vOceanVelocity
+
+ integer :: &
+ nPoints, &
+ iPoint
+
+ nPoints = size(oceanStressU,1)
+
+ do iPoint = 1, nPoints
+
+ oceanStressU(iPoint) = uOceanVelocity(iPoint) * cosOceanTurningAngle - vOceanVelocity(iPoint) * sinOceanTurningAngle
+ oceanStressV(iPoint) = uOceanVelocity(iPoint) * sinOceanTurningAngle + vOceanVelocity(iPoint) * cosOceanTurningAngle
+
+ enddo ! iPoint
+
+ end subroutine ocean_stress
+
+ !-------------------------------------------------------------
+
+ subroutine ocean_stress_coefficient(oceanStressCoeff, &
+ uOceanVelocity, vOceanVelocity, &
+ uVelocity, vVelocity, &
+ iceArea)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ oceanStressCoeff
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ uOceanVelocity, &
+ vOceanVelocity, &
+ uVelocity, &
+ vVelocity, &
+ iceArea
+
+ real(kind=RKIND), parameter :: &
+ densityWater = 1026.0_RKIND, &
+ iceOceanDragCoeff = 0.00536_RKIND
+
+ integer :: &
+ nPoints, &
+ iPoint
+
+ nPoints = size(oceanStressCoeff,1)
+
+ do iPoint = 1, nPoints
+
+ oceanStressCoeff(iPoint) = iceOceanDragCoeff * densityWater * iceArea(iPoint) * &
+ sqrt((uOceanVelocity(iPoint) - uVelocity(iPoint))**2 + &
+ (vOceanVelocity(iPoint) - vVelocity(iPoint))**2)
+
+ enddo ! iPoint
+
+ end subroutine ocean_stress_coefficient
+
+ !-------------------------------------------------------------
+
+ subroutine surface_tilt(surfaceTiltForceU, surfaceTiltForceV)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ surfaceTiltForceU, &
+ surfaceTiltForceV
+
+ integer :: &
+ nPoints, &
+ iPoint
+
+ nPoints = size(surfaceTiltForceU,1)
+
+ do iPoint = 1, nPoints
+
+ surfaceTiltForceU(iPoint) = 0.0_RKIND
+ surfaceTiltForceV(iPoint) = 0.0_RKIND
+
+ enddo ! iPoint
+
+ end subroutine surface_tilt
+
+ !-------------------------------------------------------------
+
+ subroutine principal_stresses(mesh, principalStress1, principalStress2, stress11, stress22, stress12, icePressure)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ principalStress1, &
+ principalStress2
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ stress11, &
+ stress22, &
+ stress12, &
+ icePressure
+
+ real(kind=RKIND) :: &
+ sqrtContents
+
+ integer :: &
+ iCell
+
+ open(55,file='principal_stresses.txt')
+
+ do iCell = 1, mesh % nCells
+
+ sqrtContents = (stress11(iCell) + stress22(iCell))**2 - &
+ 4.0_RKIND * stress11(iCell) * stress22(iCell) + &
+ 4.0_RKIND * stress12(iCell)**2
+
+ principalStress1(iCell) = 0.5_RKIND * (stress11(iCell) + stress22(iCell)) + 0.5_RKIND * sqrt(sqrtContents)
+ principalStress2(iCell) = 0.5_RKIND * (stress11(iCell) + stress22(iCell)) - 0.5_RKIND * sqrt(sqrtContents)
+
+ write(55,*) principalStress1(iCell) / icePressure(iCell), principalStress2(iCell) / icePressure(iCell)
+
+ enddo ! iCell
+
+ close(55)
+
+ end subroutine principal_stresses
+
+ !-------------------------------------------------------------
+ ! Testing
+ !-------------------------------------------------------------
+
+ subroutine linear_constitutive_relation(stress1, stress2, stress12, &
+ strainDivergence, strainTension, strainShearing)
+
+ use cice_testing, only: lambda
+
+ real(kind=RKIND), intent(out) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND), intent(in) :: &
+ strainDivergence, &
+ strainTension, &
+ strainShearing
+
+ stress1 = lambda * strainDivergence
+ stress2 = lambda * strainTension
+ stress12 = lambda * strainShearing * 0.5_RKIND
+
+ end subroutine linear_constitutive_relation
+
+ !-------------------------------------------------------------
+
+ subroutine test_stress_divergence(mesh, icestate, weak, boundary)
+
+ use cice_testing, only: divergence_stress_test_stress_set_weak, divergence_stress_test_velocity_set
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(weak_type), pointer :: weak
+ type(boundary_type), pointer :: boundary
+
+ real(kind=RKIND) :: &
+ dtDynamics, dtElastic
+
+ call divergence_stress_test_velocity_set(&
+ weak % uVelocity % array, &
+ weak % vVelocity % array, &
+ mesh % xVertex % array, &
+ mesh % yVertex % array, &
+ "div1")
+
+ write(*,*) "strain test"
+ call strain_tensor(mesh, &
+ weak % strain11 % array, &
+ weak % strain22 % array, &
+ weak % strain12 % array, &
+ weak % uVelocity % array, &
+ weak % vVelocity % array, &
+ weak % normalVectorPolygon % array)
+
+ call stress_tensor(mesh, &
+ weak % stress11 % array, &
+ weak % stress22 % array, &
+ weak % stress12 % array, &
+ weak % strain11 % array, &
+ weak % strain22 % array, &
+ weak % strain12 % array, &
+ weak % icePressure % array, &
+ dtDynamics, dtElastic)
+
+ write(*,*) "divergence test"
+ !call divergence_stress_test_stress_set_weak(mesh, &
+ ! weak % stress11 % array, &
+ ! weak % stress22 % array, &
+ ! weak % stress12 % array)
+
+ call stress_divergence(mesh, &
+ weak % stressDivergenceU % array, &
+ weak % stressDivergenceV % array, &
+ weak % stress11 % array, &
+ weak % stress22 % array, &
+ weak % stress12 % array, &
+ weak % normalVectorTriangle % array, &
+ boundary % interiorVertex % array)
+
+ end subroutine test_stress_divergence
+
+ !-------------------------------------------------------------
+ ! Legacy
+ !-------------------------------------------------------------
+
+ subroutine boundary_cells(mesh, boundaryCell, boundaryCell2)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:), intent(out) :: &
+ boundaryCell, &
+ boundaryCell2
+
+ integer :: &
+ iCell, &
+ iCellOnCell, &
+ iAdjacentCell
+
+ do iCell = 1, mesh % nCells
+
+ boundaryCell(iCell) = 0
+
+ do iCellOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iAdjacentCell = mesh % cellsOnCell % array(iCellOnCell,iCell)
+
+ if (iAdjacentCell > mesh % nCells) then
+
+ boundaryCell(iCell) = 1
+ exit
+
+ endif
+
+ enddo ! iCellOnCell
+
+ enddo ! iCell
+
+ do iCell = 1, mesh % nCells
+
+ boundaryCell2(iCell) = 0
+
+ if (boundaryCell(iCell) == 0) then
+
+ do iCellOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iAdjacentCell = mesh % cellsOnCell % array(iCellOnCell,iCell)
+
+ if (boundaryCell(iAdjacentCell) == 1) then
+
+ boundaryCell2(iCell) = 1
+
+ endif
+
+ enddo ! iCellOnCell
+
+ endif
+
+ enddo ! iCell
+
+ end subroutine boundary_cells
+
+ !-------------------------------------------------------------
+
+ subroutine boundary_vertices(mesh, boundaryVertex, interiorVertex, interiorVertex2, boundaryCell)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:), intent(out) :: &
+ boundaryVertex, & ! boundary vertices of boundary vertex triangles
+ interiorVertex, & ! all normal shaped interior vertices
+ interiorVertex2
+
+ integer, dimension(:), intent(in) :: &
+ boundaryCell
+
+ integer :: &
+ nInteriorAdjacentCells
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell
+
+ logical :: &
+ linterior
+
+ ! boundary vertices
+ do iVertex = 1, mesh % nVertices
+
+ boundaryVertex(iVertex) = 0
+ interiorVertex(iVertex) = 0
+
+ nInteriorAdjacentCells = 0
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ if (iCell >= 1 .and. iCell <= mesh % nCells) then
+
+ nInteriorAdjacentCells = nInteriorAdjacentCells + 1
+
+ endif
+
+ enddo ! iVertexDegree
+
+ ! boundary vertices that have triangles
+ if (nInteriorAdjacentCells == 2) then
+
+ boundaryVertex(iVertex) = 1
+
+ endif
+
+ ! vertex points we directly calculate on
+ if (nInteriorAdjacentCells == 3) then
+
+ interiorVertex(iVertex) = 1
+
+ endif
+
+ enddo ! iVertex
+
+
+ ! second interior vertices
+ do iVertex = 1, mesh % nVertices
+
+ interiorVertex2(iVertex) = 0
+
+ linterior = .true.
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ if (boundaryCell(iCell) == 1) then
+
+ linterior = .false.
+
+ endif
+
+ enddo ! iVertexDegree
+
+ if (linterior) interiorVertex2(iVertex) = 1
+
+ enddo ! iVertex
+
+ end subroutine boundary_vertices
+
+ !-------------------------------------------------------------
+
+ subroutine boundary_edges(mesh, boundaryEdge)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:), intent(out) :: &
+ boundaryEdge
+
+ integer :: &
+ nInteriorAdjacentCells
+
+ integer :: &
+ iCell, &
+ iEdge, &
+ iCellsOnEdge
+
+ ! boundary edges
+ do iEdge = 1, mesh % nEdges
+
+ do iCellsOnEdge = 1, 2
+
+ iCell = mesh % cellsOnEdge % array(iCellsOnEdge,iEdge)
+
+ if (iCell > mesh % nCells) then
+ boundaryEdge(iEdge) = 1
+ else
+ boundaryEdge(iEdge) = 0
+ endif
+
+ enddo ! iCellsOnEdge
+
+ enddo ! iEdge
+
+ end subroutine boundary_edges
+
+ !-------------------------------------------------------------
+ ! LU decomposition
+ !-------------------------------------------------------------
+
+ subroutine solve_linear_basis_system(leftMatrix, rightHandSide, solutionVector)
+
+ real(kind=RKIND), dimension(3,3), intent(in) :: &
+ leftMatrix
+ real(kind=RKIND), dimension(3), intent(in) :: &
+ rightHandSide
+ real(kind=RKIND), dimension(3), intent(out) :: &
+ solutionVector
+
+ real(kind=RKIND), dimension(3,3) :: a
+ real(kind=RKIND), dimension(3) :: b
+ integer, dimension(3) :: indx
+ real(kind=RKIND) :: d
+
+ a(:,:) = leftMatrix(:,:)
+ b(:) = rightHandSide(:)
+
+ call ludcmp(a,indx,d)
+
+ call lubksb(a,indx,b)
+
+ solutionVector(:) = b(:)
+
+ end subroutine solve_linear_basis_system
+
+ !-------------------------------------------------------------
+
+ SUBROUTINE ludcmp(a,indx,d)
+
+ IMPLICIT NONE
+ REAL(kind=RKIND), DIMENSION(:,:), INTENT(INOUT) :: a
+ INTEGER, DIMENSION(:), INTENT(OUT) :: indx
+ REAL(kind=RKIND), INTENT(OUT) :: d
+ REAL(kind=RKIND), DIMENSION(size(a,1)) :: vv
+ REAL(kind=RKIND), PARAMETER :: TINY=1.0e-20_RKIND
+ INTEGER :: j,n,imax
+ n=assert_eq(size(a,1),size(a,2),size(indx),'ludcmp')
+ d=1.0
+ vv=maxval(abs(a),dim=2)
+ if (any(vv == 0.0)) call nrerror('singular matrix in ludcmp')
+ vv=1.0_RKIND/vv
+ do j=1,n
+ imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j)))
+ if (j /= imax) then
+ call swap(a(imax,:),a(j,:))
+ d=-d
+ vv(imax)=vv(j)
+ end if
+ indx(j)=imax
+ if (a(j,j) == 0.0) a(j,j)=TINY
+ a(j+1:n,j)=a(j+1:n,j)/a(j,j)
+ a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n))
+ end do
+ END SUBROUTINE ludcmp
+
+ !-------------------------------------------------------------
+
+ SUBROUTINE lubksb(a,indx,b)
+
+ IMPLICIT NONE
+ REAL(kind=RKIND), DIMENSION(:,:), INTENT(IN) :: a
+ INTEGER, DIMENSION(:), INTENT(IN) :: indx
+ REAL(kind=RKIND), DIMENSION(:), INTENT(INOUT) :: b
+ INTEGER :: i,n,ii,ll
+ REAL(kind=RKIND) :: summ
+ n=assert_eq(size(a,1),size(a,2),size(indx),'lubksb')
+ ii=0
+ do i=1,n
+ ll=indx(i)
+ summ=b(ll)
+ b(ll)=b(i)
+ if (ii /= 0) then
+ summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1))
+ else if (summ /= 0.0) then
+ ii=i
+ end if
+ b(i)=summ
+ end do
+ do i=n,1,-1
+ b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i)
+ end do
+ END SUBROUTINE lubksb
+
+ !-------------------------------------------------------------
+
+ FUNCTION assert_eq(n1,n2,n3,string)
+ CHARACTER(LEN=*), INTENT(IN) :: string
+ INTEGER, INTENT(IN) :: n1,n2,n3
+ INTEGER :: assert_eq
+ if (n1 == n2 .and. n2 == n3) then
+ assert_eq=n1
+ else
+ write (*,*) 'nrerror: an assert_eq failed with this tag:', &
+ string
+ STOP 'program terminated by assert_eq3'
+ end if
+ END FUNCTION assert_eq
+
+ !-------------------------------------------------------------
+
+ FUNCTION imaxloc(arr)
+ REAL(kind=RKIND), DIMENSION(:), INTENT(IN) :: arr
+ INTEGER :: imaxloc
+ INTEGER, DIMENSION(1) :: imax
+ imax=maxloc(arr(:))
+ imaxloc=imax(1)
+ END FUNCTION imaxloc
+
+ !-------------------------------------------------------------
+
+ SUBROUTINE nrerror(string)
+ CHARACTER(LEN=*), INTENT(IN) :: string
+ write (*,*) 'nrerror: ',string
+ STOP 'program terminated by nrerror'
+ END SUBROUTINE nrerror
+
+ !-------------------------------------------------------------
+
+ FUNCTION outerprod(a,b)
+ REAL(kind=RKIND), DIMENSION(:), INTENT(IN) :: a,b
+ REAL(kind=RKIND), DIMENSION(size(a),size(b)) :: outerprod
+ outerprod = spread(a,dim=2,ncopies=size(b)) * &
+ spread(b,dim=1,ncopies=size(a))
+ END FUNCTION outerprod
+
+ !-------------------------------------------------------------
+
+ SUBROUTINE swap(a,b)
+ REAL(kind=RKIND), DIMENSION(:), INTENT(INOUT) :: a,b
+ REAL(kind=RKIND), DIMENSION(SIZE(a)) :: dum
+ dum=a
+ a=b
+ b=dum
+ END SUBROUTINE swap
+
+ !-------------------------------------------------------------
+
+end module cice_dynamics_shared
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_init.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_init.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_init.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,20 +0,0 @@
-module cice_init
-
- implicit none
-
-
-contains
-
- !-------------------------------------------------------------------------
-
- subroutine init_data()
-
- use cice_state
-
- call init_cice_state()
-
- end subroutine init_data
-
- !-------------------------------------------------------------------------
-
-end module cice_init
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_kinds_mod.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_kinds_mod.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_kinds_mod.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,42 +0,0 @@
-!=======================================================================
-!BOP
-!
-! !MODULE: ice_kinds_mod - defines variable precision
-!
-! !DESCRIPTION:
-!
-! Defines variable precision for all common data types \\
-! Code originally based on kinds_mod.F in POP
-!
-! !REVISION HISTORY:
-! SVN:$Id: ice_kinds_mod.F90 140 2008-07-25 20:15:53Z eclare $
-!
-! author: Elizabeth C. Hunke and William H. Lipscomb, LANL
-! 2006: ECH converted to free source form (F90)
-!
-! !INTERFACE:
-!
- module cice_kinds_mod
-!
-! !USES:
-!
-!EOP
-!=======================================================================
-
- use mpas_kind_types
-
- implicit none
-
- integer, parameter :: char_len = 80, &
- char_len_long = 256, &
- log_kind = kind(.true.), &
- int_kind = kind(1), &
- real_kind = RKIND, &
- dbl_kind = RKIND, &
- r16_kind = selected_real_kind(26)
-
-!=======================================================================
-
- end module cice_kinds_mod
-
-!=======================================================================
Modified: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_mpas_core.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_mpas_core.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_mpas_core.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -14,14 +14,13 @@
contains
- !-------------------------------------------------------------------------
+!--------------------------------------------------------------------------
subroutine mpas_core_init(domain, startTimeStamp)
use mpas_configure
use mpas_grid_types
- use cice_init
-
+
implicit none
type (domain_type), intent(inout) :: domain
@@ -46,12 +45,9 @@
current_outfile_frames = 0
- ! CICE init
- call init_data()
-
end subroutine mpas_core_init
- !-------------------------------------------------------------------------
+!--------------------------------------------------------------------------
subroutine simulation_clock_init(domain, dt, startTimeStamp)
@@ -102,35 +98,39 @@
end subroutine simulation_clock_init
- !-------------------------------------------------------------------------
+!--------------------------------------------------------------------------
subroutine mpas_init_block(block, mesh, dt)
use mpas_grid_types
use mpas_rbf_interpolation
use mpas_vector_reconstruction
-
+
+ use cice_dynamics_shared, only: init_weak_dynamics, &
+ test_stress_divergence
+
+ use cice_testing, only: init_square_test_case_weak
+
implicit none
-
+
type (block_type), intent(inout) :: block
type (mesh_type), intent(inout) :: mesh
real (kind=RKIND), intent(in) :: dt
-
- call compute_mesh_scaling(mesh)
+ integer :: iCell
+
call mpas_rbf_interp_initialize(mesh)
call mpas_init_reconstruct(mesh)
- !call mpas_reconstruct(mesh, block % state % time_levs(1) % state % u % array, &
- ! block % state % time_levs(1) % state % uReconstructX % array, &
- ! block % state % time_levs(1) % state % uReconstructY % array, &
- ! block % state % time_levs(1) % state % uReconstructZ % array, &
- ! block % state % time_levs(1) % state % uReconstructZonal % array, &
- ! block % state % time_levs(1) % state % uReconstructMeridional % array &
- ! )
+ ! dynamics
+ call init_weak_dynamics(mesh, block % weak, block % boundary, dt)
+
+ ! test case
+ call init_square_test_case_weak(mesh, block % icestate, block % weak)
+
end subroutine mpas_init_block
- !-------------------------------------------------------------------------
+!--------------------------------------------------------------------------
subroutine mpas_core_run(domain, output_obj, output_frame)
@@ -138,7 +138,6 @@
use mpas_kind_types
use mpas_io_output
use mpas_timer
- use cice_timestep
implicit none
@@ -211,8 +210,8 @@
end subroutine mpas_core_run
- !-------------------------------------------------------------------------
-
+!--------------------------------------------------------------------------
+
subroutine write_output_frame(output_obj, output_frame, domain)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Compute diagnostic fields for a domain and write model state to output file
@@ -254,7 +253,7 @@
end subroutine write_output_frame
- !-------------------------------------------------------------------------
+!--------------------------------------------------------------------------
subroutine compute_output_diagnostics(state, grid)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -277,9 +276,36 @@
integer :: iEdge, k
end subroutine compute_output_diagnostics
+
+!--------------------------------------------------------------------------
+
+ subroutine mpas_timestep(domain, itimestep, dt, timeStamp)
+
+ use mpas_grid_types
+ use cice_time_integration
+ use mpas_timer
+
+ implicit none
+
+ type (domain_type), intent(inout) :: domain
+ integer, intent(in) :: itimestep
+ real (kind=RKIND), intent(in) :: dt
+ character(len=*), intent(in) :: timeStamp
- !-------------------------------------------------------------------------
+ type (block_type), pointer :: block
+
+ write(*,*) trim(timeStamp)
+
+ block => domain % blocklist
+ do while (associated(block))
+ call cice_timestep(block, dt)
+ block => block % next
+ end do
+
+ end subroutine mpas_timestep
+!--------------------------------------------------------------------------
+
subroutine mpas_core_finalize(domain)
use mpas_grid_types
@@ -293,39 +319,6 @@
end subroutine mpas_core_finalize
- !-------------------------------------------------------------------------
+!--------------------------------------------------------------------------
- subroutine compute_mesh_scaling(mesh)
-
- use mpas_grid_types
-
- implicit none
-
- type (mesh_type), intent(inout) :: mesh
-
- integer :: iEdge, cell1, cell2
- real (kind=RKIND), dimension(:), pointer :: meshDensity, meshScalingDel2, meshScalingDel4
-
- meshDensity => mesh % meshDensity % array
- meshScalingDel2 => mesh % meshScalingDel2 % array
- meshScalingDel4 => mesh % meshScalingDel4 % array
-
- !
- ! Compute the scaling factors to be used in the del2 and del4 dissipation
- !
- !meshScalingDel2(:) = 1.0
- !meshScalingDel4(:) = 1.0
- !if (config_h_ScaleWithMesh) then
- ! do iEdge=1,mesh%nEdges
- ! cell1 = mesh % cellsOnEdge % array(1,iEdge)
- ! cell2 = mesh % cellsOnEdge % array(2,iEdge)
- ! meshScalingDel2(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/12.0)
- ! meshScalingDel4(iEdge) = 1.0 / ( (meshDensity(cell1) + meshDensity(cell2) )/2.0)**(5.0/6.0)
- ! end do
- !end if
-
- end subroutine compute_mesh_scaling
-
- !-------------------------------------------------------------------------
-
end module mpas_core
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_state.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_state.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_state.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,263 +0,0 @@
-!=======================================================================
-!BOP
-!
-! !MODULE: ice_state - primary state variables
-!
-! !DESCRIPTION:
-!
-! Primary state variables in various configurations
-! Note: other state variables are at the end of this...
-! The primary state variable names are:
-!-------------------------------------------------------------------
-! for each category aggregated over units
-! categories
-!-------------------------------------------------------------------
-! aicen(i,j,n) aice(i,j) ---
-! vicen(i,j,n) vice(i,j) m
-! vsnon(i,j,n) vsno(i,j) m
-! trcrn(i,j,it,n) trcr(i,j,it)
-!
-! Area is dimensionless because aice is the fractional area
-! (normalized so that the sum over all categories, including open
-! water, is 1.0). That is why vice/vsno have units of m instead of m^3.
-!
-! Variable names follow these rules:
-!
-! (1) For 3D variables (indices i,j,n), write 'ice' or 'sno' or
-! 'sfc' and put an 'n' at the end.
-! (2) For 2D variables (indices i,j) aggregated over all categories,
-! write 'ice' or 'sno' or 'sfc' without the 'n'.
-! (3) For 2D variables (indices i,j) associated with an individual
-! category, write 'i' or 's' instead of 'ice' or 'sno' and put an 'n'
-! at the end: e.g. hin, hsn. These are not declared here
-! but in individual modules (e.g., ice_therm_vertical).
-!
-! !REVISION HISTORY:
-! SVN:$Id: ice_state.F90 567 2013-01-07 02:57:36Z eclare $
-!
-! authors C. M. Bitz, UW
-! Elizabeth C. Hunke and William H. Lipscomb, LANL
-!
-! 2004: Block structure added by William Lipscomb
-! 2006: Converted to free form source (F90) by Elizabeth Hunke
-!
-! !INTERFACE:
-!
- module cice_state
-!
-! !USES:
-!
- use cice_kinds_mod
- use cice_domain_size
- !use ice_blocks
-!
-!EOP
-!
- implicit none
- save
-
- !-----------------------------------------------------------------
- ! state of the ice aggregated over all categories
- !-----------------------------------------------------------------
-
- !real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
- ! aice , & ! concentration of ice
- ! vice , & ! volume per unit area of ice (m)
- ! vsno ! volume per unit area of snow (m)
-
- !real (kind=dbl_kind), &
- ! dimension(nx_block,ny_block,max_ntrcr,max_blocks) :: &
- ! trcr ! ice tracers
- ! 1: surface temperature of ice/snow (C)
-
- !-----------------------------------------------------------------
- ! state of the ice for each category
- !-----------------------------------------------------------------
-
- !real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks):: &
- ! aice0 ! concentration of open water
-
- !real (kind=dbl_kind), &
- ! dimension (nx_block,ny_block,ncat,max_blocks) :: &
- ! aicen , & ! concentration of ice
- ! vicen , & ! volume per unit area of ice (m)
- ! vsnon ! volume per unit area of snow (m)
-
- !real (kind=dbl_kind), &
- ! dimension (nx_block,ny_block,max_ntrcr,ncat,max_blocks) :: &
- ! trcrn ! tracers
- ! 1: surface temperature of ice/snow (C)
-
- !integer (kind=int_kind), dimension (max_ntrcr) :: &
- ! trcr_depend ! = 0 for ice area tracers
- ! = 1 for ice volume tracers
- ! = 2 for snow volume tracers
-
- integer (kind=int_kind) :: &
- ntrcr ! number of tracers in use
-
- !-----------------------------------------------------------------
- ! indices and flags for tracers
- !-----------------------------------------------------------------
-
- integer (kind=int_kind) :: &
- ntraceb , & ! number of bio layer tracers in use &
- ntrace_start ! index of first bio tracer
-
- integer (kind=int_kind) :: &
- nt_Tsfc , & ! ice/snow temperature
- nt_qice , & ! volume-weighted ice enthalpy (in layers)
- nt_qsno , & ! volume-weighted snow enthalpy (in layers)
- nt_sice , & ! volume-weighted ice bulk salinity (CICE grid layers)
- nt_fbri , & ! volume fraction of ice with dynamic salt (hinS/vicen*aicen)
- nt_iage , & ! volume-weighted ice age
- nt_FY , & ! area-weighted first-year ice area
- nt_alvl , & ! level ice area fraction
- nt_vlvl , & ! level ice volume fraction
- nt_apnd , & ! melt pond area fraction
- nt_hpnd , & ! melt pond depth
- nt_ipnd , & ! melt pond refrozen lid thickness
- nt_aero , & ! starting index for aerosols in ice
- nt_bgc_N_sk, & ! algae (skeletal layer)
- nt_bgc_C_sk, & !
- nt_bgc_chl_sk, & !
- nt_bgc_Nit_sk, & ! nutrients (skeletal layer)
- nt_bgc_Am_sk, & !
- nt_bgc_Sil_sk, & !
- nt_bgc_DMSPp_sk, & ! trace gases (skeletal layer)
- nt_bgc_DMSPd_sk, & !
- nt_bgc_DMS_sk, & !
- nt_bgc_Nit_ml, & ! nutrients (ocean mixed layer)
- nt_bgc_Am_ml, & !
- nt_bgc_Sil_ml, & !
- nt_bgc_DMSP_ml, & ! trace gases (ocean mixed layer)
- nt_bgc_DMS_ml, & !
- nt_bgc_N, & ! algae: bulk quantities are tracers (volume preserved)
- nt_bgc_C, & !
- nt_bgc_chl, & !
- nt_bgc_NO, & ! nutrients
- nt_bgc_NH, & !
- nt_bgc_Sil, & !
- nt_bgc_DMSPp, & ! trace gases (skeletal layer)
- nt_bgc_DMSPd, & !
- nt_bgc_DMS, & !
- nt_bgc_PON, & ! zooplankton and detritus
- nt_bgc_S ! Bulk salinity in fraction ice with dynamic salinity (Bio grid)
-
- logical (kind=log_kind) :: &
- tr_iage, & ! if .true., use age tracer
- tr_FY, & ! if .true., use first-year area tracer
- tr_lvl, & ! if .true., use level ice tracer
- tr_pond, & ! if .true., use melt pond tracer
- tr_pond_cesm,& ! if .true., use cesm pond tracer
- tr_pond_lvl, & ! if .true., use level-ice pond tracer
- tr_pond_topo,& ! if .true., use explicit topography-based ponds
- tr_aero ! if .true., use aerosol tracers
-
- !-----------------------------------------------------------------
- ! dynamic variables closely related to the state of the ice
- !-----------------------------------------------------------------
-
- !real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
- ! uvel , & ! x-component of velocity (m/s)
- ! vvel , & ! y-component of velocity (m/s)
- ! divu , & ! strain rate I component, velocity divergence (1/s)
- ! shear , & ! strain rate II component (1/s)
- ! strength ! ice strength (N/m)
-
- !-----------------------------------------------------------------
- ! ice state at start of time step, saved for later in the step
- !-----------------------------------------------------------------
-
- !real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
- ! aice_init ! initial concentration of ice, for diagnostics
-
- !real (kind=dbl_kind), &
- ! dimension(nx_block,ny_block,ncat,max_blocks) :: &
- ! aicen_init , & ! initial ice concentration, for linear ITD
- ! vicen_init ! initial ice volume (m), for linear ITD
-
- ! logical (kind=log_kind) :: &
- ! hbrine ! if .true., brine height differs from ice thickness
-
-!=======================================================================
-
- contains
-
-!=======================================================================
-
- subroutine init_cice_state()
-
- tr_iage = .false.
- tr_FY = .false.
- tr_lvl = .false.
- tr_pond = .false.
- tr_pond_cesm = .false.
- tr_pond_lvl = .false.
- tr_pond_topo = .false.
- tr_aero = .false.
-
- nt_Tsfc = 1 ! index tracers, starting with Tsfc = 1
- ntrcr = 1 ! count tracers, starting with Tsfc = 1
-
- nt_qice = ntrcr + 1
- ntrcr = ntrcr + nilyr ! qice in nilyr layers
-
- nt_qsno = ntrcr + 1
- ntrcr = ntrcr + nslyr ! qsno in nslyr layers
-
- nt_sice = ntrcr + 1
- ntrcr = ntrcr + nilyr ! sice in nilyr layers
-
- nt_iage = 0
- if (tr_iage) then
- nt_iage = ntrcr + 1 ! chronological ice age
- ntrcr = ntrcr + 1
- endif
-
- nt_FY = 0
- if (tr_FY) then
- nt_FY = ntrcr + 1 ! area of first year ice
- ntrcr = ntrcr + 1
- endif
-
- nt_alvl = 0
- nt_vlvl = 0
- if (tr_lvl) then
- nt_alvl = ntrcr + 1
- ntrcr = ntrcr + 1
- nt_vlvl = ntrcr + 1
- ntrcr = ntrcr + 1
- endif
-
- nt_apnd = 0
- nt_hpnd = 0
- nt_ipnd = 0
- if (tr_pond) then ! all explicit melt pond schemes
- nt_apnd = ntrcr + 1
- ntrcr = ntrcr + 1
- nt_hpnd = ntrcr + 1
- ntrcr = ntrcr + 1
- if (tr_pond_lvl) then
- nt_ipnd = ntrcr + 1 ! refrozen pond ice lid thickness
- ntrcr = ntrcr + 1 ! on level-ice ponds (if frzpnd='hlid')
- endif
- if (tr_pond_topo) then
- nt_ipnd = ntrcr + 1 ! refrozen pond ice lid thickness
- ntrcr = ntrcr + 1 !
- endif
- endif
-
- nt_aero = 0
- if (tr_aero) then
- nt_aero = ntrcr + 1
- ntrcr = ntrcr + 4*n_aero ! 4 dEdd layers, n_aero species
- endif
-
- end subroutine init_cice_state
-
-!=======================================================================
-
- end module cice_state
-
-!=======================================================================
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_test.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_test.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_test.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,262 +0,0 @@
-module cice_test
-
- use cice_kinds_mod
- use cice_domain_size
- use cice_constants
-
-contains
-
-!=======================================================================
-!BOP
-!
-! !ROUTINE: frzmlt_bottom_lateral - bottom and lateral heat fluxes
-!
-! !DESCRIPTION:
-!
-! Adjust frzmlt to account for changes in fhocn since from_coupler.
-! Compute heat flux to bottom surface.
-! Compute fraction of ice that melts laterally.
-!
-! !REVISION HISTORY:
-!
-! authors C. M. Bitz, UW
-! William H. Lipscomb, LANL
-! Elizabeth C. Hunke, LANL
-!
-! !INTERFACE:
-!
- subroutine frzmlt_bottom_lateral (nx_block, ny_block, &
- ilo, ihi, jlo, jhi, &
- ntrcr, dt, &
- aice, frzmlt, &
- vicen, vsnon, &
- trcrn, &
- sst, Tf, &
- strocnxT, strocnyT, &
- Tbot, fbot, &
- rside)
-!
-! !USES:
-!
-! !INPUT/OUTPUT PARAMETERS:
-
- integer (kind=int_kind), intent(in) :: &
- nx_block, ny_block, & ! block dimensions
- ilo,ihi,jlo,jhi , & ! beginning and end of physical domain
- ntrcr ! number of tracers
-
- real (kind=dbl_kind), intent(in) :: &
- dt ! time step
-
- real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
- aice , & ! ice concentration
- frzmlt , & ! freezing/melting potential (W/m^2)
- sst , & ! sea surface temperature (C)
- Tf , & ! freezing temperature (C)
- strocnxT, & ! ice-ocean stress, x-direction
- strocnyT ! ice-ocean stress, y-direction
-
- real (kind=dbl_kind), dimension(nx_block,ny_block,ncat), &
- intent(in) :: &
- vicen , & ! ice volume (m)
- vsnon ! snow volume (m)
-
- real (kind=dbl_kind), dimension(nx_block,ny_block,ntrcr,ncat), &
- intent(in) :: &
- trcrn ! tracer array
-
- real (kind=dbl_kind), dimension(nx_block,ny_block), &
- intent(out) :: &
- Tbot , & ! ice bottom surface temperature (deg C)
- fbot , & ! heat flux to ice bottom (W/m^2)
- rside ! fraction of ice that melts laterally
-!
-!EOP
-!
- integer (kind=int_kind) :: &
- i, j , & ! horizontal indices
- n , & ! thickness category index
- k , & ! layer index
- ij , & ! horizontal index, combines i and j loops
- imelt ! number of cells with ice melting
-
- integer (kind=int_kind), dimension (nx_block*ny_block) :: &
- indxi, indxj ! compressed indices for cells with ice melting
-
- real (kind=dbl_kind), dimension (:), allocatable :: &
- etot , & ! total energy in column
- fside ! lateral heat flux (W/m^2)
-
- real (kind=dbl_kind) :: &
- deltaT , & ! SST - Tbot >= 0
- ustar , & ! skin friction velocity for fbot (m/s)
- wlat , & ! lateral melt rate (m/s)
- xtmp ! temporary variable
-
- ! Parameters for bottom melting
-
- ! 0.006 = unitless param for basal heat flx ala McPhee and Maykut
-
- real (kind=dbl_kind), parameter :: &
- cpchr = -cp_ocn*rhow*0.006_dbl_kind
-
-
- ! Parameters for lateral melting
-
- real (kind=dbl_kind), parameter :: &
- floediam = 300.0_dbl_kind, & ! effective floe diameter (m)
- alpha = 0.66_dbl_kind , & ! constant from Steele (unitless)
- m1 = 1.6e-6_dbl_kind , & ! constant from Maykut & Perovich
- ! (m/s/deg^(-m2))
- m2 = 1.36_dbl_kind ! constant from Maykut & Perovich
- ! (unitless)
-
- do j = 1, ny_block
- do i = 1, nx_block
- rside(i,j) = c0
- Tbot (i,j) = Tf(i,j)
- fbot (i,j) = c0
- enddo
- enddo
-
- !-----------------------------------------------------------------
- ! Identify grid cells where ice can melt.
- !-----------------------------------------------------------------
-
- imelt = 0
- do j = jlo, jhi
- do i = ilo, ihi
-#if (defined notz_experiment || defined notz_fieldwork)
- if (aice(i,j) > puny) then
-#else
- if (aice(i,j) > puny .and. frzmlt(i,j) < c0) then ! ice can melt
-#endif
- imelt = imelt + 1
- indxi(imelt) = i
- indxj(imelt) = j
- endif
- enddo ! i
- enddo ! j
-
- allocate(etot (imelt))
- allocate(fside(imelt))
-
- do ij = 1, imelt ! cells where ice can melt
- i = indxi(ij)
- j = indxj(ij)
-
- fside(ij) = c0
-
- !-----------------------------------------------------------------
- ! Use boundary layer theory for fbot.
- ! See Maykut and McPhee (1995): JGR, 100, 24,691-24,703.
- !-----------------------------------------------------------------
-
- deltaT = max((sst(i,j)-Tbot(i,j)),c0)
-
- ! strocnx has units N/m^2 so strocnx/rho has units m^2/s^2
- ustar = sqrt (sqrt(strocnxT(i,j)**2+strocnyT(i,j)**2)/rhow)
- ustar = max (ustar,ustar_min)
-
-#if defined notz_experiment
- fbot(i,j) = c0
-!#elif (defined pond_barrow || defined snowice_maksym)
-! fbot(i,j) = c0
-#elif defined notz_fieldwork
- fbot(i,j) = cpchr * deltaT * ustar ! < 0
-#else
- fbot(i,j) = cpchr * deltaT * ustar ! < 0
- fbot(i,j) = max (fbot(i,j), frzmlt(i,j)) ! frzmlt < fbot < 0
-#endif
-
-!!! uncomment to use all frzmlt for standalone runs
- ! fbot(i,j) = min (c0, frzmlt(i,j))
-
- !-----------------------------------------------------------------
- ! Compute rside. See these references:
- ! Maykut and Perovich (1987): JGR, 92, 7032-7044
- ! Steele (1992): JGR, 97, 17,729-17,738
- !-----------------------------------------------------------------
-
- wlat = m1 * deltaT**m2 ! Maykut & Perovich
- rside(i,j) = wlat*dt*pi/(alpha*floediam) ! Steele
- rside(i,j) = max(c0,min(rside(i,j),c1))
-
- enddo ! ij
-
- !-----------------------------------------------------------------
- ! Compute heat flux associated with this value of rside.
- !-----------------------------------------------------------------
-
- do n = 1, ncat
-
- do ij = 1, imelt
- etot(ij) = c0
- enddo
-
- ! melting energy/unit area in each column, etot < 0
-
- do k = 1, nslyr
-!DIR$ CONCURRENT !Cray
-!cdir nodep !NEC
-!ocl novrec !Fujitsu
- do ij = 1, imelt
- i = indxi(ij)
- j = indxj(ij)
- etot(ij) = etot(ij) + trcrn(i,j,nt_qsno+k-1,n) &
- * vsnon(i,j,n)/real(nslyr,kind=dbl_kind)
- enddo ! ij
- enddo
-
- do k = 1, nilyr
-!DIR$ CONCURRENT !Cray
-!cdir nodep !NEC
-!ocl novrec !Fujitsu
- do ij = 1, imelt
- i = indxi(ij)
- j = indxj(ij)
- etot(ij) = etot(ij) + trcrn(i,j,nt_qice+k-1,n) &
- * vicen(i,j,n)/real(nilyr,kind=dbl_kind)
- enddo ! ij
- enddo ! nilyr
-
-!DIR$ CONCURRENT !Cray
-!cdir nodep !NEC
-!ocl novrec !Fujitsu
- do ij = 1, imelt
- i = indxi(ij)
- j = indxj(ij)
- ! lateral heat flux
- fside(ij) = fside(ij) + rside(i,j)*etot(ij)/dt ! fside < 0
- enddo ! ij
-
- enddo ! n
-
- !-----------------------------------------------------------------
- ! Limit bottom and lateral heat fluxes if necessary.
- !-----------------------------------------------------------------
-
-!DIR$ CONCURRENT !Cray
-!cdir nodep !NEC
-!ocl novrec !Fujitsu
- do ij = 1, imelt
- i = indxi(ij)
- j = indxj(ij)
-
- xtmp = frzmlt(i,j)/(fbot(i,j) + fside(ij) + puny)
- xtmp = min(xtmp, c1)
-#if !(defined notz_experiment || defined notz_fieldwork || defined snowice_maksym || defined pond_barrow)
- fbot(i,j) = fbot(i,j) * xtmp
-#endif
- rside(i,j) = rside(i,j) * xtmp
- enddo ! ij
-
- deallocate(etot)
- deallocate(fside)
-
- end subroutine frzmlt_bottom_lateral
-
-!=======================================================================
-
-end module cice_test
-
Added: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_testing.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_testing.F (rev 0)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_testing.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -0,0 +1,1876 @@
+module cice_testing
+
+ use mpas_grid_types
+
+ implicit none
+
+ private
+ public :: divergence_stress_test_velocity_set, &
+ divergence_stress_test_stress_set_hex, &
+ divergence_stress_test_stress_set_tri, &
+ divergence_stress_test_stress_set_weak, &
+ gnuplot_cell, &
+ gnuplot_triangle, &
+ matlab_jet, &
+ plot_boundary_triangles, &
+ boundary_locations, &
+ labels_verticesOnCell, &
+ labels_verticesDegreeOnVertex, &
+ writeout_minmax, &
+ writeout_on_cell_real, &
+ writeout_on_vertex_real, &
+ writeout_on_vertex_on_cell_real, &
+ writeout_array, &
+ writeout_edgecell, &
+ writeout_edgeedgecell, &
+ rotate_ninety, &
+ reverse_grids, &
+ init_square_test_case_weak, &
+ gnuplot_vertexvector, &
+ find_nearest_cell
+
+ integer :: &
+ iObject = 1, &
+ iLabel = 1, &
+ iArrow = 1
+
+ real(kind=RKIND) :: &
+ xmin = 1.0e30_RKIND, &
+ xmax = -1.0e30_RKIND, &
+ ymin = 1.0e30_RKIND, &
+ ymax = -1.0e30_RKIND
+
+ ! parameter in constant stress strain relation
+ real(kind=RKIND), parameter, public :: lambda = 1.0e9_RKIND
+
+contains
+
+ !-------------------------------------------------------------
+ ! Square test case
+ !-------------------------------------------------------------
+
+ subroutine init_square_test_case_weak(mesh, icestate, weak)
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(weak_type), pointer :: weak
+
+ call init_atmos_velocity(weak % uAirVelocity % array, &
+ weak % vAirVelocity % array, &
+ mesh % xVertex % array, &
+ mesh % yVertex % array, &
+ 0.0_RKIND)
+
+ call init_ocean_velocity(weak % uOceanVelocity % array, &
+ weak % vOceanVelocity % array, &
+ mesh % xVertex % array, &
+ mesh % yVertex % array)
+
+ call init_ice_state(mesh, icestate, weak)
+
+ end subroutine init_square_test_case_weak
+
+ !-------------------------------------------------------------
+
+ subroutine init_square_test_case_tri(mesh, icestate, forcing)
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(forcing_type), pointer :: forcing
+
+ call init_atmos_velocity(forcing % uAirVelocity % array, &
+ forcing % vAirVelocity % array, &
+ mesh % xCell % array, &
+ mesh % yCell % array, &
+ 0.0_RKIND)
+
+ call init_ocean_velocity(forcing % uOceanVelocity % array, &
+ forcing % vOceanVelocity % array, &
+ mesh % xCell % array, &
+ mesh % yCell % array)
+
+ !call init_ice_state(mesh, icestate)
+
+ end subroutine init_square_test_case_tri
+
+ !-------------------------------------------------------------
+
+ subroutine init_square_test_case_hex(mesh, icestate, hexfor)
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(hexfor_type), pointer :: hexfor
+
+ call init_atmos_velocity(hexfor % uAirVelocity % array, &
+ hexfor % vAirVelocity % array, &
+ mesh % xVertex % array, &
+ mesh % yVertex % array, &
+ 0.0_RKIND)
+
+ call init_ocean_velocity(hexfor % uOceanVelocity % array, &
+ hexfor % vOceanVelocity % array, &
+ mesh % xVertex % array, &
+ mesh % yVertex % array)
+
+ !call init_ice_state(mesh, icestate)
+
+ end subroutine init_square_test_case_hex
+
+ !-------------------------------------------------------------
+
+ subroutine init_ocean_velocity(uOceanVelocity, vOceanVelocity, x, y)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ uOceanVelocity, &
+ vOceanVelocity
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ x, y
+
+ real(kind=RKIND), parameter :: a = 0.1_RKIND
+
+ real(kind=RKIND), parameter :: &
+ Lx = 1.28e6_RKIND, &
+ Ly = 1.28e6_RKIND
+
+ integer :: iPoint, nPoints
+
+ nPoints = size(uOceanVelocity,1)
+
+ do iPoint = 1, nPoints
+
+ uOceanVelocity(iPoint) = a * ((2.0_RKIND * y(iPoint) - Ly) / Ly)
+
+ vOceanVelocity(iPoint) = -a * ((2.0_RKIND * x(iPoint) - Lx) / Lx)
+
+ enddo ! iCell
+
+ end subroutine init_ocean_velocity
+
+ !-------------------------------------------------------------
+
+ subroutine init_atmos_velocity(uAirVelocity, vAirVelocity, xin, yin, time)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ uAirVelocity, &
+ vAirVelocity
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ xin, yin
+
+ real(kind=RKIND), intent(in) :: &
+ time
+
+ real(kind=RKIND), parameter :: pi = &
+ 3.142_RKIND
+
+ real(kind=RKIND), parameter :: a = 5.0_RKIND
+ real(kind=RKIND), parameter :: b = 3.0_RKIND
+
+ real(kind=RKIND), parameter :: theta = 4.0_RKIND * 24.0_RKIND * 3600.0_RKIND
+
+ real(kind=RKIND), parameter :: &
+ Lx = 1.28e6_RKIND, &
+ Ly = 1.28e6_RKIND
+
+ real(kind=RKIND) :: &
+ x, y, &
+ xmin, xmax, ymin, ymax
+
+ integer :: iPoint, nPoints
+
+ xmin = minval(xin)
+ xmax = maxval(xin)
+ ymin = minval(yin)
+ ymax = maxval(yin)
+
+ nPoints = size(uAirVelocity,1)
+
+ do iPoint = 1, nPoints
+
+ !x = ((xin(iPoint) - xmin) / (xmax - xmin)) * Lx
+ !y = ((yin(iPoint) - ymin) / (ymax - ymin)) * Ly
+
+ x = xin(iPoint)
+ y = yin(iPoint)
+
+ uAirVelocity(iPoint) = &!a * (y / Ly)
+ a + (sin((2.0_RKIND * pi * time) / theta) - b) * sin(2.0_RKIND * pi * (x / Lx)) * sin(pi * (y / Ly))
+ vAirVelocity(iPoint) = &!0.0_RKIND!&
+ a + (sin((2.0_RKIND * pi * time) / theta) - b) * sin(2.0_RKIND * pi * (y / Ly)) * sin(pi * (x / Lx))
+
+ !uAirVelocity(iPoint) = a
+ !vAirVelocity(iPoint) = a
+
+ enddo ! iPoint
+
+ end subroutine init_atmos_velocity
+
+ !-------------------------------------------------------------
+
+ subroutine init_ice_state(mesh, icestate, weak)
+
+ type(mesh_type), intent(in) :: mesh
+
+ type(icestate_type), pointer :: icestate
+ type(weak_type), pointer :: weak
+
+ real(kind=RKIND) :: iceThickness
+
+ real(kind=RKIND) :: x
+
+ real(kind=RKIND), parameter :: &
+ Lx = 1.28e6_RKIND, &
+ rhoi = 917.0_RKIND
+
+ integer :: iCell
+
+ iceThickness = 2.0_RKIND
+
+ do iCell = 1, mesh % nCells
+
+ x = mesh % xCell % array(iCell)
+
+ icestate % iceAreaCell % array(iCell) = &!0.95_RKIND!&
+ max(min(x / Lx, 1.0_RKIND), 0.0_RKIND)
+
+ icestate % iceVolumeCell % array(iCell) = &
+ iceThickness * icestate % iceAreaCell % array(iCell)
+
+ icestate % totalMassCell % array(iCell) = icestate % iceVolumeCell % array(iCell) * rhoi
+
+ weak % uVelocity % array = weak % uOceanVelocity % array
+ weak % vVelocity % array = weak % vOceanVelocity % array
+
+ enddo ! iCell
+
+ end subroutine init_ice_state
+
+ !-------------------------------------------------------------
+ ! stress divergence operator test velocities
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_test_velocity_set(uVelocity, vVelocity, x, y, type)
+
+ real(kind=RKIND), dimension(:), intent(out) :: &
+ uVelocity, &
+ vVelocity
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ x, y
+
+ character(len=*), intent(in) :: &
+ type
+
+ integer :: iPoint, nPoints
+
+ real(kind=RKIND), parameter :: &
+ velocityConstantU = 112.87654_RKIND, &
+ velocityConstantV = -34.5678_RKIND, &
+ velocityScale = 1.0_RKIND
+
+ nPoints = size(uVelocity,1)
+
+ select case (type)
+ case ("zero")
+ write(*,*) "zero velocities"
+
+ ! zero velocities
+ do iPoint = 1, nPoints
+ uVelocity(iPoint) = 0.0_RKIND
+ vVelocity(iPoint) = 0.0_RKIND
+ enddo ! iPoint
+
+ case ("constant")
+ write(*,*) "constant velocities"
+
+ ! constant velocities
+ do iPoint = 1, nPoints
+ uVelocity(iPoint) = velocityConstantU
+ vVelocity(iPoint) = velocityConstantV
+ enddo ! iPoint
+
+ case ("linearx")
+ write(*,*) "linearx velocities"
+
+ ! constant velocities
+ do iPoint = 1, nPoints
+ uVelocity(iPoint) = x(iPoint)
+ vVelocity(iPoint) = 0.0_RKIND
+ enddo ! iPoint
+
+ case ("lineary")
+ write(*,*) "lineary velocities"
+
+ ! constant velocities
+ do iPoint = 1, nPoints
+ uVelocity(iPoint) = 0.0_RKIND
+ vVelocity(iPoint) = y(iPoint)
+ enddo ! iPoint
+
+ case ("constantsig12")
+ write(*,*) "constant sigma_12"
+
+ ! constant velocities
+ do iPoint = 1, nPoints
+ uVelocity(iPoint) = y(iPoint)
+ vVelocity(iPoint) = x(iPoint)
+ enddo ! iPoint
+
+ case ("div1")
+ write(*,*) "div1 velocities"
+
+ ! velocity field that gives constant divergence
+ do iPoint = 1, nPoints
+ uVelocity(iPoint) = x(iPoint)**2 / (4.0_RKIND * lambda) + (x(iPoint) * y(iPoint)) / lambda
+ vVelocity(iPoint) = y(iPoint)**2 / (4.0_RKIND * lambda) + (x(iPoint) * y(iPoint)) / lambda
+ uVelocity(iPoint) = uVelocity(iPoint) * velocityScale
+ vVelocity(iPoint) = vVelocity(iPoint) * velocityScale
+ enddo ! iPoint
+
+ end select
+
+ end subroutine divergence_stress_test_velocity_set
+
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_test_stress_set_hex(mesh, stress1, stress2, stress12)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(inout) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND) :: &
+ xpy, &
+ x, y
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ iVertex
+
+ do iCell = 1, mesh % nCells
+
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ x = mesh % xVertex % array(iVertex)
+ y = mesh % yVertex % array(iVertex)
+
+ xpy = x + y
+
+ ! divu = 1 ; divv = 1
+ stress1(iVertexOnCell,iCell) = 1.5_RKIND * xpy
+ stress2(iVertexOnCell,iCell) = -0.5_RKIND * xpy
+ stress12(iVertexOnCell,iCell) = 0.5_RKIND * xpy
+
+ ! divu = 1 ; divv = 0
+ !stress1(iVertexOnCell,iCell) = 0.5_RKIND * x
+ !stress2(iVertexOnCell,iCell) = 0.5_RKIND * x
+ !stress12(iVertexOnCell,iCell) = 0.5_RKIND * y
+
+ ! divu = 0 ; divv = 1
+ !stress1(iVertexOnCell,iCell) = 0.5_RKIND * y
+ !stress2(iVertexOnCell,iCell) = -0.5_RKIND * y
+ !stress12(iVertexOnCell,iCell) = 0.5_RKIND * x
+
+ ! others
+ !stress1(iVertexOnCell,iCell) = 0.5_RKIND * x
+ !stress2(iVertexOnCell,iCell) = 0.5_RKIND * x
+ !stress12(iVertexOnCell,iCell) = 0.0_RKIND
+
+ !stress1(iVertexOnCell,iCell) = 0.0_RKIND
+ !stress2(iVertexOnCell,iCell) = 0.0_RKIND
+ !stress12(iVertexOnCell,iCell) = xpy
+
+ !stress1(iVertexOnCell,iCell) = y
+ !stress2(iVertexOnCell,iCell) = -y
+ !stress12(iVertexOnCell,iCell) = 0.0_RKIND
+
+ stress1(iVertexOnCell,iCell) = 0.5_RKIND * x + y
+ stress2(iVertexOnCell,iCell) = 0.5_RKIND * y + x
+ stress12(iVertexOnCell,iCell) = 0.5_RKIND * (x + y)
+
+ stress1(iVertexOnCell,iCell) = x
+ stress2(iVertexOnCell,iCell) = 0.0_RKIND
+ stress12(iVertexOnCell,iCell) = 0.0_RKIND
+
+
+ enddo ! iVertexOnCell
+
+ enddo ! iCell
+
+ end subroutine divergence_stress_test_stress_set_hex
+
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_test_stress_set_tri(mesh, stress1, stress2, stress12)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(inout) :: &
+ stress1, &
+ stress2, &
+ stress12
+
+ real(kind=RKIND) :: &
+ xpy, &
+ x, y
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell
+
+ do iVertex = 1, mesh % nVertices
+
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+
+ xpy = x + y
+
+ ! divu = 1 ; divv = 1
+ !stress1(iVertexDegree,iVertex) = 1.5_RKIND * xpy
+ !stress2(iVertexDegree,iVertex) = -0.5_RKIND * xpy
+ !stress12(iVertexDegree,iVertex) = 0.5_RKIND * xpy
+
+ ! divu = 1 ; divv = 0
+ !stress1(iVertexDegree,iVertex) = 0.5_RKIND * x
+ !stress2(iVertexDegree,iVertex) = 0.5_RKIND * x
+ !stress12(iVertexDegree,iVertex) = 0.5_RKIND * y
+
+ ! divu = 0 ; divv = 1
+ !stress1(iVertexDegree,iVertex) = 0.5_RKIND * y
+ !stress2(iVertexDegree,iVertex) = -0.5_RKIND * y
+ !stress12(iVertexDegree,iVertex) = 0.5_RKIND * x
+
+ ! others
+ !stress1(iVertexDegree,iVertex) = 0.5_RKIND * x
+ !stress2(iVertexDegree,iVertex) = 0.5_RKIND * x
+ !stress12(iVertexDegree,iVertex) = 0.0_RKIND
+
+ !stress1(iVertexDegree,iVertex) = 0.0_RKIND
+ !stress2(iVertexDegree,iVertex) = 0.0_RKIND
+ !stress12(iVertexDegree,iVertex) = xpy
+
+ !stress1(iVertexDegree,iVertex) = y
+ !stress2(iVertexDegree,iVertex) = -y
+ !stress12(iVertexDegree,iVertex) = 0.0_RKIND
+
+ !stress1(iVertexDegree,iVertex) = 0.5_RKIND * x + y
+ !stress2(iVertexDegree,iVertex) = 0.5_RKIND * y + x
+ !stress12(iVertexDegree,iVertex) = 0.5_RKIND * (x + y)
+
+ stress1(iVertexDegree,iVertex) = x
+ stress2(iVertexDegree,iVertex) = 0.0_RKIND
+ stress12(iVertexDegree,iVertex) = 0.0_RKIND
+
+
+ enddo ! iVertexOnCell
+
+ enddo ! iCell
+
+ end subroutine divergence_stress_test_stress_set_tri
+
+ !-------------------------------------------------------------
+
+ subroutine divergence_stress_test_stress_set_weak(mesh, stress11, stress22, stress12)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(inout) :: &
+ stress11, &
+ stress22, &
+ stress12
+
+ real(kind=RKIND) :: &
+ xpy, &
+ x, y
+
+ integer :: &
+ iCell
+
+ do iCell = 1, mesh % nCells
+
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+
+ xpy = x + y
+
+ ! divu = 1 ; divv = 0
+ !stress11(iCell) = x
+ !stress22(iCell) = 0.0_RKIND
+ !stress12(iCell) = 0.0_RKIND
+
+ ! divu = 0 ; divv = 1
+ !stress11(iCell) = 0.0_RKIND
+ !stress22(iCell) = y
+ !stress12(iCell) = 0.0_RKIND
+
+ ! divu = 1 ; divv = 1
+ !stress11(iCell) = x
+ !stress22(iCell) = y
+ !stress12(iCell) = 0.0_RKIND
+
+ ! divu = 1 ; divv = 1
+ !stress11(iCell) = 0.5_RKIND * xpy
+ !stress22(iCell) = 0.5_RKIND * xpy
+ !stress12(iCell) = 0.5_RKIND * xpy
+
+ ! divu = 1 ; divv = 1
+ stress11(iCell) = 100.0_RKIND
+ stress22(iCell) = -1000.0_RKIND
+ stress12(iCell) = 1.0_RKIND * xpy
+
+ enddo ! iCell
+
+ end subroutine divergence_stress_test_stress_set_weak
+
+ !-------------------------------------------------------------
+ ! general plotting
+ !-------------------------------------------------------------
+
+ subroutine gnuplot_cell(mesh, cellArray, filename, append, nofill)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ cellArray
+
+ character(len=*), intent(in) :: &
+ filename
+
+ logical, optional, intent(in) :: &
+ append, &
+ nofill
+
+ real(kind=RKIND) :: &
+ x, y, x0, y0, xc, yc
+
+ real(kind=RKIND) :: &
+ fMin, fMax, &
+ fValue
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ iVertex
+
+ character(len=400) :: &
+ stroutvertex, &
+ strout, &
+ stroutlabel, &
+ stroutint
+
+ integer :: &
+ red, &
+ green, &
+ blue
+
+ character(len=7) :: &
+ color
+
+ ! optional arguments
+ logical :: &
+ lappend, &
+ lnofill
+
+ ! optional arguments
+ lappend = .false.
+ if (present(append)) lappend = append
+
+ lnofill = .false.
+ if (present(nofill)) lnofill = nofill
+
+ fMin = minval(cellArray(1:mesh % nCells))
+ fMax = maxval(cellArray(1:mesh % nCells))
+ !fMin = 1
+ !fMax = mesh % nCells
+
+ if (lappend) then
+ open(55,file=filename,action="write",position="append")
+ else
+ open(55,file=filename,action="write")
+ end if
+
+ open(56,file="border.txt",action="write")
+
+ write(56,*) minval(mesh % cellsOnCell % array), maxval(mesh % cellsOnCell % array)
+ write(56,*) minloc(mesh % cellsOnCell % array), maxloc(mesh % cellsOnCell % array)
+ write(56,*) size(mesh % cellsOnCell % array,1), size(mesh % cellsOnCell % array,2)
+
+ do iCell = 1, mesh % nCells
+
+ xc = mesh % xCell % array(iCell)
+ yc = mesh % yCell % array(iCell)
+
+ xmin = min(xmin,xc) ; xmax = max(xmax,xc)
+ ymin = min(ymin,yc) ; ymax = max(ymax,yc)
+
+ !write(stroutint,fmt='(i5)') iCell
+ !write(stroutlabel,fmt='(a,i5,a,a,a,f14.2,a,f14.2,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at ',xc, ",",yc, " center"
+ !iLabel = iLabel + 1
+ !write(55,*) trim(stroutlabel)
+
+ fValue = (cellArray(iCell) - fMin) / (fMax - fMin)
+ !fValue = (real(iCell,RKIND) - fMin) / (fMax - fMin)
+
+ write(strout,fmt='(a,i5,a)') "set object ",iObject," polygon from "
+
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ !write(*,*) iCell, iVertexOnCell
+
+ write(56,*) iCell, mesh % nCells, iVertexOnCell, mesh % cellsOnCell % array(iVertexOnCell, iCell)
+
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell, iCell)
+
+ x = mesh % xVertex % array(iVertex)
+ y = mesh % yVertex % array(iVertex)
+
+ xmin = min(xmin,x) ; xmax = max(xmax,x)
+ ymin = min(ymin,y) ; ymax = max(ymax,y)
+
+ if (iVertexOnCell == 1) then
+ x0 = x
+ y0 = y
+ endif
+
+ write(stroutvertex,fmt='(f14.2,a,f14.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ enddo ! iVertex
+
+ write(stroutvertex,fmt='(f14.2,a,f14.2)') x0, ",", y0
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ write(55,*) trim(strout)
+
+ color = matlab_jet(fValue)
+ !if (iCell == 3272) color = "#000000"
+
+ if (lnofill) then
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iObject,' fillstyle empty border lt -1'
+ else
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iObject,' fc rgb "', color, '" fillstyle solid'
+ endif
+
+ write(55,*) trim(strout)
+
+ iObject = iObject + 1
+
+ enddo ! iCell
+
+ write(stroutint,fmt='(e20.8)') fMax
+ write(stroutlabel,fmt='(a,i5,a,a,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at graph 1.02, 0.9'
+ iLabel = iLabel + 1
+ write(55,*) trim(stroutlabel)
+
+
+ write(stroutint,fmt='(e20.8)') fMin
+ write(stroutlabel,fmt='(a,i5,a,a,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at graph 1.02, 0.8'
+ iLabel = iLabel + 1
+ write(55,*) trim(stroutlabel)
+
+ close(55)
+ close(56)
+
+ !call gnuplot_triangle(mesh, mesh % xVertex % array, filename, .true., mesh % nCells, .true.)
+!stop
+ end subroutine gnuplot_cell
+
+ !-------------------------------------------------------------
+
+ subroutine gnuplot_triangle(mesh, vertexArray, filename, append, nofill)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ vertexArray
+
+ character(len=*), intent(in) :: &
+ filename
+
+ logical, optional, intent(in) :: &
+ append, &
+ nofill
+
+ real(kind=RKIND) :: &
+ x, y, x0, y0, xv, yv
+
+ real(kind=RKIND) :: &
+ fMin, fMax, &
+ fValue
+
+ integer :: &
+ iVertex, &
+ iVertexDegree, &
+ iCell
+
+ character(len=400) :: &
+ stroutvertex, &
+ strout, &
+ stroutlabel, &
+ stroutint
+
+ integer :: &
+ red, &
+ green, &
+ blue
+
+ character(len=7) :: &
+ color
+
+ logical :: atBoundary
+
+ ! optional arguments
+ logical :: &
+ lappend, &
+ lnofill
+
+ ! optional arguments
+ lappend = .false.
+ if (present(append)) lappend = append
+
+ lnofill = .false.
+ if (present(nofill)) lnofill = nofill
+
+ ! min/maxes
+ fMin = 1e30
+ fMax = -1e30
+
+ ! loop over triangles
+ do iVertex = 1, mesh % nVertices
+ atBoundary = .false.
+
+ ! loop over triangle vertices
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ if (iCell > mesh % nCells) then
+ atBoundary = .true.
+ end if
+ enddo ! iVertexDegree
+
+ if (.not. atBoundary) then
+ fMin = min(fMin,vertexArray(iVertex))
+ fMax = max(fMax,vertexArray(iVertex))
+ endif
+
+ enddo ! iVertex
+ !fMin = minval(vertexArray)
+ !fMax = maxval(vertexArray)
+ !fMin = 1
+ !fMax = mesh % nVertices
+
+ if (lappend) then
+ open(55,file=filename,action="write",position="append")
+ else
+ open(55,file=filename,action="write")
+ end if
+
+ ! loop over triangles
+ do iVertex = 1, mesh % nVertices
+
+ xv = mesh % xVertex % array(iVertex)
+ yv = mesh % yVertex % array(iVertex)
+
+ xmin = min(xmin,xv) ; xmax = max(xmax,xv)
+ ymin = min(ymin,yv) ; ymax = max(ymax,yv)
+
+ !write(stroutint,fmt='(i5)') iVertex
+ !write(stroutlabel,fmt='(a,i5,a,a,a,f14.2,a,f14.2,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at ',xv, ",",yv, " center"
+ !iLabel = iLabel + 1
+ !write(55,*) trim(stroutlabel)
+
+ atBoundary = .false.
+
+ fValue = (vertexArray(iVertex) - fMin) / (fMax - fMin)
+ !fValue = (real(iVertex,RKIND) - fMin) / (fMax - fMin)
+ !write(*,*) iVertex, vertexArray(iVertex), fMin, fMax
+
+ write(strout,fmt='(a,i5,a)') "set object ",iObject," polygon from "
+
+ ! loop over triangle vertices
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ if (iCell > mesh % nCells) then
+ atBoundary = .true.
+ end if
+
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+
+ if (.not. atBoundary) then
+ xmin = min(xmin,x) ; xmax = max(xmax,x)
+ ymin = min(ymin,y) ; ymax = max(ymax,y)
+ endif
+
+ if (iVertexDegree == 1) then
+ x0 = x
+ y0 = y
+ endif
+
+ write(stroutvertex,fmt='(f14.2,a,f14.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ enddo ! iVertexDegree
+
+ write(stroutvertex,fmt='(f14.2,a,f14.2)') x0, ",", y0
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ if (.not. atBoundary) write(55,*) trim(strout)
+
+ color = matlab_jet(fValue)
+
+ if (lnofill) then
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iObject,' fillstyle empty border lt -1'
+ else
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iObject,' fc rgb "', color, '" fillstyle solid'
+ endif
+
+ if (.not. atBoundary) iObject = iObject + 1
+
+ if (.not. atBoundary) write(55,*) trim(strout)
+
+ enddo ! iVertex
+
+ write(stroutint,fmt='(e20.8)') fMax
+ write(stroutlabel,fmt='(a,i5,a,a,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at graph 1.02, 0.9'
+ iLabel = iLabel + 1
+ write(55,*) trim(stroutlabel)
+
+
+ write(stroutint,fmt='(e20.8)') fMin
+ write(stroutlabel,fmt='(a,i5,a,a,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at graph 1.02, 0.8'
+ iLabel = iLabel + 1
+ write(55,*) trim(stroutlabel)
+
+ close(55)
+
+ end subroutine gnuplot_triangle
+
+ !-------------------------------------------------------------
+
+ subroutine gnuplot_vertexvector(mesh, vertexArrayU, vertexArrayV, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ vertexArrayU, &
+ vertexArrayV
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer, parameter :: plot_every = 10
+
+ real(kind=RKIND), parameter :: &
+ plotLength = 100000.0_RKIND
+
+ character(len=200) :: strout
+
+ real(kind=RKIND) :: &
+ x1, x2, y1, y2, &
+ n1, n2, &
+ arrowLength, &
+ arrowLengthMax, &
+ arrowLengthPlot, &
+ random
+
+ integer :: iVertex
+
+ ! find maximum vector
+ arrowLengthMax = -1e30
+
+ do iVertex = 1, mesh % nVertices
+
+ n1 = vertexArrayU(iVertex)
+ n2 = vertexArrayV(iVertex)
+ arrowLength = sqrt(n1**2 + n2**2)
+
+ arrowLengthMax = Max(arrowLengthMax,arrowLength)
+
+ enddo ! iVertex
+
+ write(*,*) arrowLengthMax
+
+
+ open(55,file=filename)
+
+ do iVertex = 1, mesh % nVertices
+
+ random = rand()
+
+ if (random <= 1.0_RKIND / real(plot_every,RKIND)) then
+
+
+ n1 = vertexArrayU(iVertex)
+ n2 = vertexArrayV(iVertex)
+ arrowLength = sqrt(n1**2 + n2**2)
+ n1 = n1 / arrowLength
+ n2 = n2 / arrowLength
+
+ arrowLengthPlot = (arrowLength / arrowLengthMax) * plotLength
+
+ x1 = mesh % xVertex % array(iVertex)
+ y1 = mesh % yVertex % array(iVertex)
+
+ x2 = x1 + n1 * arrowLengthPlot
+ y2 = y1 + n2 * arrowLengthPlot
+
+ write(*,*) arrowLengthPlot, n1, n2, x1, y1, x2, y2
+
+ write(strout,fmt='(a,i5,a,f14.2,a,f14.2,a,f14.2,a,f14.2)') "set arrow ",iArrow," from ",x1,",",y1," to ",x2,",",y2
+ iArrow = iArrow + 1
+
+ write(55,*) trim(strout)
+
+ endif
+
+ enddo ! iVertex
+
+ close(55)
+
+ end subroutine gnuplot_vertexvector
+
+ !-------------------------------------------------------------
+
+ function hexstring(red,green,blue) result(hex)
+
+ integer, intent(in) :: red,green,blue
+ character(len=7) :: hex
+
+ integer :: hex1, hex2
+
+ hex = "#"//dectohex(red)//dectohex(green)//dectohex(blue)
+
+ end function hexstring
+
+ !-------------------------------------------------------------
+
+ function dectohex(dec) result(hex)
+
+ integer, intent(in) :: dec
+ character(len=2) :: hex
+
+ integer :: hex1, hex2
+
+ character(len=1), dimension(16), parameter :: hexchars = &
+ (/'0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F'/)
+
+ hex1 = floor(real(dec, RKIND) / 16.0_RKIND)
+ hex2 = dec - hex1 * 16
+
+ hex = hexchars(hex1+1)//hexchars(hex2+1)
+
+ end function dectohex
+
+ !-------------------------------------------------------------
+
+ function matlab_jet(f) result(color)
+
+ real(kind=RKIND), intent(in) :: f
+ character(len=7) :: color
+
+ real(kind=RKIND) :: r_red, r_green, r_blue
+ integer :: red, green, blue
+
+ r_blue = max(min(min(4.0_RKIND * f + 0.5_RKIND, -4.0_RKIND * f + 2.5_RKIND),1.0_RKIND),0.0_RKIND)
+
+ r_green = max(min(min(4.0_RKIND * f - 0.5_RKIND, -4.0_RKIND * f + 3.5_RKIND),1.0_RKIND),0.0_RKIND)
+
+ r_red = max(min(min(4.0_RKIND * f - 1.5_RKIND, -4.0_RKIND * f + 4.5_RKIND),1.0_RKIND),0.0_RKIND)
+
+ red = nint(r_red * 255)
+ green = nint(r_green * 255)
+ blue = nint(r_blue * 255)
+
+ color = hexstring(red,green,blue)
+
+ end function matlab_jet
+
+ !-------------------------------------------------------------
+
+ subroutine boundary_locations(mesh, boundaryCell, boundaryVertex, interiorVertex, boundaryEdge)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:), intent(in) :: &
+ boundaryCell, &
+ boundaryVertex, &
+ interiorVertex, &
+ boundaryEdge
+
+ integer :: &
+ iCell, &
+ iVertex, &
+ iEdge
+
+ ! identify border cells
+ open(55,file="border_cell.txt",action="write")
+
+ do iCell = 1, mesh % nCells
+
+ if (boundaryCell(iCell) == 1) then
+
+ write(55,*) mesh % xCell % array(iCell), mesh % yCell % array(iCell)
+
+ endif
+
+ enddo ! iCell
+
+ close(55)
+
+ ! identify border vertices
+ open(55,file="border_vertex.txt",action="write")
+
+ do iVertex = 1, mesh % nVertices
+
+ if (boundaryVertex(iVertex) == 1) then
+
+ write(55,*) mesh % xVertex % array(iVertex), mesh % yVertex % array(iVertex)
+
+ endif
+
+ enddo ! iVertex
+
+ close(55)
+
+ ! identify interior vertices
+ open(55,file="interior_vertex.txt",action="write")
+
+ do iVertex = 1, mesh % nVertices
+
+ if (interiorVertex(iVertex) == 1) then
+
+ write(55,*) mesh % xVertex % array(iVertex), mesh % yVertex % array(iVertex)
+
+ endif
+
+ enddo ! iVertex
+
+ close(55)
+
+ ! identify border edges
+ open(55,file="border_edge.txt",action="write")
+
+ do iEdge = 1, mesh % nEdges
+
+ if (boundaryEdge(iEdge) == 1) then
+
+ write(55,*) mesh % xEdge % array(iEdge), mesh % yEdge % array(iEdge)
+
+ endif
+
+ enddo ! iVertex
+
+ close(55)
+
+ end subroutine boundary_locations
+
+ !-------------------------------------------------------------
+
+ subroutine plot_boundary_triangles(mesh, boundaryVertex, interiorVertex, boundaryEdge)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:), intent(in) :: &
+ boundaryVertex, &
+ interiorVertex, &
+ boundaryEdge
+
+ real(kind=RKIND) :: &
+ x, y, &
+ x0, y0
+
+ character(len=200) :: &
+ stroutvertex, &
+ strout
+
+ integer :: &
+ iEdgeBoundary, &
+ iEdge, &
+ iCellOnEdge, &
+ iCell, &
+ iVertexOnEdge, &
+ iVertex, &
+ iVertexDegree
+
+ integer :: start_index
+
+ start_index = mesh % nCells + mesh % nVertices + 10
+
+ open(55,file="boundary_triangles_edge.txt",action="write")
+
+ ! edge triangles
+ iEdgeBoundary = 0
+
+ do iEdge = 1, mesh % nEdges
+
+ write(0,*) iEdge, boundaryEdge(iEdge)
+
+ if (boundaryEdge(iEdge) == 1) then
+ ! am on boundary
+ iEdgeBoundary = iEdgeBoundary + 1
+
+ write(strout,fmt='(a,i5,a)') "set object ",iEdgeBoundary+start_index," polygon from "
+
+ ! find the interior cell
+ do iCellOnEdge = 1, 2
+
+ iCell = mesh % cellsOnEdge %array(iCellOnEdge,iEdge)
+
+ if (iCell >= 1 .and. iCell <= mesh % nCells) then
+
+ ! this cell is interior to domain
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+ x0 = x
+ y0 = y
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ endif
+
+ enddo ! iCellOnEdge
+
+ ! find the adjacent vertices
+ do iVertexOnEdge = 1, 2
+
+ iVertex = mesh % verticesOnEdge % array(iVertexOnEdge,iEdge)
+
+ x = mesh % xVertex % array(iVertex)
+ y = mesh % yVertex % array(iVertex)
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ enddo ! iVertexOnEdge
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2)') x0, ",", y0
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ write(55,*) trim(strout)
+
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iEdgeBoundary+start_index,' fillstyle empty border lt -1'
+
+ write(55,*) trim(strout)
+
+ endif
+
+ enddo ! iEdge
+
+ close(55)
+
+ open(55,file="boundary_triangles_vertex.txt",action="write")
+
+ ! loop over vertices
+ do iVertex = 1, mesh % nVertices
+
+ ! boundary vertex adjacent to two interior cells
+ if (boundaryVertex(iVertex) == 1) then
+
+ iEdgeBoundary = iEdgeBoundary + 1
+
+ write(strout,fmt='(a,i5,a)') "set object ",iEdgeBoundary+start_index," polygon from "
+
+ ! starting vertex point
+ x = mesh % xVertex % array(iVertex)
+ y = mesh % yVertex % array(iVertex)
+ x0 = x
+ y0 = y
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ ! the two interior cells
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ ! adjacent cell
+ iCell = mesh % cellsOnVertex % array(iVertexDegree, iVertex)
+
+ if (iCell >= 1 .and. iCell <= mesh % nCells) then
+ ! interior adjacent cell
+
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ endif
+
+ enddo ! iVertexDegree
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2)') x0, ",", y0
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ write(55,*) trim(strout)
+
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iEdgeBoundary+start_index,' fillstyle empty border lt -1'
+
+ write(55,*) trim(strout)
+
+ endif
+
+ enddo ! iVertex
+
+ close(55)
+
+ end subroutine plot_boundary_triangles
+
+ !-------------------------------------------------------------
+
+ subroutine labels_verticesOnCell(mesh, verticesOnCell, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:,:), intent(in) :: &
+ verticesOnCell
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: &
+ iCell
+
+ do iCell = 1, mesh % nCells
+
+ call label_verticesOnCell(mesh, iCell, verticesOnCell(:,iCell), filename, .true.)
+
+ enddo ! iCell
+
+ end subroutine labels_verticesOnCell
+
+ !-------------------------------------------------------------
+
+ subroutine label_verticesOnCell(mesh, iCell, verticesOnCell, filename, append)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, intent(in) :: &
+ iCell
+
+ integer, dimension(:), intent(in) :: &
+ verticesOnCell
+
+ character(len=*), intent(in) :: &
+ filename
+
+ logical, intent(in), optional :: &
+ append
+
+ real(kind=RKIND) :: &
+ x, y, x0, y0, xc, yc, xl, yl
+
+ real(kind=RKIND) :: &
+ fMin, fMax
+
+ character(len=200) :: &
+ stroutvertex, &
+ strout, &
+ stroutlabel, &
+ stroutint
+
+ integer :: &
+ iVertexOnCell, &
+ iVertex
+
+ logical :: &
+ lappend
+
+ lappend = .false.
+ if (present(append)) lappend = .true.
+
+ if (lappend) then
+ open(55,file=filename,position="append")
+ else
+ open(55,file=filename)
+ endif
+
+ write(strout,fmt='(a,i5,a)') "set object ",iObject," polygon from "
+
+ xc = mesh % xCell % array(iCell)
+ yc = mesh % yCell % array(iCell)
+
+ ! draw cell
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ x = mesh % xVertex % array(iVertex)
+ y = mesh % yVertex % array(iVertex)
+
+ if (iVertexOnCell == 1) then
+ x0 = x
+ y0 = y
+ endif
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ xl = 0.3_RKIND * xc + 0.7_RKIND * x
+ yl = 0.3_RKIND * yc + 0.7_RKIND * y
+
+ !write(stroutint,fmt='(i5)') iVertexOnCell
+ write(stroutint,fmt='(i5)') verticesOnCell(iVertexOnCell)
+ write(stroutlabel,fmt='(a,i5,a,a,a,f10.2,a,f10.2,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at ',xl, ",",yl, ' center font "Helvetica,8"'
+
+ write(55,*) trim(stroutlabel)
+
+ iLabel = iLabel + 1
+
+ enddo ! iVertexOnCell
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2)') x0, ",", y0
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ write(55,*) trim(strout)
+
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iObject,' fillstyle empty border lt -1'
+
+ write(55,*) trim(strout)
+
+ iObject = iObject + 1
+
+ close(55)
+
+ end subroutine label_verticesOnCell
+
+ !-------------------------------------------------------------
+
+ subroutine labels_verticesDegreeOnVertex(mesh, verticesDegreeOnVertex, interiorVertex, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, dimension(:,:), intent(in) :: &
+ verticesDegreeOnVertex
+
+ integer, dimension(:), intent(in) :: &
+ interiorVertex
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: &
+ iVertex
+
+ do iVertex = 1, mesh % nVertices
+
+ if (interiorVertex(iVertex) == 1) then
+
+ call label_verticesDegreeOnVertex(mesh, iVertex, verticesDegreeOnVertex(:,iVertex), filename, .true.)
+
+ endif
+
+ enddo ! iVertex
+
+ end subroutine labels_verticesDegreeOnVertex
+
+ !-------------------------------------------------------------
+
+ subroutine label_verticesDegreeOnVertex(mesh, iVertex, verticesDegreeOnVertex, filename, append)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer, intent(in) :: &
+ iVertex
+
+ integer, dimension(:), intent(in) :: &
+ verticesDegreeOnVertex
+
+ character(len=*), intent(in) :: &
+ filename
+
+ logical, intent(in), optional :: &
+ append
+
+ real(kind=RKIND) :: &
+ x, y, x0, y0, xv, yv, xl, yl
+
+ real(kind=RKIND) :: &
+ fMin, fMax
+
+ character(len=200) :: &
+ stroutvertex, &
+ strout, &
+ stroutlabel, &
+ stroutint
+
+ integer :: &
+ iVertexDegree, &
+ iCell
+
+ logical :: &
+ lappend
+
+ lappend = .false.
+ if (present(append)) lappend = .true.
+
+ if (lappend) then
+ open(55,file=filename,position="append")
+ else
+ open(55,file=filename)
+ endif
+
+ write(strout,fmt='(a,i5,a)') "set object ",iObject," polygon from "
+
+ xv = mesh % xVertex % array(iVertex)
+ yv = mesh % yVertex % array(iVertex)
+
+ ! draw cell
+ do iVertexDegree = 1, mesh % vertexDegree
+
+ iCell = mesh % cellsOnVertex % array(iVertexDegree,iVertex)
+
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+
+ if (iVertexDegree == 1) then
+ x0 = x
+ y0 = y
+ endif
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2,a)') x, ",", y, " to "
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ xl = 0.3_RKIND * xv + 0.7_RKIND * x
+ yl = 0.3_RKIND * yv + 0.7_RKIND * y
+
+ !write(stroutint,fmt='(i5)') iVertexDegree
+ write(stroutint,fmt='(i5)') verticesDegreeOnVertex(iVertexDegree)
+ write(stroutlabel,fmt='(a,i5,a,a,a,f10.2,a,f10.2,a)') "set label ",iLabel,' "',trim(adjustl(stroutint)),'" at ',xl, ",",yl, ' center font "Helvetica,8"'
+ iLabel = iLabel + 1
+
+ write(55,*) trim(stroutlabel)
+
+ enddo ! iVertexOnCell
+
+ write(stroutvertex,fmt='(f10.2,a,f10.2)') x0, ",", y0
+
+ strout = trim(strout)//trim(stroutvertex)
+
+ write(55,*) trim(strout)
+
+ write(strout,fmt='(a,i5,a,a,a)') "set object ",iObject,' fillstyle empty border lt -1'
+
+ write(55,*) trim(strout)
+
+ iObject = iObject + 1
+
+ close(55)
+
+ end subroutine label_verticesDegreeOnVertex
+
+ !-------------------------------------------------------------
+
+ subroutine writeout_minmax()
+
+ character(len=200) :: strout
+ real(kind=RKIND), parameter :: border = 7000.0_RKIND
+
+ real(kind=RKIND) :: xming, xmaxg, yming, ymaxg
+
+ xming = min(xmin,ymin)-border ; yming = xming
+ xmaxg = max(xmax,ymax)+border ; ymaxg = xmaxg
+
+ write(strout,fmt='(a,f10.2,a,f10.2,a,f10.2,a,f10.2,a)') &
+ "set size square ; set xrange [", xming, ":", xmaxg, "] ; set yrange [", yming, ":", ymaxg, "]"
+ write(*,*) trim(strout)
+
+ end subroutine writeout_minmax
+
+ !-------------------------------------------------------------
+ ! write out variables routines
+ !-------------------------------------------------------------
+
+ subroutine writeout_on_cell_real(mesh, arrayOnCell, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ arrayOnCell
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: iCell
+
+ open(55,file=filename)
+
+ do iCell = 1, mesh % nCells
+
+ write(55,*) iCell, arrayOnCell(iCell), &
+ mesh % xCell % array(iCell), mesh % yCell % array(iCell)
+
+ enddo ! iCell
+
+ close(55)
+
+ end subroutine writeout_on_cell_real
+
+ !-------------------------------------------------------------
+
+ subroutine writeout_on_vertex_real(mesh, arrayOnVertex, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:), intent(in) :: &
+ arrayOnVertex
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: iVertex
+
+ open(55,file=filename)
+
+ do iVertex = 1, mesh % nVertices
+
+ write(55,*) iVertex, arrayOnVertex(iVertex), &
+ mesh % xVertex % array(iVertex), mesh % yVertex % array(iVertex)
+
+ enddo ! iVertex
+
+ close(55)
+
+ end subroutine writeout_on_vertex_real
+
+ !-------------------------------------------------------------
+
+ subroutine writeout_on_vertex_on_cell_real(mesh, arrayVertexOnCell, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ arrayVertexOnCell
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: &
+ iCell, &
+ iVertexOnCell, &
+ iVertex
+
+ open(55,file=filename)
+
+ do iCell = 1, mesh % nCells
+
+ do iVertexOnCell = 1, mesh % nEdgesOnCell % array(iCell)
+
+ iVertex = mesh % verticesOnCell % array(iVertexOnCell,iCell)
+
+ write(55,*) iCell, iVertexOnCell, arrayVertexOnCell(iVertexOnCell,iCell), &
+ mesh % xVertex % array(iVertex), mesh % yVertex % array(iVertex), &
+ arrayVertexOnCell(iVertexOnCell,iCell) / (mesh % xVertex % array(iVertex) + mesh % yVertex % array(iVertex))
+
+ enddo ! iVertexOnCell
+
+ enddo ! iCell
+
+ close(55)
+
+ end subroutine writeout_on_vertex_on_cell_real
+
+ !-------------------------------------------------------------
+
+ subroutine writeout_array(n, arr, filename)
+
+ integer, intent(in) :: n
+
+ real(kind=RKIND), dimension(n), intent(in) :: &
+ arr
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: i
+
+ open(55,file=filename)
+
+ do i = 1, n
+
+ write(55,*) i, arr(i)
+
+ enddo ! i
+
+ close(55)
+
+ end subroutine writeout_array
+
+ !-------------------------------------------------------------
+
+ subroutine writeout_edgecell(mesh, arrayEdgeCell, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:), intent(in) :: &
+ arrayEdgeCell
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: &
+ iCell, &
+ iEdge
+
+ open(55,file=filename)
+
+ do iCell = 1, mesh % nCells
+
+ do iEdge = 1, mesh % nEdgesOnCell % array(iCell)
+
+ write(55,*) iCell, iEdge, arrayEdgeCell(iEdge,iCell)
+
+ enddo ! iEdge
+
+ enddo ! iCell
+
+ close(55)
+
+ end subroutine writeout_edgecell
+
+ !-------------------------------------------------------------
+
+ subroutine writeout_edgeedgecell(mesh, arrayEdgeEdgeCell, filename)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), dimension(:,:,:), intent(in) :: &
+ arrayEdgeEdgeCell
+
+ character(len=*), intent(in) :: &
+ filename
+
+ integer :: &
+ iCell, &
+ iEdge, &
+ jEdge
+
+ open(55,file=filename)
+
+ do iCell = 1, mesh % nCells
+
+ do iEdge = 1, mesh % nEdgesOnCell % array(iCell)
+
+ do jEdge = 1, mesh % nEdgesOnCell % array(iCell)
+
+ write(55,*) iCell, iEdge, jEdge, arrayEdgeEdgeCell(jEdge,iEdge,iCell)
+
+ enddo ! jEdge
+
+ enddo ! iEdge
+
+ enddo ! iCell
+
+ close(55)
+
+ end subroutine writeout_edgeedgecell
+
+ !-------------------------------------------------------------
+
+ function find_nearest_cell(mesh,x,y) result(iNear)
+
+ type(mesh_type), intent(in) :: mesh
+
+ real(kind=RKIND), intent(in) :: &
+ x, y
+
+ integer :: iNear
+
+ real(kind=RKIND) :: &
+ distance, &
+ min_distance
+
+ integer :: &
+ iCell
+
+ iNear = -1
+ min_distance = 1e30_RKIND
+
+ do iCell = 1, mesh % nCells
+
+ distance = sqrt((x - mesh % xCell % array(iCell))**2 + (y - mesh % yCell % array(iCell))**2)
+
+ if (distance < min_distance) then
+
+ min_distance = distance
+
+ iNear = iCell
+
+ endif
+
+ enddo ! iCell
+
+ end function find_nearest_cell
+
+ !-------------------------------------------------------------
+
+ subroutine rotate_ninety(mesh)
+
+ type(mesh_type), intent(in) :: mesh
+
+ integer :: &
+ iCell, &
+ iVertex, &
+ iEdge
+
+ real(kind=RKIND) :: x, y
+
+ do iCell = 1, mesh % nCells
+
+ x = mesh % xCell % array(iCell)
+ y = mesh % yCell % array(iCell)
+
+ mesh % xCell % array(iCell) = y
+ mesh % yCell % array(iCell) = -x
+
+ enddo ! iCell
+
+ do iVertex = 1, mesh % nVertices
+
+ x = mesh % xVertex % array(iVertex)
+ y = mesh % yVertex % array(iVertex)
+
+ mesh % xVertex % array(iVertex) = y
+ mesh % yVertex % array(iVertex) = -x
+
+ enddo ! iVertex
+
+ do iEdge = 1, mesh % nEdges
+
+ x = mesh % xEdge % array(iEdge)
+ y = mesh % yEdge % array(iEdge)
+
+ mesh % xEdge % array(iEdge) = y
+ mesh % yEdge % array(iEdge) = -x
+
+ enddo ! iVertex
+
+ end subroutine rotate_ninety
+
+ !-------------------------------------------------------------
+
+ subroutine reverse_grids(mesh)
+
+ type(mesh_type), intent(inout) :: mesh
+
+ real(kind=RKIND), dimension(:), allocatable :: tmp_ncells, tmp_nvertices
+ integer, dimension(:,:), allocatable :: tmp_maxedgesncells, tmp_vertexdegreenvertices
+
+
+ integer :: nCells0, nVertices0, vertexDegree0, maxEdges0
+
+ !write(*,*) lbound(mesh % xCell % array), ubound(mesh % xCell % array), size(mesh % xCell % array,1), mesh % nCells
+ !stop
+
+
+ nCells0 = mesh % nCells
+ nVertices0 = mesh % nVertices
+ vertexDegree0 = mesh % vertexDegree
+ maxEdges0 = mesh % maxEdges
+
+ allocate(tmp_ncells(1:mesh % nCells+1))
+ allocate(tmp_nvertices(1:mesh % nVertices+1))
+ allocate(tmp_maxedgesncells(1:mesh % maxEdges,1:mesh % nCells+1))
+ allocate(tmp_vertexdegreenvertices(1:mesh % vertexDegree,1:mesh % nVertices+1))
+
+ ! x position
+ tmp_ncells = mesh % xCell % array
+ tmp_nvertices = mesh % xVertex % array
+ deallocate(mesh % xCell % array)
+ deallocate(mesh % xVertex % array)
+ allocate(mesh % xCell % array(1:nVertices0+1))
+ allocate(mesh % xVertex % array(1:nCells0+1))
+ mesh % xCell % array = tmp_nvertices
+ mesh % xVertex % array = tmp_ncells
+
+ ! x position
+ tmp_ncells = mesh % yCell % array
+ tmp_nvertices = mesh % yVertex % array
+ deallocate(mesh % yCell % array)
+ deallocate(mesh % yVertex % array)
+ allocate(mesh % yCell % array(1:nVertices0+1))
+ allocate(mesh % yVertex % array(1:nCells0+1))
+ mesh % yCell % array = tmp_nvertices
+ mesh % yVertex % array = tmp_ncells
+
+ ! z position
+ tmp_ncells = mesh % zCell % array
+ tmp_nvertices = mesh % zVertex % array
+ deallocate(mesh % zCell % array)
+ deallocate(mesh % zVertex % array)
+ allocate(mesh % zCell % array(1:nVertices0+1))
+ allocate(mesh % zVertex % array(1:nCells0+1))
+ mesh % zCell % array = tmp_nvertices
+ mesh % zVertex % array = tmp_ncells
+
+ ! nEdgesOnCell
+ deallocate(mesh % nEdgesOnCell % array)
+ allocate(mesh % nEdgesOnCell % array(nVertices0+1))
+ mesh % nEdgesOnCell % array = 3
+
+
+ ! verticesOnCell / cellsOnVertex
+ tmp_maxedgesncells = mesh % verticesOnCell % array
+ tmp_vertexdegreenvertices = mesh % cellsOnVertex % array
+ deallocate(mesh % verticesOnCell % array)
+ deallocate(mesh % cellsOnVertex % array)
+ allocate(mesh % verticesOnCell % array(vertexDegree0,nVertices0+1))
+ allocate(mesh % cellsOnVertex % array(maxEdges0,nCells0+1))
+ mesh % verticesOnCell % array = tmp_vertexdegreenvertices
+ mesh % cellsOnVertex % array = tmp_maxedgesncells
+
+ ! area
+ deallocate(mesh % areaCell % array)
+ allocate(mesh % areaCell % array(nVertices0+1))
+ mesh % areaCell % array = mesh % areaTriangle % array
+
+ ! sizes
+ mesh % nCells = nVertices0
+ mesh % nVertices = nCells0
+ mesh % vertexDegree = maxEdges0
+ mesh % maxEdges = vertexDegree0
+
+ deallocate(tmp_ncells)
+ deallocate(tmp_nvertices)
+ deallocate(tmp_maxedgesncells)
+ deallocate(tmp_vertexdegreenvertices)
+
+ !stop
+ end subroutine reverse_grids
+
+ !-------------------------------------------------------------
+
+end module cice_testing
+
Added: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_time_integration.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_time_integration.F (rev 0)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_time_integration.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -0,0 +1,24 @@
+module cice_time_integration
+
+ use mpas_grid_types
+
+contains
+
+ !--------------------------------------------------------------------------
+
+ subroutine cice_timestep(block, dt)
+
+ use cice_dynamics_shared, only: run_dynamics
+
+ implicit none
+
+ type(block_type), intent(inout) :: block
+ real(kind=RKIND), intent(in) :: dt
+
+ call run_dynamics(block % mesh, block % icestate, block % weak, block % boundary, dt)
+
+ end subroutine cice_timestep
+
+ !--------------------------------------------------------------------------
+
+end module cice_time_integration
Deleted: branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_timestep.F
===================================================================
--- branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_timestep.F 2013-04-19 17:49:14 UTC (rev 2779)
+++ branches/cice_projects/initial_cice_core/src/core_cice/mpas_cice_timestep.F 2013-04-19 18:44:55 UTC (rev 2780)
@@ -1,139 +0,0 @@
-module cice_timestep
-
- use mpas_grid_types
-
- implicit none
-
-
-contains
-
- !-------------------------------------------------------------------------
-
- subroutine mpas_timestep(domain, itimestep, dt, timeStamp)
-
- type (domain_type), intent(inout) :: domain
- integer, intent(in) :: itimestep
- real (kind=RKIND), intent(in) :: dt
- character(len=*), intent(in) :: timeStamp
-
- call step_therm1(domain, dt)
-
-
-
- end subroutine mpas_timestep
-
- !-------------------------------------------------------------------------
-
- subroutine step_therm1(domain, dt)
-
- use cice_test
- use cice_state, only: ntrcr, nilyr, nt_qice, nt_qsno
-
- type (domain_type), intent(inout) :: domain
- real (kind=RKIND), intent(in) :: dt
- integer :: nCells, nCategories
- integer :: iCell
-
- type (block_type), pointer :: block
- type (states_type), pointer :: states
- type (tracer_type), pointer :: tracer
- type (flux_type), pointer :: flux
-
- ! original non-local variables
- real(kind=RKIND), dimension(1,1,1) :: &
- aice , &
- sst , &
- Tf , &
- strocnxT , &
- strocnyT , &
- rside , &
- frzmlt
-
- real(kind=RKIND), dimension(:,:,:,:), allocatable :: &
- vicen , &
- vsnon
-
- real(kind=RKIND), dimension(:,:,:,:,:), allocatable :: &
- trcrn
-
- ! original step_therm1 variables
- real(kind=RKIND), dimension(1,1) :: &
- Tbot, &
- fbot
-
- ! allocate local cice variables
- allocate(vicen(1,1,nCategories,1))
- allocate(vsnon(1,1,nCategories,1))
- allocate(trcrn(1,1,ntrcr,nCategories,1))
-
- ! loop over blocks
- block => domain % blocklist
- do while (associated(block))
-
- ! loop over all cells
- do iCell = 1, nCells
-
- ! structure pointers
- states => block % states
- tracer => block % tracer
- flux => block % flux
-
- ! state variable conversion
- aice(1,1,1) = states % iceAreaCell % array(iCell)
- vicen(1,1,:,1) = states % iceVolumeCategory % array(iCell,:)
- vsnon(1,1,:,1) = states % snowVolumeCategory % array(iCell,:)
-
- ! flux variable conversion
- sst(1,1,1) = flux % seaSurfaceTemperature % array(iCell)
- Tf(1,1,1) = flux % seaFreezingTemperature % array(iCell)
- rside(1,1,1) = flux % iceLateralMeltFraction % array(iCell)
- frzmlt(1,1,1) = flux % seaFreezeMeltPotential % array(iCell)
-
- strocnxT(1,1,1) = 0.0_RKIND
- strocnyT(1,1,1) = 0.0_RKIND
-
- ! tracer conversion
- trcrn(1,1,nt_qice:nt_qice+nilyr-1,:,1) = &
- tracer % iceVolumeLayerTracers % array(tracer % index_iceEnthalpy,:,iCell,:)
- trcrn(1,1,nt_qsno:nt_qsno+nilyr-1,:,1) = &
- tracer % snowVolumeLayerTracers % array(tracer % index_snowEnthalpy,:,iCell,:)
-
-!!! BEGIN ORIGINAL CICE CALL
- call frzmlt_bottom_lateral(1, 1, &
- 1, 1, 1, 1, &
- ntrcr, dt, &
- aice (:,:,:), frzmlt (:,:,:), &
- vicen (:,:,:,1), vsnon (:,:,:,:), &
- trcrn (:,:,1:ntrcr,:,:), &
- sst (:,:,:), Tf (:,:,:), &
- strocnxT(:,:,:), strocnyT(:,:,:), &
- Tbot, fbot, &
- rside (:,:,:))
-!!! END ORIGINAL CICE CALL
-
- ! state variable back conversion
- states % iceAreaCell % array(iCell) = aice(1,1,1)
- states % iceVolumeCategory % array(iCell,:) = vicen(1,1,:,1)
- states % snowVolumeCategory % array(iCell,:) = vsnon(1,1,:,1)
-
- ! tracer back conversion
- tracer % iceVolumeLayerTracers % array(tracer % index_iceEnthalpy,:,iCell,:) = &
- trcrn(1,1,nt_qice:nt_qice+nilyr-1,:,1)
- tracer % snowVolumeLayerTracers % array(tracer % index_snowEnthalpy,:,iCell,:) = &
- trcrn(1,1,nt_qsno:nt_qsno+nilyr-1,:,1)
-
- enddo ! iCell
-
- block => block % next
- end do ! blocks
-
- ! clean up original variables
- deallocate(vicen)
- deallocate(vsnon)
- deallocate(trcrn)
-
- end subroutine step_therm1
-
- !-------------------------------------------------------------------------
-
-end module cice_timestep
</font>
</pre>