1!*****************************************************************************/
2! *
3! *  Elmer, A Finite Element Software for Multiphysical Problems
4! *
5! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6! *
7! *  This program is free software; you can redistribute it and/or
8! *  modify it under the terms of the GNU General Public License
9! *  as published by the Free Software Foundation; either version 2
10! *  of the License, or (at your option) any later version.
11! *
12! *  This program is distributed in the hope that it will be useful,
13! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15! *  GNU General Public License for more details.
16! *
17! *  You should have received a copy of the GNU General Public License
18! *  along with this program (in file fem/GPL-2); if not, write to the
19! *  Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
20! *  Boston, MA 02110-1301, USA.
21! *
22! *****************************************************************************/
23!
24!/******************************************************************************
25! *
26! *  Solver to predict calving events in 2D based on crevasse depth calving
27! *  criterion (Benn et al. 2007). Description of algorithm in
28! *  Todd & Christoffersen (2014) [doi:10.5194/tc-8-2353-2014]
29! *
30! *  Relies on TwoMeshes.F90 for mesh migration and interpolation
31! *  following calving, and FrontDisplacement.F90 for mesh update computation
32! *
33! ******************************************************************************
34! *
35! *  Authors: Joe Todd
36! *  Email:   jat39@st-andrews.ac.uk
37! *  Web:     http://www.csc.fi/elmer
38! *  Address: CSC - IT Center for Science Ltd.
39! *           Keilaranta 14
40! *           02101 Espoo, Finland
41! *
42! *  Original Date: 10.11.2014
43! *
44! ****************************************************************************/
45SUBROUTINE Find_Calving (Model, Solver, dt, TransientSimulation )
46
47   USE types
48   USE CoordinateSystems
49   USE SolverUtils
50   USE ElementDescription
51   USE DefUtils
52
53   IMPLICIT NONE
54
55   TYPE CrevasseGroups_t
56      LOGICAL, ALLOCATABLE :: NotEmpty(:),Valid(:)
57      INTEGER, ALLOCATABLE :: NodeIndexes(:,:), NoNodes(:)
58      INTEGER :: CurrentGroup
59      !This derived type is for storing groups of connected nodes which
60      !have surface or basal crevasses.  I can't think of a sensible strategy
61      !to determine the number of groups required, nor the number of nodes in
62      !each group.  However, given that they only contain integers and logicals,
63      !declaring them massive probably isn't an issue.
64   END TYPE CrevasseGroups_t
65
66   TYPE(Model_t) :: Model
67   TYPE(Solver_t) :: Solver
68   Type(Nodes_t) :: CornerElementNodes, CurrentElementNodes, &
69        TargetNodes, CalvedNodes
70   Type(Nodes_t), POINTER :: Nodes0
71   Type(Nodes_t), TARGET :: EvalNodes
72   TYPE(Matrix_t), POINTER :: Vel1Matrix !For finding neighbours
73   TYPE(ValueList_t), POINTER :: Material, SolverParams
74   TYPE(Element_t), POINTER :: CurrentElement, CornerElement
75   TYPE(CrevasseGroups_t) :: SurfaceCrevasseGroups, BasalCrevasseGroups
76   TYPE(Mesh_t), POINTER :: Mesh, EvalMesh, Mesh0 => Null()
77   TYPE(Variable_t), POINTER :: DepthSol, StressSol, Stress1Sol, Stress4Sol, &
78        Vel1Sol, DistanceSol, FlowSol, &
79        WPressureSol, TimeVar, CSurfIndexSol, CBasalIndexSol, Var, &
80        CrevasseGroupVar
81
82   REAL (KIND=dp) :: t, dt, xcoord, ycoord, NodeDepth, NodeStress, NodeStress1, NodeStress4, &
83        NodeWPressure, CalvingSurfIndex, CalvingBasalIndex, CalvingCoordHolder = 1.0E20, &
84        FreeConnect, FreeConnected, Length, WaterDepth, Dw, NodeLength, RhoWF, RhoWS, Rho, &
85        g, OverlapCalvingCoordinate, SeaLevel, BasalCalvingCoordinate, EvalResolution, &
86        STuningParam, BTuningParam, DwStart, DwStop, season, Te, sign, Normal(3), &
87        CornerNormal(3), BedSecond, BedSecondDiff, beddiff, BedToler, dx, dy, &
88        LocalDist, LocalDistNode, PropDistNode, normalcond, ForceCalveSize, &
89        localM(2,2), EigValues(2), EI(2), dumy, dumy2, work(8),YieldStress, &
90#ifdef USE_ISO_C_BINDINGS
91        rt0, rt
92#else
93        rt0, rt, RealTime
94#endif
95   REAL (KIND=DP), POINTER :: DepthValues(:), StressValues(:), Stress1Values(:), Stress4Values(:), &
96        Calving1Values(:), DistanceValues(:), WPressureValues(:), FrontValues(:), &
97        CSurfIndexValues(:), CBasalIndexValues(:), CIndexValues(:), CrevasseGroupValues(:), &
98        Calving2Values(:), EvalPoints(:,:), Field(:), WorkReal(:)
99   REAL (KIND=dp), ALLOCATABLE :: CumDist(:), PropCumDist(:), TargetCumDist(:),TargetPropDist(:)
100   INTEGER :: DIM, i, j, NoNodes, MaxN, FrontNodes, BotNodes, TopNodes, &
101        NofEvalPoints, NextSlot, BotCornerIndex, &
102        BotSecondIndex, county, GoToNode, PrevNode, NextNode, &
103        TimeSinceLast, NoNeighbours, MaxNeighbours, DOFs, ierr, StressDOFs
104
105   INTEGER, POINTER :: DepthPerm(:), StressPerm(:), Stress1Perm(:), Stress4Perm(:),&
106        Vel1Perm(:), Vel1InvPerm(:), DistancePerm(:), OrderPerm(:),&
107        WPressurePerm(:), FrontPerm(:)=>NULL(), InvFrontPerm(:), &
108        TopPerm(:)=>NULL(), BotPerm(:)=>NULL(), Permutation(:), &
109        CSurfIndexPerm(:), CBasalIndexPerm(:), CIndexPerm(:), FieldPerm(:), &
110        NodeNeighbours(:,:), NumNeighbours(:), CrevasseGroupPerm(:)
111
112   INTEGER, ALLOCATABLE :: ThisNodeNeighbours(:)
113   LOGICAL :: FirstTime = .TRUE., CalvingOccurs, RemeshOccurs, &
114        BasalCalvingOccurs = .FALSE., BasalCrevasseModel, OverlapOccurs, DwSeason, &
115        CornerCalving, KeepLooking, TransientSimulation, Found, Debug = .FALSE., &
116        ForceCalving = .FALSE., PlaneStressOnly, Cauchy, OldWay, BasalFS, CornerBadBed, &
117        CornerBadSlope
118   CHARACTER(LEN=MAX_NAME_LEN) :: SolverName, FrontMaskName, BotMaskName, TopMaskName, &
119        DwMode, CalvingModel, FlowVarName,BasalFSVarName
120
121   SAVE :: FirstTime, DIM, Permutation, Calving1Values, Calving2Values, &
122        NoNodes, Mesh, FrontMaskName, FreeConnect, WaterDepth, Rho, RhoWS, RhoWF, g, &
123        FrontPerm, FrontValues, FrontNodes, BotPerm, BotNodes, TopPerm, &
124        TopNodes, CSurfIndexSol, &
125        CBasalIndexSol, CSurfIndexPerm, CBasalIndexPerm, CSurfIndexValues, CBasalIndexValues, &
126        EvalMesh, Material, EvalNodes, NodeNeighbours, NumNeighbours, Mesh0, Nodes0, &
127        TimeSinceLast, BasalCrevasseModel, PlaneStressOnly, Cauchy, SolverParams, OldWay, &
128        BasalFS, BasalFSVarName
129
130   INTERFACE
131      SUBROUTINE InterpolateMeshToMesh( OldMesh, NewMesh, OldVariables, &
132           NewVariables, UseQuadrantTree, Projector, MaskName, UnfoundNodes )
133        !------------------------------------------------------------------------------
134        USE Lists
135        USE SParIterComm
136        USE Interpolation
137        USE CoordinateSystems
138        !-------------------------------------------------------------------------------
139        TYPE(Mesh_t), TARGET  :: OldMesh, NewMesh
140        TYPE(Variable_t), POINTER, OPTIONAL :: OldVariables, NewVariables
141        LOGICAL, OPTIONAL :: UseQuadrantTree
142        LOGICAL, POINTER, OPTIONAL :: UnfoundNodes(:)
143        TYPE(Projector_t), POINTER, OPTIONAL :: Projector
144        CHARACTER(LEN=*),OPTIONAL :: MaskName
145      END SUBROUTINE InterpolateMeshToMesh
146   END INTERFACE
147
148   rt0 = RealTime() !time it!
149
150   SolverName = "Calving"
151
152   Debug = .FALSE.
153
154   IF(ParEnv % Pes > 1) CALL Fatal(SolverName, "Solver does not work in parallel!")
155
156   Timevar => VariableGet( Model % Variables,'Time', UnfoundFatal=.TRUE.)
157   t = TimeVar % Values(1)
158   dt = Model % Solver % dt
159   season = t - floor(t)
160
161   RemeshOccurs = .FALSE. !For undercutting
162
163   !*****************************Get Variables*******************************
164
165   !This is the main solver variable, which only has a value on the calving
166   !face of the glacier, and corresponds to the magnitude of retreat.
167   IF(.NOT. ASSOCIATED(Solver % Variable)) CALL Fatal('Calving', 'Variable not associated!')
168
169   !This is the Calving Solver's Surface Crevasse exported variable, which has
170   !a positive value where a crevasse is predicted to exist, and a negative value elsewhere.
171   CSurfIndexSol => VariableGet( Solver % Mesh % Variables, 'Calving Surface Index', &
172        UnfoundFatal=.TRUE.)
173   CSurfIndexPerm => CSurfIndexSol % Perm
174   CSurfIndexValues => CSurfIndexSol % Values
175
176   !This is the Calving Solver's Basal crevasse exported variable.
177   !Positive value where basal crevasse
178   !is predicted to exist, negative otherwise, and zero at the tip.
179
180   !From Faezeh: C_basal = LongitudinalDevStress(Rxx)-Rho*g*(IceThickness-
181   !HeightAboveGlacierBase)+RhoOcean*g*(PiezometricHeight-HeightAboveGlacierBase)
182
183   CBasalIndexSol => VariableGet( Solver % Mesh % Variables, 'Calving Basal Index', &
184        UnfoundFatal=.TRUE.)
185   CBasalIndexPerm => CBasalIndexSol % Perm
186   CBasalIndexValues => CBasalIndexSol % Values
187
188   CrevasseGroupVar => VariableGet( Solver % Mesh % Variables, 'Crevasse Group ID', &
189        UnfoundFatal=.TRUE.)
190   CrevasseGroupPerm => CrevasseGroupVar % Perm
191   CrevasseGroupValues => CrevasseGroupVar % Values
192
193   DepthSol => VariableGet( Model % Variables, 'Depth', UnfoundFatal=.TRUE.)
194   DepthPerm => DepthSol % Perm
195   DepthValues => DepthSol % Values
196
197   StressSol => VariableGet( Model % Variables, "Stress", UnfoundFatal=.TRUE.)
198   StressPerm => StressSol % Perm
199   StressValues => StressSol % Values
200   StressDOFs = StressSol % DOFs
201
202   Stress1Sol => VariableGet( Model % Variables, "Stress 1", &
203        UnfoundFatal=.TRUE.)
204   Stress1Perm => Stress1Sol % Perm
205   Stress1Values => Stress1Sol % Values
206
207   Stress4Sol => VariableGet( Model % Variables, "Stress 4", &
208        UnfoundFatal=.TRUE.)
209   Stress4Perm => Stress4Sol % Perm
210   Stress4Values => Stress4Sol % Values
211
212   DistanceSol => VariableGet( Model % Variables, "Distance", &
213        UnfoundFatal=.TRUE.)
214   DistancePerm => DistanceSol % Perm
215   DistanceValues => DistanceSol % Values
216
217   WPressureSol => VariableGet( Model % Variables, "Water Pressure", &
218        UnfoundFatal=.TRUE.)
219   WPressurePerm => WPressureSol % Perm
220   WPressureValues => WPressureSol % Values
221
222   IF(FirstTime) THEN
223
224      DIM = CoordinateSystemDimension()
225      IF(DIM /= 2) CALL Fatal('Calving','Solver only works in 2D!')
226      MaxN = Model % Mesh % MaxElementNodes
227      !To track time since last calving event. Probably not needed.
228      TimeSinceLast = 0
229
230      Mesh => Solver % Mesh
231      NoNodes = Mesh % NumberOfNodes
232
233      FrontMaskName = 'Calving Front Mask'
234      TopMaskName = 'Top Surface Mask'
235      BotMaskName = 'Bottom Surface Mask'
236      ALLOCATE( FrontPerm(NoNodes), TopPerm(NoNodes), BotPerm(NoNodes))
237
238      CALL MakePermUsingMask( Model,Solver,Mesh,FrontMaskName, &
239           .FALSE., FrontPerm, FrontNodes )
240      CALL MakePermUsingMask( Model,Solver,Mesh,TopMaskName, &
241           .FALSE., TopPerm, TopNodes )
242      CALL MakePermUsingMask( Model,Solver,Mesh,BotMaskName, &
243           .FALSE., BotPerm, BotNodes )
244
245      !Holds the variable values
246      ALLOCATE( FrontValues(FrontNodes * 2) )
247
248      Mesh0 => AllocateMesh()
249      Mesh0 = Mesh
250      Mesh0 % Name = TRIM(Mesh % Name)//'_initial'
251      CALL Info('Calving','Created initial reference mesh to remap front, maintaining quality')
252      ALLOCATE( Nodes0 )
253      ALLOCATE( Nodes0 % x(NoNodes), Nodes0 % y(NoNodes), Nodes0 % z(NoNodes) )
254      Nodes0 % x = Mesh % Nodes % x
255      Nodes0 % y = Mesh % Nodes % y
256      Nodes0 % z = Mesh % Nodes % z
257      Mesh0 % Nodes => Nodes0
258
259      Permutation => FrontPerm
260      Calving1Values => FrontValues(1::2)
261      Calving2Values => FrontValues(2::2)
262
263      !Initialize
264      Calving1Values = 0.0_dp
265      Calving2Values = 0.0_dp
266
267      !Potential to force an initial calving event
268      !Specified by LOGICAL and REAL in SIF
269      SolverParams => GetSolverParams()
270      ForceCalving = GetLogical(SolverParams, 'Force Calving', Found)
271      !TODO: This doesn't work - Calving1Values isn't associated with Solver % Variable % Values yet...
272      IF(ForceCalving) THEN
273         ForceCalveSize = GetConstReal (SolverParams, 'Force Calving Size', Found)
274         IF(.NOT.Found) CALL Fatal(SolverName, "Force Calving Size not found.")
275         Calving1Values = -ForceCalveSize
276         CalvingOccurs =.TRUE.
277         CALL ListAddLogical( Model % Simulation, 'CalvingOccurs', CalvingOccurs )
278         RETURN
279      END IF
280   ENDIF
281
282   Solver % Variable % Values => FrontValues
283   Solver % Variable % Perm => FrontPerm
284
285   Var => VariableGet(Solver % Mesh % Variables, ComponentName(Solver % Variable % Name, 1), .TRUE.)
286   Var % Values => Calving1Values
287   Var % Perm => Permutation
288   Var => VariableGet(Solver % Mesh % Variables, ComponentName(Solver % Variable % Name, 2), .TRUE.)
289   Var % Values => Calving2Values
290   Var % Perm => Permutation
291
292   IF(FirstTime .OR. Solver % Mesh % Changed) THEN
293      FirstTime = .FALSE.
294
295      !STRATEGY: Finding neighbours on the fly works fine UNLESS you are in a recursive subroutine
296      !Then it messes up, because at each level, ThisNodeNeighbours is deallocate and reallocated,
297      !meaning that when you jump back up, info is already overwritten. SO:
298      !Keep the structure as it was with CalvingNeighbours, cycle as below to fill it, and ALSO
299      !create an array to hold the number of neighbours for each node
300
301      !Get the Matrix of the N-S Solver
302      Vel1Sol => VariableGet( Solver % Mesh % Variables, 'Velocity 1', UnfoundFatal=.TRUE.)
303      Vel1Perm => Vel1Sol % Perm
304      Vel1Matrix => Vel1Sol % Solver % Matrix
305
306      !Vel * DIM + Pressure...
307      DOFs = DIM + 1
308      MaxNeighbours = DIM * 10  !totally arbitrary...
309
310      !Create inverse perm to lookup Matrix later
311      ALLOCATE(Vel1InvPerm(MAXVAL(Vel1Perm)*DOFs)) !TODO DEALLOCATE
312      !2D array to hold each nodes neighbours
313      ALLOCATE(NodeNeighbours(NoNodes,MaxNeighbours))
314      !1D array to hold number of neighbours for each node
315      ALLOCATE(NumNeighbours(NoNodes))
316      NodeNeighbours = 0
317      NumNeighbours = 0
318      Vel1InvPerm = 0
319
320      j = 0
321      DO i=1,SIZE(Vel1Perm)
322         IF(Vel1Perm(i) == 0) CYCLE
323         j = j + 1
324         Vel1InvPerm( (Vel1Perm(i)*DOFs-2) : (Vel1Perm(i)*DOFs) ) = j !The 2 here is suspect...
325      END DO
326
327      DO i = 1,NoNodes
328         CALL FindNodeNeighbours(i) !Updates the allocatable array 'ThisNodeNeighbours'
329         NumNeighbours(i) = SIZE(ThisNodeNeighbours)
330         NodeNeighbours(i,1:NumNeighbours(i)) = ThisNodeNeighbours
331      END DO
332   END IF
333
334   PRINT *, '****Calving Timestep : ',t
335
336   rt = RealTime() - rt0
337   IF(ParEnv % MyPE == 0) &
338        PRINT *, 'Calving, time taken for variable loading etc:', rt
339   rt0 = RealTime()
340
341   MaxN = Model % Mesh % MaxElementNodes
342   OldWay = ListGetLogical(SolverParams, "Old Way", Found)
343   IF(.NOT. Found) OldWay = .FALSE.
344
345   !Identify the basal freesurface variable
346   BasalFS = ListGetLogical(SolverParams, "Basal FreeSurface", Found)
347   IF(.NOT. Found) BasalFS = .TRUE.
348   IF(BasalFS) THEN
349     BasalFSVarName = ListGetString(SolverParams, "Basal FreeSurface Variable Name",Found)
350     IF(.NOT. Found) CALL Fatal(SolverName, &
351          "Basal FreeSurface = True but no Basal FreeSurface Variable Name found!")
352   END IF
353
354   Material => GetMaterial(Model % Mesh % Elements(Solver % ActiveElements(1)))
355   Cauchy = ListGetLogical( Material , 'Cauchy', Found )
356   IF(.NOT. Found) THEN
357     Cauchy = .FALSE.
358     CALL Warn(SolverName, "Couldn't find 'cauchy' logical in material, &
359          &assuming deviatoric stress.")
360   END IF
361
362   FlowVarName = ListGetString(SolverParams, "Flow Solution Variable Name",Found)
363   IF(.NOT. Found) FlowVarName = "Flow Solution"
364
365   FlowSol => VariableGet(Model % Variables, FlowVarName, .TRUE., UnfoundFatal=.TRUE.)
366
367   !! FreeConnect: distance from calving front where free connection to proglacial
368   !!water body exists. Declare in SIF Constants!
369   FreeConnect = GetConstReal( Model % Constants, "FreeConnect", Found )
370   IF(.NOT.Found) CALL Fatal(SolverName, "Unable to find FreeConnect")
371   Dw = GetConstReal( Model % Constants, "Water Depth", Found )
372   IF(.NOT.Found) THEN
373     CALL Warn(SolverName, "Unable to find Water Depth, assuming zero.")
374     Dw = 0.0_dp
375   END IF
376
377   Rho = GetConstReal( Model % Constants, "Rho", Found )
378   IF(.NOT.Found) CALL Fatal(SolverName, "Unable to find Rho.")
379   RhoWS = GetConstReal( Model % Constants, "RhoWS", Found )
380   IF(.NOT.Found) CALL Fatal(SolverName, "Unable to find RhoW.")
381   RhoWF = GetConstReal( Model % Constants, "RhoWF", Found )
382   IF(.NOT.Found) CALL Fatal(SolverName, "Unable to find RhoWF.")
383   g = GetConstReal( Model % Constants, "g", Found )
384   IF(.NOT.Found) CALL Fatal(SolverName, "Unable to find Gravity.")
385   IF(g < 0.0_dp) g = -1.0_dp * g
386
387   SeaLevel = GetConstReal( Model % Constants, "Sea Level", Found )
388   IF(.NOT.Found) THEN
389     WRITE(Message,'(A)') 'Variable Sea level not found. &
390          &Setting to 0.0'
391     CALL Info('Calving', Message, level=2)
392     SeaLevel = 0.0_dp
393   END IF
394
395   !This next bit stolen from MaterialModels.src
396   !Only needed if calvingmodel=planestrain, could put logical?
397   STuningParam = GetConstReal( Model % Constants, "Surface Crevasse Tuning Parameter", Found )
398   IF(.NOT.Found) THEN
399     WRITE(Message,'(A)') 'Surface Crevasse Tuning Parameter not found. &
400          Setting to 1.0'
401     CALL Info('Calving', Message, level=2)
402     STuningParam = 1.0_dp
403   END IF
404
405   BTuningParam = GetConstReal( Model % Constants, "Basal Crevasse Tuning Parameter", Found )
406   IF(.NOT.Found) THEN
407     WRITE(Message,'(A)') 'Basal Crevasse Tuning Parameter not found. &
408          Setting to 1.0'
409     CALL Info('Calving', Message, level=2)
410     BTuningParam = 1.0_dp
411   END IF
412
413   !Get stuff for variable water depth in crevasses.
414   !Don't think this is the logical place for this, but whatever.
415   SolverParams => GetSolverParams()
416
417   DwMode = GetString (SolverParams, 'Dw Mode', Found)
418   IF(.NOT.Found) THEN
419     CALL Warn(SolverName, "Dw Mode not found, assuming off.")
420     DwMode = "off"
421   END IF
422   DwMode = TRIM(DwMode)
423   IF(DwMode /= "off") THEN
424     DwStart = GetConstReal (SolverParams, 'Dw Start', Found)
425     IF(.NOT.Found) CALL Fatal(SolverName, "Dw Start not found.")
426
427     DwStop = GetConstReal (SolverParams, 'Dw Stop', Found)
428     IF(.NOT.Found) CALL Fatal(SolverName, "Dw Stop not found.")
429   END IF
430
431   YieldStress = ListGetConstReal(SolverParams, "Yield Stress", Found)
432   IF(.NOT. Found) THEN
433     YieldStress = 0.0_dp
434     CALL Warn(SolverName, "No yield stress found, assuming 0.0")
435   ELSE
436     WRITE(Message,'(A,f8.2)') 'Yield Stress = ',YieldStress
437     CALL Info(SolverName, Message)
438   END IF
439
440   BasalCrevasseModel = ListGetLogical(SolverParams, "Basal Crevasse Model", Found)
441   IF(.NOT. Found) THEN
442     BasalCrevasseModel = .TRUE.
443     CALL Warn(SolverName,"No 'Basal Crevasse Model' found, assuming True")
444   END IF
445   IF(BasalCrevasseModel) THEN
446     CALL Info(SolverName,"Using Basal Crevasse Model, &
447          &calving occurs when surface and basal crevasses meet", Level=2)
448   ELSE
449     CALL Info(SolverName,"Using Surface Crevasse Model, &
450          &calving occurs when surface crevasses meet sea level",Level=2)
451   END IF
452
453   PlaneStressOnly = ListGetLogical(SolverParams, "Plane Stress Only", Found)
454   IF(.NOT. Found) PlaneStressOnly = .FALSE.
455
456   SELECT CASE(DwMode)
457   CASE("off")
458     WaterDepth = 0.0
459   CASE ("constant")
460
461     WaterDepth = Dw
462
463   CASE("binary")
464
465     DwSeason = .FALSE.
466     IF(DwStart .LT. DwStop) THEN !normal case, dw season doesn't pass winter
467       IF((season .GT. DwStart) .AND. (season .LT. DwStop)) DwSeason = .TRUE.
468     ELSE !this is unlikely to be needed
469       IF((season .GT. DwStart ) .OR. (season .LT. DwStop )) DwSeason = .TRUE.
470     END IF
471     IF(DwSeason) THEN
472       WaterDepth = Dw
473     ELSE
474       WaterDepth = 0.0
475     ENDIF
476
477   CASE DEFAULT
478     CALL Fatal(SolverName,"Invalid Dw mode selection")
479   END SELECT
480
481   !Finds the maximum value of Coordinate 1 (i.e. the calving face) and
482   !FreeConnected (the minimum Coordinate 1 which should be checked for calving.
483
484   !First find the length of the glacier
485   Length = 0.0
486   DO i = 1, NoNodes
487     IF(Permutation(i) == 0) CYCLE
488     NodeLength = Model % Nodes % x(i)
489     IF(NodeLength > Length) Length = NodeLength
490   END DO
491   PRINT *, '**** Calving Glacier Length = ',Length
492   CalvingCoordHolder = Length
493   FreeConnected = Length - FreeConnect
494
495   !---------------------------------------
496   !
497   ! Evaluate the crevasse criteria for basal and surface crevasses
498   !
499   !---------------------------------------
500
501   DO i = 1, NoNodes
502     xcoord = Mesh % Nodes % x(i)
503     ycoord = Mesh % Nodes % y(i)
504
505     !TODO - Tuning Parameters aren't used except in specific PlaneStressOnly case.
506     !SORT THIS OUT - but what sense does the tuning parameter have in the cauchy case?
507     IF(.NOT. OldWay) THEN
508
509       !Compute maximum extensive principal stress
510       localM=0.0_dp
511       !xx,yy,zz,xy,yz,zx
512       localM(1,1)=StressValues( StressDOFs * (StressPerm(i) - 1 ) + 1) !xx
513       localM(2,2)=StressValues( StressDOFs * (StressPerm(i) - 1 ) + 2) !yy
514       localM(1,2)=StressValues( StressDOFs * (StressPerm(i) - 1 ) + 4) !xy
515       localM(2,1)=localM(1,2)
516
517       IF(.NOT. Cauchy) THEN
518         DO j=1,2                                       !dim+1, only runs in 3D
519           localM(j,j)=localM(j,j) - FlowSol % Values(FlowSol % Perm(i)*FlowSol % DOFs)
520         END DO
521       END IF
522
523       CALL DGEEV('N','N',2,localM,2,EigValues,EI,Dumy,1,Dumy2,1,Work,8,ierr )
524       IF (ierr/=0) THEN
525         WRITE(Message,'(A,i0)') 'Failed to compute EigenValues, error code: ',ierr
526         CALL FATAL(SolverName, Message)
527       END IF
528
529       NodeStress = MAXVAL(EigValues) !Maximum principalstress
530
531       NodeWPressure = -WPressureValues(WPressurePerm(i))
532       IF(NodeWPressure < 0.0_dp) NodeWPressure = 0.0_dp
533
534       CalvingSurfIndex = NodeStress + (WaterDepth * RhoWF * g) - YieldStress
535
536       CalvingBasalIndex = NodeStress + NodeWPressure - YieldStress
537     ELSE
538       CALL Info('Calving', 'Calculating calving criterion the old way', level=2)
539       NodeDepth = DepthValues(DepthPerm(i))
540
541       NodeStress1 = Stress1Values(Stress1Perm(i))
542       IF(Cauchy) &
543            NodeStress1 = NodeStress1 + FlowSol % Values(FlowSol % Perm(i)*FlowSol % DOFs)
544
545
546       NodeWPressure = -WPressureValues(WPressurePerm(i))
547       IF(NodeWPressure < 0.0_dp) NodeWPressure = 0.0_dp
548
549       !NodeWPressure is the water pressure in a basal crevasse (piezometric - height) basically
550       !From Faezeh:
551       !C_surface = LongitudinalDevStress(Rxx) - Rho*g*NodeDepth + WaterDepth*RhoW*g
552
553       !C_basal = LongitudinalDevStress(Rxx)-Rho*g*(IceThickness-HeightAboveGlacierBase)+RhoOcean*g*(PiezometricHeight-HeightAboveGlacierBase)
554       !Water pressure (Rho_w * g * Piezo) is PressureSol,PressureValues,PressurePerm
555
556       IF(PlaneStressOnly) THEN
557
558         CalvingSurfIndex = ( STuningParam * 2.0 * NodeStress1 ) - &
559              (NodeDepth * g * Rho) + (WaterDepth * RhoWF * g) - YieldStress
560         CalvingBasalIndex = (  BTuningParam * 2.0 * NodeStress1 ) - &
561              (NodeDepth * g * Rho) + NodeWPressure - YieldStress
562         !Last term appears if you split the brackets on the last term in C_Basal from Faezeh
563       ELSE
564         !Shear stress for Te calc
565         NodeStress4 = Stress4Values(Stress4Perm(i))
566         !2(Te^2) = Txx^2 + Tyy^2 + 2*Txy^2
567         !Txx = Tyy thus:
568         !Te^2 = (Txx^2) + (Txy^2)
569         !Te = ((Txx^2) + (Txy^2)) ^ 0.5
570         Te = ((NodeStress1**2) + (NodeStress4**2)) ** 0.5
571         sign = NodeStress1/abs(NodeStress1)
572
573         CalvingSurfIndex = ( STuningParam * 2.0 * sign * Te ) - &
574              (NodeDepth * g * Rho) + (WaterDepth * RhoWF * g) - YieldStress
575         CalvingBasalIndex = (  BTuningParam * 2.0 * sign * Te ) - &
576              (NodeDepth * g * Rho) + NodeWPressure - YieldStress
577       END IF
578
579     END IF
580
581     CSurfIndexValues(CSurfIndexPerm(i)) = CalvingSurfIndex
582     CBasalIndexValues(CBasalIndexPerm(i)) = CalvingBasalIndex
583   END DO
584
585   rt = RealTime() - rt0
586   IF(ParEnv % MyPE == 0) &
587        PRINT *, 'Calving, time taken for evaluating CIndex:', rt
588   rt0 = RealTime()
589
590   !---------------------------------------
591   !
592   ! Find connected crevassed regions and check for calving
593   !
594   !---------------------------------------
595
596   !TODO, allow both
597   IF(BasalCrevasseModel) THEN
598     CrevasseGroupValues = 0
599
600     !Find groups of nodes which have surface crevasses
601     CIndexValues => CSurfIndexValues
602     CIndexPerm => CSurfIndexPerm
603     CALL FindCrevasseGroups(SurfaceCrevasseGroups,.TRUE., TopPerm)
604
605     !As above, basal crevasses
606     CIndexValues => CBasalIndexValues
607     CIndexPerm => CBasalIndexPerm
608     CALL FindCrevasseGroups(BasalCrevasseGroups,.TRUE., BotPerm)
609
610     !Look for touching/almost touching crevasse groups
611     CALL FindCalvingBasal(SurfaceCrevasseGroups, BasalCrevasseGroups, BasalCalvingCoordinate)
612
613   ELSE !Surface Crevasse Model
614
615     !This dictates the resolution of the mesh for interpolating calving index
616     !Only relevant for surface crevasse mode, as opposed to surf and basal
617     !TODO: Unhardcode this
618     EvalResolution = 1.0_dp
619
620     !Create the points (at the waterline) at which Calving Index will be evaluated
621     !and assign them to a new mesh for interpolation
622
623     NofEvalPoints = FLOOR(FreeConnect/EvalResolution)
624     EvalMesh => AllocateMesh()
625     EvalMesh % Name = "Eval_Mesh"
626     EvalMesh % NumberOfNodes = NofEvalPoints
627     EvalMesh % Nodes => EvalNodes
628     ALLOCATE(Field(NofEvalPoints), &
629          FieldPerm(NofEvalPoints), &
630          EvalPoints(NofEvalPoints,2), &
631          EvalNodes % x(NofEvalPoints), &
632          EvalNodes % y(NofEvalPoints), &
633          EvalNodes % z(NofEvalPoints))
634
635     Field = -1.0_dp
636     EvalPoints(:,2) = SeaLevel
637     DO i=1,NofEvalPoints
638       EvalPoints(i,1) = Length - i*EvalResolution
639       FieldPerm(i) = i
640     END DO
641     EvalMesh % Nodes % x = EvalPoints(:,1)
642     EvalMesh % Nodes % y = SeaLevel
643     EvalMesh % Nodes % z = 0.0_dp
644
645     CALL VariableAdd( EvalMesh % Variables, EvalMesh, Solver, "Calving Surface Index", 1, Field, FieldPerm )
646     !TEST - should be true
647     CALL InterpolateMeshToMesh( Mesh, EvalMesh, Mesh % Variables, EvalMesh % Variables, .FALSE. )
648
649     CalvingOccurs = .FALSE.
650     Var => VariableGet( EvalMesh % Variables, 'Calving Surface Index', UnfoundFatal=.TRUE.)
651     DO i=1,NofEvalPoints
652       IF(Var % Values(Var % Perm(i)) < 0.0_dp ) CYCLE
653       CalvingOccurs = .TRUE.
654       PRINT *,'Surface Calving Event'
655       IF(EvalMesh % Nodes % x(i) < CalvingCoordHolder) &
656            CalvingCoordHolder = EvalMesh % Nodes % x(i)
657     END DO
658
659   END IF !Surface crevasse model
660
661   IF(BasalCrevasseModel) THEN
662      CalvingOccurs = .FALSE.
663      IF(BasalCalvingOccurs) THEN
664         CalvingOccurs = .TRUE.
665         DO j = 1, NoNodes
666            IF(Permutation(j) == 0) CYCLE
667            Calving1Values(Permutation(j)) = BasalCalvingCoordinate - &
668                 Mesh % Nodes % x(j)
669            IF(Calving1Values(Permutation(j)) .GE. 0.0_dp) &
670                 Calving1Values(Permutation(j)) = 0.0_dp
671         END DO
672         PRINT *, 'Basal Calving Event at timestep ',t
673         PRINT *, 'Calving Coordinate :',BasalCalvingCoordinate
674      ELSE
675         Calving1Values = 0.0_dp
676      END IF
677   ELSE
678      !Surface-only calving model
679      IF(CalvingOccurs) THEN
680         DO j = 1, NoNodes
681            IF(Permutation(j) == 0) CYCLE
682            Calving1Values(Permutation(j)) = CalvingCoordHolder - &
683                 Mesh % Nodes % x(j)
684            IF(Calving1Values(Permutation(j)) .GE. 0.0_dp) &
685                 Calving1Values(Permutation(j)) = 0.0_dp
686         END DO
687         PRINT *, 'Calving Event at timestep ',t
688         PRINT *, 'Calving Coordinate :',CalvingCoordHolder
689      ELSE
690         Calving1Values = 0.0_dp
691      END IF
692   END IF
693
694   rt = RealTime() - rt0
695   IF(ParEnv % MyPE == 0) &
696        PRINT *, 'Calving, time taken for calving connectivity:', rt
697   rt0 = RealTime()
698
699   !---------------------------------------
700   !   CALVING DONE
701   !---------------------------------------
702
703   !At this point, the 'calving' solution is done.  However...
704   !Here we solve a problem to do with undercutting:
705   !Progressive undercutting by melting of the calving front
706   !can lead to a situation where 'Front' nodes start to look like
707   !basal nodes.  However, they are 'officially' front nodes and so
708   !don't have a friction law, grounding dynamics OR (most importantly
709   !I think), a bed constraint.  Thus, it is necessary to check for this
710   !occurring, and shift the bed nodes appropriately.
711   !
712   !Strategy:
713   !-Identify the corner node by BotPerm and FrontPerm
714   !
715   !-Use BCelement connections to find the second to bottom
716   !    node on the calving front
717   !
718   !-Check a condition: either
719   !    second node is 'near' bed OR
720   !    BCelement slope is below some critical level
721   !
722   !-If condition is met (i.e. we need to take action)
723   !    Calving1Values @CornerNode = (X@2nd - X@corner)
724   !    CalvingOccurs = .TRUE.  <-- will this have any unforeseen consequences?
725   !
726   !NOTE: This works in tandem with a section of TwoMeshes.f90 which does the actual
727   !deformation
728
729   !Get the node index of the bottom corner
730   !NOTE: this could be 'FirstTime'd if it was also 'SAVE'd
731   DO i=1,NoNodes
732      IF(BotPerm(i) > 0 .AND. FrontPerm(i) > 0) THEN
733         BotCornerIndex = i
734      END IF
735   END DO
736
737   !Loop boundary elements, we're looking for the BCelement
738   !containing BotCornerIndex and ANOTHER FrontPerm node...
739   DO i=Mesh % NumberOfBulkElements+1,&
740        Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
741      CurrentElement => Mesh % Elements(i)
742      IF(.NOT.(ANY(CurrentElement % nodeindexes == BotCornerIndex))) CYCLE
743      IF(ALL(FrontPerm(CurrentElement % nodeindexes) .GT. 0)) THEN
744         !We have a winner
745         CornerElement => Mesh % Elements(i)
746         DO j=1,2
747            IF(CurrentElement % NodeIndexes(j) .NE. BotCornerIndex) THEN
748               BotSecondIndex = CurrentElement % NodeIndexes(j)
749            END IF
750         END DO
751      END IF
752   END DO
753
754   !Check corner node isn't already calving
755   IF(Calving1Values(Permutation(BotCornerIndex)) .LT. 0.0_dp) THEN
756      CornerCalving = .TRUE.
757   ELSE
758      CornerCalving = .FALSE.
759   END IF
760
761   !Get normal vector:
762   CALL GetElementNodes(CornerElementNodes, CornerElement)
763   CornerNormal = NormalVector(CornerElement, CornerElementNodes)
764
765   IF(Debug) PRINT *, 'Debug Calving, corner normal is: ' , &
766        CornerNormal(1), CornerNormal(2), CornerNormal(3)
767
768   IF(BasalFS) THEN
769     BedSecond = ListGetRealAtNode( Material, "Min "//BasalFSVarName, &
770          BotSecondIndex, UnfoundFatal=.TRUE. )
771
772     IF(Debug) PRINT *, 'Debug Calving, second node bed is: ',&
773          BedSecond,' and y coord is: ', Model % Nodes % y(BotSecondIndex)
774
775     PRINT *, 'Debug Calving, second node bed is: ',&
776          BedSecond,' and y coord is: ', Model % Nodes % y(BotSecondIndex)
777
778     BedSecondDiff = Model % Nodes % y(BotSecondIndex) - BedSecond
779   END IF
780
781   !TODO - unhardcode these
782   BedToler = 2.0_dp
783   normalcond = 0.95_dp
784
785   CornerBadBed = BasalFS .AND. (BedSecondDiff < BedToler)
786   CornerBadSlope = ABS(CornerNormal(2)) > normalcond
787
788   !If the slope normal is above threshold, or the second node is too close to the bed,
789   !move the corner node to the second, via 'calving'
790   IF((CornerBadSlope .OR. CornerBadBed) .AND. (.NOT. CornerCalving)) THEN
791
792      IF(Debug) PRINT *,'Debug Calving, migrating mesh'
793      county = 1
794      GoToNode = BotSecondIndex
795      PrevNode = BotCornerIndex
796      KeepLooking = .TRUE.
797      DO WHILE (KeepLooking)
798         !Check if we should shift more than one node forward...
799         IF(Debug) PRINT *, 'Debug calving: looking!'
800         KeepLooking = .FALSE.
801         DO i=Mesh % NumberOfBulkElements+1,&
802              Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements
803            CurrentElement => Mesh % Elements(i)
804            IF(.NOT.(ALL(FrontPerm(CurrentElement % nodeindexes) .GT. 0))) CYCLE
805            !This point reached by Front BC elements, stupid way of doing it, but whatever
806            IF(.NOT.(ANY(CurrentElement % nodeindexes == GoToNode))) CYCLE
807            IF(ANY(CurrentElement % nodeindexes == PrevNode)) CYCLE
808            !We only get here if element is the next one along from previous
809            CALL GetElementNodes(CurrentElementNodes, Currentelement)
810            Normal = NormalVector(CurrentElement, CurrentElementNodes)
811            DO j=1,2
812              IF(CurrentElement % NodeIndexes(j) .NE. GoToNode) &
813                   NextNode = CurrentElement % NodeIndexes(j)
814            END DO
815
816            IF(BasalFS) THEN
817              beddiff = Model % Nodes % y(NextNode) - ListGetRealAtNode( Material, &
818                   "Min "//BasalFSVarName, NextNode, UnfoundFatal=.TRUE. )
819            END IF
820
821            CornerBadBed = BasalFS .AND. (beddiff < BedToler)
822            CornerBadSlope = ABS(CornerNormal(2)) > normalcond
823
824            IF(CornerBadBed .OR. CornerBadSlope) THEN
825               PrevNode = GoToNode
826               GoToNode = NextNode
827               county = county + 1
828               IF(Debug) PRINT *, 'Debug calving: Found another shift'
829               KeepLooking = .TRUE.
830               EXIT
831            END IF
832         END DO
833      END DO
834
835      RemeshOccurs = .TRUE.
836   ELSE
837      county = 0
838   END IF
839
840   rt = RealTime() - rt0
841   IF(ParEnv % MyPE == 0) &
842        PRINT *, 'Calving, time taken for corner problems:', rt
843   rt0 = RealTime()
844
845   IF(CalvingOccurs .OR. RemeshOccurs) THEN
846
847      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
848      !Find the FrontDisplacement (= Calving 2 <-sif) for each frontal node
849      !resulting from the shift in the corner node
850      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
851
852      !Use FrontPerm to construct ordered node list
853      !I think MakeMaskUsingPerm already ordered the nodes, from the comments:
854      !! The bandwidth optimization for lines results to perfectly ordered
855      !! permutations. If there is only one line the 1st node should be the
856      !! lower left corner.
857      ALLOCATE(InvFrontPerm(FrontNodes),&
858           CumDist(FrontNodes),&
859           PropCumDist(FrontNodes),&
860           TargetCumDist(FrontNodes),&
861           TargetPropDist(FrontNodes),&
862           TargetNodes % x(FrontNodes),&
863           TargetNodes % y(FrontNodes),&
864           TargetNodes % z(FrontNodes),&
865           CalvedNodes % x(FrontNodes))
866
867      !InvFrontPerm(FrontNodes) points to nodeindexes in order they appear...
868      !InvFrontPerm(1) points to 'lower left' corner, according to MakePermUsingMask
869      DO i=1,NoNodes
870         IF(FrontPerm(i) .GT. 0) THEN
871            InvFrontPerm(FrontPerm(i)) = i
872         END IF
873      END DO
874
875      ALLOCATE(OrderPerm(FrontNodes), WorkReal(FrontNodes))
876      OrderPerm = [(i,i=1,FrontNodes)]
877      DO i=1,FrontNodes
878         WorkReal(i) = Mesh0 % Nodes % y(InvFrontPerm(i))
879      END DO
880      CALL SortD( FrontNodes, WorkReal, OrderPerm )
881      DEALLOCATE(WorkReal)
882
883      DO i=1,FrontNodes
884         j = InvFrontPerm(OrderPerm(i))
885         CalvedNodes % x(i) = Mesh % Nodes % x(j) + Calving1Values(Permutation(j))
886         IF(Debug) PRINT *,'Debug Calving, CalvedNodes node: ',j,' is ',&
887              CalvedNodes % x(i),' init coord: ',Mesh % Nodes % x(j)
888      END DO
889
890      !cycle through in order
891      !First get target distribution from Mesh0
892      TargetCumDist(1) = 0.0_dp
893      DO i=2,FrontNodes
894        dx = Mesh0 % Nodes % x(InvFrontPerm(OrderPerm(i))) - Mesh0 % Nodes % x(InvFrontPerm(OrderPerm(i-1)))
895        dy = Mesh0 % Nodes % y(InvFrontPerm(OrderPerm(i))) - Mesh0 % Nodes % y(InvFrontPerm(OrderPerm(i-1)))
896
897        TargetCumDist(i) = TargetCumDist(i-1) + (((dx**2) + (dy**2)) ** 0.5_dp)
898      END DO
899      TargetPropDist = TargetCumDist / MAXVAL(TargetCumDist)
900
901      !Now find the length segments of our current line
902      !If RemeshOccurs (because of bad corner node), county dictates
903      !the offset from the previous bottom node of the new calving front
904      CumDist(1:county+1) = 0.0_dp
905      DO i=county+2,FrontNodes
906        !sum coord magnitude from base upwards to give front 'length'
907        !keep cumulative total
908        !allocate proporitional y (and x) distances (i.e. out of 1)
909        dx = CalvedNodes % x(i) - CalvedNodes % x(i-1)
910        dy = Mesh % Nodes % y(InvFrontPerm(OrderPerm(i))) - Mesh % Nodes % y(InvFrontPerm(OrderPerm(i-1)))
911
912        CumDist(i) = CumDist(i-1) + (((dx**2) + (dy**2)) ** 0.5_dp)
913        !Remember first one is corner node...
914        IF(Debug) PRINT *, 'Debug Calving: CumDist at node: ',&
915             InvFrontPerm(OrderPerm(i)),' is ',CumDist(i)
916        IF(Debug) PRINT *, 'Debug Calving: TargetDist at node: ',&
917             InvFrontPerm(OrderPerm(i)),' is ',TargetCumDist(i)
918      END DO
919      PropCumDist = CumDist / MAXVAL(CumDist)
920
921      !Loop each front node
922      TargetNodes % x(1) = CalvedNodes % x(county+1)
923      TargetNodes % y(1) = Mesh % Nodes % y(InvFrontPerm(OrderPerm(county+1)))
924      TargetNodes % x(FrontNodes) = CalvedNodes % x(FrontNodes)
925      TargetNodes % y(FrontNodes) = Mesh % Nodes % y(InvFrontPerm(OrderPerm(FrontNodes)))
926
927      DO i=2,FrontNodes-1
928        !and find nearest two nodes to interpolate
929        DO j=county+2,FrontNodes
930          IF(PropCumDist(j) .GT. TargetPropDist(i)) THEN
931            !lin interp between j and j-1
932            LocalDist = PropCumDist(j) - PropCumDist(j-1)
933            LocalDistNode = TargetPropDist(i) - PropCumDist(j-1)
934
935            PropDistNode = LocalDistNode / LocalDist
936            IF(Debug) PRINT *, 'Debug Calving: PropDist at node: ',&
937                 InvFrontPerm(OrderPerm(i)),' is ',PropDistNode
938
939            TargetNodes % x(i) = ((1 - PropDistNode) * CalvedNodes % x(j-1))  + &
940                 (PropDistNode * CalvedNodes % x(j))
941            TargetNodes % y(i) = ((1 - PropDistNode) * Mesh % Nodes % y(InvFrontPerm(OrderPerm(j-1))))  + &
942                 (PropDistNode * Mesh % Nodes % y(InvFrontPerm(OrderPerm(j))))
943            EXIT
944          END IF
945        END DO
946      END DO
947
948      !At this point, we have obtained, for each FrontNode, a TargetNode % x and y
949      !Thus, it simply remains to compute the two components of the displacement
950
951      !Calving 1 = Diff X  (New % x - Old % x)
952      !Calving 2 = Diff Y  (New % y - Old % y)
953
954      DO i=1,FrontNodes
955         Calving1Values(Permutation(InvFrontPerm(OrderPerm(i)))) = TargetNodes % x(i) &
956              - Mesh % Nodes % x(InvFrontPerm(OrderPerm(i)))
957
958         Calving2Values(Permutation(InvFrontPerm(OrderPerm(i)))) = TargetNodes % y(i) &
959              - Mesh % Nodes % y(InvFrontPerm(OrderPerm(i)))
960
961         IF(Debug) THEN
962            PRINT *,'Debug Calving: Node: ',InvFrontPerm(OrderPerm(i)),' pos x: ',&
963                 Mesh % nodes % x(InvFrontPerm(OrderPerm(i))),&
964                 ' pos y: ',Mesh % nodes % y(InvFrontPerm(OrderPerm(i)))
965            PRINT *,'Moving to: x: ',TargetNodes % x(i),' y: ',TargetNodes % y(i)
966            PRINT *,'Displacement 1: ',Calving1Values(Permutation(InvFrontPerm(OrderPerm(i)))),&
967                 'Displacement 2: ',Calving2Values(Permutation(InvFrontPerm(OrderPerm(i))))
968         END IF
969       END DO
970     END IF
971
972   CALL ListAddLogical( Model % Simulation, 'CalvingOccurs', CalvingOccurs )
973   CALL ListAddLogical( Model % Simulation, 'RemeshOccurs', RemeshOccurs )
974
975   rt = RealTime() - rt0
976   IF(ParEnv % MyPE == 0) &
977        PRINT *, 'Calving, time taken for calculating displacements (end):', rt
978   rt0 = RealTime()
979
980 CONTAINS
981
982   SUBROUTINE FindNodeNeighbours(NodeNumber)
983     INTEGER :: NodeNumber, i, count
984
985     NoNeighbours = Vel1Matrix % Rows((Vel1Perm(NodeNumber)*DOFs)+1) - Vel1Matrix % Rows(Vel1Perm(NodeNumber)*DOFs)
986     IF(MOD(NoNeighbours, DOFs).NE. 0) CALL Fatal(SolverName,"This shouldn't have happened...")
987     !Each neighbour appears once per DOF, and there's also the current node thus: (x/DOFS) - 1...
988     NoNeighbours = (NoNeighbours / DOFs) - 1
989     IF(NoNeighbours .GT. MaxNeighbours) CALL Fatal(SolverName,"Need more neighbour slots!")
990
991     IF(ALLOCATED(ThisNodeNeighbours)) DEALLOCATE(ThisNodeNeighbours)
992     ALLOCATE(ThisNodeNeighbours(NoNeighbours))
993     ThisNodeNeighbours = 0
994
995     count = 0
996     DO i=Vel1Matrix % Rows(Vel1Perm(NodeNumber)*DOFs),&
997          (Vel1Matrix % Rows((Vel1Perm(NodeNumber)*DOFs)+1)-1)
998        IF(MOD(i,DOFs) .NE. 0) CYCLE !Stored DOF1, DOF2, DOF3, only need every 3rd
999        IF(Vel1InvPerm(Vel1Matrix % Cols(i)) == NodeNumber) CYCLE !Not our own neighbour
1000        count = count + 1
1001        ThisNodeNeighbours(count) = &
1002             Vel1InvPerm(Vel1Matrix % Cols(i))
1003     END DO
1004
1005   END SUBROUTINE FindNodeNeighbours
1006
1007   !After finding groups, pass them between partitions and join...
1008   !Can just be done in one partition and then passed back?
1009   SUBROUTINE FindCrevasseGroups(CrevasseGroups, CheckValidity, MaskPerm)
1010
1011     TYPE(CrevasseGroups_t), INTENT(INOUT) :: CrevasseGroups
1012     INTEGER :: MaxGroups = 100, MaxNodesPerGroup = 10000, locali
1013     LOGICAL :: CheckValidity, FoundIt
1014     INTEGER, POINTER, OPTIONAL, INTENT(IN) :: MaskPerm(:)
1015
1016     ALLOCATE( &
1017          CrevasseGroups % NodeIndexes(MaxGroups,MaxNodesPerGroup), &
1018          CrevasseGroups % NotEmpty(MaxGroups), &
1019          CrevasseGroups % NoNodes(MaxGroups), &
1020          CrevasseGroups % Valid(MaxGroups))
1021
1022     CrevasseGroups % NodeIndexes = 0
1023     CrevasseGroups % NoNodes = 0
1024     CrevasseGroups % NotEmpty = .FALSE.
1025     CrevasseGroups % Valid = .FALSE.
1026     CrevasseGroups % CurrentGroup = 0
1027
1028     DO i=1,NoNodes
1029        IF(DistanceValues(DistancePerm(i))>FreeConnect) CYCLE
1030        IF(CIndexValues(CIndexPerm(i))<0) CYCLE
1031        IF(FrontPerm(i) > 0) CYCLE
1032
1033        !Check if node is already in a group
1034        !This aint ideal, but I put it in to prevent a bug whereby, when surface and
1035        !basal groups join, all surface group nodes are ALSO basal group nodes, and so
1036        !'calving' is detected everywhere...
1037        FoundIt = .FALSE.
1038        DO locali = 1, SurfaceCrevasseGroups % CurrentGroup
1039          IF(ANY(SurfaceCrevasseGroups % NodeIndexes( &
1040               locali,1:SurfaceCrevasseGroups % NoNodes(locali)) .EQ. i)) THEN
1041            FoundIt = .TRUE.
1042            EXIT
1043          END IF
1044        END DO
1045        !Only check basal groups if they exist already...
1046        IF(ALLOCATED(BasalCrevasseGroups % NodeIndexes)) THEN
1047          DO locali = 1, BasalCrevasseGroups % CurrentGroup
1048            IF(ANY(BasalCrevasseGroups % NodeIndexes( &
1049                 locali,1:BasalCrevasseGroups % NoNodes(locali)) .EQ. i)) THEN
1050              FoundIt = .TRUE.
1051              EXIT
1052            END IF
1053          END DO
1054        END IF
1055        IF(FoundIt) CYCLE
1056
1057        !This point is only reached by nodes which ARE crevassing
1058        !and not already in a group, so it takes the first slot of
1059        !active group
1060        CrevasseGroups % CurrentGroup = CrevasseGroups % CurrentGroup + 1
1061        CrevasseGroups % NodeIndexes(CrevasseGroups % CurrentGroup,1) = i
1062        CrevasseGroups % NotEmpty(CrevasseGroups % CurrentGroup) = .TRUE.
1063        CrevasseGroups % NoNodes(CrevasseGroups % CurrentGroup) = &
1064             CrevasseGroups % NoNodes(CrevasseGroups % CurrentGroup) + 1
1065        CrevasseGroupValues(CrevasseGroupPerm(i)) = CrevasseGroups % CurrentGroup
1066        NextSlot = 2
1067        CALL SearchNeighbours(i, CrevasseGroups)
1068        NextSlot = 1
1069     END DO
1070
1071     !Now we have all the groups.  If requested, check they are valid.
1072     IF(CheckValidity) THEN
1073        DO i=1,CrevasseGroups % CurrentGroup !Cycle all groups
1074           DO j=1,SIZE(CrevasseGroups % NodeIndexes,2) !Cycle all nodes in group
1075              IF(CrevasseGroups % NodeIndexes(i,j)==0) EXIT !End of group
1076              IF(MaskPerm(CrevasseGroups % NodeIndexes(i,j))>0) THEN
1077                 !At least 1 node in group is on relevant surface
1078                 CrevasseGroups % Valid(i) = .TRUE.
1079                 EXIT
1080              END IF
1081           END DO
1082        END DO
1083     END IF
1084
1085   END SUBROUTINE FindCrevasseGroups
1086
1087   !Recursively looks in neighbouring nodes to construct contiguous groups of
1088   !crevassing nodes.
1089   RECURSIVE SUBROUTINE SearchNeighbours(nodenum, CrevasseGroups)
1090     INTEGER :: nodenum, neighbourindex
1091     INTEGER :: localj,locali
1092     TYPE(CrevasseGroups_t) :: CrevasseGroups
1093     LOGICAL :: FoundIt
1094
1095     IF(Debug) PRINT *,'Debug, nextslot, nodenum: ', nextslot, nodenum
1096     DO localj = 1,NumNeighbours(nodenum)
1097        neighbourindex = NodeNeighbours(nodenum,localj)
1098        IF(DistanceValues(DistancePerm(neighbourindex))>FreeConnect) CYCLE
1099        IF(CIndexValues(CIndexPerm(neighbourindex))<0) CYCLE
1100        IF(FrontPerm(localj) > 0) CYCLE
1101
1102        !Check if node is already in a group
1103        !NOTE: is this actually possible?
1104        FoundIt = .FALSE.
1105        DO locali = 1, SurfaceCrevasseGroups % CurrentGroup
1106          IF(ANY(SurfaceCrevasseGroups % NodeIndexes(&
1107               locali,1:SurfaceCrevasseGroups % NoNodes(locali)) &
1108               .EQ. neighbourindex)) THEN
1109            FoundIt = .TRUE.
1110            EXIT
1111          END IF
1112        END DO
1113        IF(ALLOCATED(BasalCrevasseGroups % NodeIndexes)) THEN
1114          DO locali = 1, BasalCrevasseGroups % CurrentGroup
1115            IF(ANY(BasalCrevasseGroups % NodeIndexes(&
1116                 locali,1:BasalCrevasseGroups % NoNodes(locali)) &
1117                 .EQ. neighbourindex)) THEN
1118              FoundIt = .TRUE.
1119              EXIT
1120            END IF
1121          END DO
1122        END IF
1123        IF(FoundIt) CYCLE
1124
1125        !Only get here if node is valid new calving node, not already in a group
1126        CrevasseGroups % NodeIndexes(CrevasseGroups % CurrentGroup, NextSlot) = neighbourindex
1127        CrevasseGroups % NoNodes(CrevasseGroups % CurrentGroup) = &
1128             CrevasseGroups % NoNodes(CrevasseGroups % CurrentGroup) + 1
1129        CrevasseGroupValues(CrevasseGroupPerm(neighbourindex)) = CrevasseGroups % CurrentGroup
1130        NextSlot = NextSlot + 1
1131        IF(NextSlot > SIZE(CrevasseGroups % NodeIndexes, 2)) CALL Fatal(SolverName, &
1132             "More than 10000 nodes in a crevasse group? Almost certainly an error in setup")
1133
1134        CALL SearchNeighbours(neighbourindex,CrevasseGroups)
1135        !Add a check for group validity (i.e. at least one node is on the relevant boundary
1136     END DO
1137   END SUBROUTINE SearchNeighbours
1138
1139   SUBROUTINE FindCalvingBasal(SurfaceCrevasseGroups, BasalCrevasseGroups, BasalCalvingCoordinate)
1140     TYPE(CrevasseGroups_t), INTENT(IN) :: SurfaceCrevasseGroups, BasalCrevasseGroups
1141     REAL (KIND=dp), INTENT(OUT) :: BasalCalvingCoordinate
1142     INTEGER :: i,j,k,m, BasalNode, SurfaceNode
1143
1144     BasalCalvingOccurs = .FALSE.
1145     BasalCalvingCoordinate = HUGE(BasalCalvingCoordinate)
1146     !Cycle surface crevasse groups
1147     DO i = 1, COUNT(SurfaceCrevasseGroups % NotEmpty)
1148        IF(.NOT.(SurfaceCrevasseGroups % Valid(i))) CYCLE
1149
1150        !Cycle basal crevasse groups
1151        DO j = 1, COUNT(BasalCrevasseGroups % NotEmpty)
1152           IF(.NOT.(BasalCrevasseGroups % Valid(j))) CYCLE
1153
1154           !Cycle Nodes in Surface Crevasse Group
1155           DO k = 1,SIZE(SurfaceCrevasseGroups % NodeIndexes, 2)
1156              IF(SurfaceCrevasseGroups % NodeIndexes(i,k)==0) EXIT
1157              SurfaceNode = SurfaceCrevasseGroups % NodeIndexes(i,k)
1158
1159              !Cycle Nodes in Basal Crevasse Group
1160              DO m = 1,SIZE(BasalCrevasseGroups % NodeIndexes, 2)
1161                 IF(BasalCrevasseGroups % NodeIndexes(j,m)==0) EXIT
1162                 BasalNode = BasalCrevasseGroups % NodeIndexes(j,m)
1163                 !Is it the same node? i.e. do the groups touch?
1164                 IF(SurfaceNode == BasalNode) THEN
1165                    !Yes, so calving coord
1166                    IF(Mesh % Nodes % x(SurfaceNode) .LT. BasalCalvingCoordinate) THEN
1167                       BasalCalvingCoordinate = Mesh % Nodes % x(SurfaceNode)
1168                       BasalCalvingOccurs = .TRUE.
1169                       !TEST
1170                       IF(Debug) THEN
1171                          PRINT *, "Debugging Calving, Surface node x: ", &
1172                               Mesh % Nodes % x(SurfaceNode)," y: ",&
1173                               Mesh % Nodes % y(SurfaceNode)
1174                          PRINT *, "Debugging Calving, Basal node x: ", &
1175                               Mesh % Nodes % x(BasalNode)," y: ",&
1176                               Mesh % Nodes % y(BasalNode)
1177                       END IF
1178                    END IF
1179                 ELSE
1180                    !No, but are they neighbours?
1181                    IF(ANY(NodeNeighbours(BasalNode,:)==SurfaceNode)) THEN
1182                       !yes, neighbours
1183                       CALL CheckCIndexOverlap(SurfaceNode, BasalNode, OverlapCalvingCoordinate, OverlapOccurs)
1184                       IF(OverlapOccurs .AND. (OverlapCalvingCoordinate .LT. BasalCalvingCoordinate)) THEN
1185                          BasalCalvingCoordinate = OverlapCalvingCoordinate
1186                          BasalCalvingOccurs = .TRUE.
1187                       END IF
1188                    END IF
1189                    !not neighbours, cycle
1190                 END IF
1191
1192              END DO
1193           END DO
1194        END DO
1195     END DO
1196   END SUBROUTINE FindCalvingBasal
1197
1198   SUBROUTINE CheckCIndexOverlap(SurfaceNode, BasalNode, OverlapCoord, OverlapOccurs)
1199
1200     IMPLICIT NONE
1201     REAL(KIND=dp) :: CSurfSurf, CSurfBasal, CBasalSurf, &
1202          CBasalBasal, XSurf, YSurf, XBasal, YBasal, dx, dy,dxdy, dCSurf, &
1203          dCBasal, xzerobasal, yzerobasal, xzerosurf, yzerosurf, &
1204          dbindexdy,dsindexdy
1205     REAL(KIND=dp), INTENT(OUT) :: OverlapCoord
1206     INTEGER :: SurfaceNode, BasalNode
1207     LOGICAL, INTENT(OUT) :: OverlapOccurs
1208
1209     OverlapOccurs = .FALSE.
1210
1211
1212     !4 values at 2 nodes
1213     !Format is CIndexLocation.  So CSurfBasal is the value of
1214     !CSurfIndex at the Basal node.
1215
1216     !This used to be simple linear interpolation between two nodes, but this was problematic:
1217     !For surface crevassing nodes, basal crevasse index is also positive, because they are
1218     !governed by the same equations.  This lead to situations where both CBasalSurf and CBasalBasal
1219     !would be small positive, with an OBVIOUS massive gap inbetween, but because both were positive,
1220     !linear interp didn't see the gap.
1221
1222     !Strategy to overcome the previous:
1223     !No problem with CSurf* because it decreases with depth as expected.
1224     !Need to address CBasal*: at the lower node (C*Basal), assuming negligible upward changes in
1225     !stress_xx (any better way?), then the y-gradient of CBasal* is predictable, as the remaining
1226     !components decrease linearly with depth:
1227     ! rho_w * g * dy   -   BTuningParam * rho_i * g * dy
1228
1229     !So, we replace the part of this subroutine which used to calculate the zero level of basal
1230     !crevassing.  We need:
1231     !g, rho_i, rho_w, BTuningParam
1232
1233     !the various indices
1234     CSurfSurf = CSurfIndexValues(CSurfIndexPerm(SurfaceNode))
1235     CBasalSurf = CBasalIndexValues(CBasalIndexPerm(SurfaceNode))
1236     CSurfBasal = CSurfIndexValues(CSurfIndexPerm(BasalNode))
1237     CBasalBasal = CBasalIndexValues(CBasalIndexPerm(BasalNode))
1238
1239     !the coords
1240     XSurf = Mesh % Nodes % x(SurfaceNode)
1241     YSurf = Mesh % Nodes % y(SurfaceNode)
1242     XBasal = Mesh % Nodes % x(BasalNode)
1243     YBasal = Mesh % Nodes % y(BasalNode)
1244
1245     !The rest of the maths assumes that the surface node is ABOVE the basal node, so check:
1246     !If the node in the basal crevasse field is *above* the node in the surface crevasse field
1247     !we have calving and we assume its x-coord is between the two.  else check...
1248
1249     IF(YSurf .LE. YBasal) THEN
1250        IF(Debug) PRINT *, 'DEBUG Calving: Basal Node above Surf Node, calving...'
1251        OverlapOccurs = .TRUE.
1252        OverlapCoord = (XSurf + XBasal)/2.0
1253     ELSE
1254        !Gradients
1255        dx = XSurf - XBasal
1256        dy = YSurf - YBasal !+ve
1257        dxdy = dx/dy
1258
1259        dCSurf = CSurfSurf - CSurfBasal  !+ve
1260        dCBasal = CBasalSurf - CBasalBasal  !-ve
1261
1262        dsindexdy = dCSurf/dy !+ve
1263        dbindexdy = -g*(RhoWF - Rho) !-ve
1264
1265        yzerobasal = YBasal - (CBasalBasal/dbindexdy)
1266        yzerosurf = YSurf - (CSurfSurf/dsindexdy)
1267
1268        IF(yzerosurf .LT. yzerobasal) THEN
1269           OverlapOccurs = .TRUE.
1270
1271           xzerobasal = xbasal + ((yzerobasal - YBasal)*dxdy)
1272           xzerosurf = xsurf + ((yzerosurf - YSurf)*dxdy)
1273
1274           OverlapCoord = MIN(xzerosurf, xzerobasal)
1275           OverlapCoord = MAX(OverlapCoord,MIN(XSurf, XBasal))
1276        ELSE
1277           OverlapOccurs = .FALSE.
1278        ENDIF
1279     END IF
1280
1281     IF(OverlapOccurs .AND. Debug) THEN
1282        PRINT *, "Overlap occurs!, OverlapCoord: ",OverlapCoord
1283        PRINT *, "Surf: ",SurfaceNode," x: ",Mesh % Nodes % x(SurfaceNode),&
1284             "y: ",Mesh % Nodes % y(SurfaceNode)
1285        PRINT *, "Base: ",BasalNode," x: ",Mesh % Nodes % x(BasalNode),&
1286             "y: ",Mesh % Nodes % y(BasalNode)
1287
1288     END IF
1289
1290   END SUBROUTINE CheckCIndexOverlap
1291END SUBROUTINE Find_Calving
1292