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