!!!!!!!!!!!!!!!!!!! ! Returns the year index. !!------------------------------------------------------------------------------!! FUNCTION getYear(Time) RESULT(year) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE Real(KIND=dp) :: Time INTEGER :: monthIndex, getMonthIndex, year monthIndex = getMonthIndex(Time); IF(monthIndex <=12) THEN year = 1 ELSEIF(monthIndex > 12 .AND. monthIndex <= 24) THEN year = 2 ELSEIF(monthIndex > 24 .AND. monthIndex <= 36) THEN year = 3 ELSEIF(monthIndex > 36 .AND. monthIndex <= 48) THEN year = 4 ELSEIF(monthIndex > 48 .AND. monthIndex <= 60) THEN year = 5 ELSEIF(monthIndex > 60 .AND. monthIndex <= 72) THEN year = 6 ELSEIF(monthIndex > 72 .AND. monthIndex <= 84) THEN year = 7 ELSEIF(monthIndex > 84 .AND. monthIndex <= 96) THEN year = 8 ELSEIF(monthIndex > 96 .AND. monthIndex <= 108) THEN year = 9 ELSEIF(monthIndex > 108 .AND. monthIndex <= 120) THEN year = 10 ELSEIF(monthIndex > 120 .AND. monthIndex <= 132) THEN year = 11 ELSEIF(monthIndex > 132 .AND. monthIndex <= 144) THEN year = 12 ELSE CALL Fatal("getYear", "Time greater than 12 years") END IF END FUNCTION getYear !!!!!!!!!!!!!!!!!!! ! Returns the month index. !!------------------------------------------------------------------------------!! FUNCTION getMonthIndex(Time) RESULT(monthIndex) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE Real(KIND=dp) :: Time INTEGER :: monthIndex monthIndex = Time / (86400*30); END FUNCTION getMonthIndex !!!!!!!!!!!!!!!!!!! ! Checks whether to change simulation parameters. The function returns 1 if the given time ! is the first month of the new year, 0 if it is the middle of the year, -1 otherwise. !!------------------------------------------------------------------------------!! FUNCTION checkTime(Time, SwitchPeriod) RESULT(status) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model Real(KIND=dp) :: Time, SwitchPeriod LOGICAL :: Found INTEGER :: status, TimeInMonths, getMonthIndex, currentyear, TimeInMonths1 ! Time in months TimeInMonths = getMonthIndex(Time); currentyear = TimeInMonths / 12 TimeInMonths1 = TimeInMonths - (currentyear * 12) ! Check whether the new time instant is the beggining of a new year. ! Get the time and remove 1 month. If the rem is 0 it is the new year. ! Example 33696000 secs = 13*86400*365 => First month of new year. ! Returns status=1 IF (TimeInMonths1 == 1 .OR. MOD(TimeInMonths1 - 1, 12) == 0) THEN status = 1 ! Check if the new time instant is the middle of the year. Returns status=0 ELSE IF(MOD(TimeInMonths1, 7) == 0) THEN status = 0 ! Nothing to change ELSE status = -1 END IF END FUNCTION checkTime !!!!!!!!!!!!!!!!!!! ! Returns the flux balance from given files to be used in the free surface equation. !!------------------------------------------------------------------------------!! FUNCTION getBalance ( Model, nodenumber) RESULT(An) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver TYPE(Variable_t), POINTER :: TimeVar INTEGER :: nodenumber Real(KIND=dp), DIMENSION(:), ALLOCATABLE :: x, y, A INTEGER :: i, Nb, idx, timeStatus, checkTime, yearIndex, getYear CHARACTER(LEN=200) :: file, prevFile, WinterAccumulationFile, SummerAccumulationFile, BalanceFolder, format_string Real(KIND=dp) :: xn, yn, An, dist, distmin Real(KIND=dp) :: Time, SwitchPeriod LOGICAL :: Found LOGICAL :: FirstTime=.True. SAVE FirstTime, file SAVE x, y, A SAVE prevFile ! number of lines in balance files Nb = 336420 ! get current time TimeVar => VariableGet( Model % Variables, "Time") Time = TimeVar % Values(1) SwitchPeriod = GetConstReal(Model % Constants, "SwitchPeriod", Found) IF(.NOT.Found) CALL Fatal("getBalance", "Unable to find SwitchPeriod") BalanceFolder = GetString(Model % Constants, "BalanceFolder", Found) ! Use mean balance IF (.NOT.Found) THEN WinterAccumulationFile = GetString(Model % Constants, "WinterAccumulationFile", Found) IF(.NOT.Found) CALL Fatal("getBalance", "Unable to find WinterAccumulationFile") SummerAccumulationFile = GetString(Model % Constants, "SummerAccumulationFile", Found) IF(.NOT.Found) CALL Fatal("getBalance", "Unable to find SummerAccumulationFile") ! Use yearly balance ELSE yearIndex = getYear(Time); IF (yearIndex < 10) THEN format_string = "(A, I1, A)" ELSE format_string = "(A, I2, A)" END IF WRITE(WinterAccumulationFile, format_string) "winter-summer 1997-2007/winter_", yearIndex, ".dat" WRITE(SummerAccumulationFile, format_string) "winter-summer 1997-2007/summer_", yearIndex, ".dat" END IF ! Check whether to apply the accumulation timeStatus = checkTime(Time, SwitchPeriod); IF (timeStatus == 1) THEN file = WinterAccumulationFile ELSE IF (timeStatus == 0) THEN file = SummerAccumulationFile END IF IF(prevFile /= file) THEN FirstTime = .True. CALL INFO("getBalance", "Set accumulation file to " // file, Level=3) END IF prevFile = file ! get (x, y) for the point xn = Model % Nodes % x (nodenumber) yn = Model % Nodes % y (nodenumber) ! determine total number of lines in file IF (FirstTime) THEN IF (ALLOCATED(x)) DEALLOCATE(x) IF (ALLOCATED(y)) DEALLOCATE(y) IF (ALLOCATED(A)) DEALLOCATE(A) FirstTime = .False. OPEN(10, file=file) ALLOCATE(x(Nb), y(Nb), A(Nb)) READ(10, *)(x(i), y(i), A(i), i=1, Nb) CLOSE(10) END IF ! brutal search for nearest point idx = -1 DO i = 1, SIZE(x) dist = (x(i) - xn) ** 2 + (y(i) -yn) ** 2 IF ( idx .EQ. -1 .OR. dist .LT. distmin ) THEN idx = i distmin = dist END IF END DO An = A(idx); ! IF (An < 0) THEN ! An = 1.5 * An; ! END IF An = An / 86400 / 365 * 2 !OPEN(UNIT=12, FILE="aoutput.txt", ACTION="write", STATUS="old", POSITION="append") ! WRITE(12, *) " #", nodenumber, " => X: ", xn, " Y: ", yn, " A: ", A(idx), " An: ", An !WRITE(12, *) xn, yn, A(idx), An !CLOSE(12); END FUNCTION getBalance !!!!!!!!!!!!!!!!!!! ! Interpolates the node coordinates from DEM points. !!------------------------------------------------------------------------------!! FUNCTION InterpolateDEM (x, y, xb, yb, zb, Nbx, Nby, xb0, yb0, lbx, lby) RESULT(zbed) USE DefUtils IMPLICIT NONE INTEGER :: imin, Npt, t INTEGER :: NMAX, i, j, Nb, Nbx, Nby, ib, ix, iy REAL(KIND=dp) :: x, y, zbed, xb0, yb0, x1, x2, y1, y2, zi(2,2) REAL(KIND=dp) :: R, Rmin, lbx, lby, dbx, dby REAL(KIND=dp) :: xb(Nbx*Nby), yb(Nbx*Nby), zb(Nbx*Nby) ! Find zbed for that point from the Bedrock MNT dbx = lbx / (Nbx-1.0) dby = lby / (Nby-1.0) Nb = Nbx*Nby ix = INT((x-xb0)/dbx)+1 iy = INT((y-yb0)/dbx)+1 ib = Nbx * (iy - 1) + ix x1 = xb(ib) x2 = xb(ib+1) y1 = yb(ib) y2 = yb(ib + Nbx) zi(1,1) = zb(ib) zi(2,1) = zb(ib+1) zi(2,2) = zb(ib + Nbx + 1) zi(1,2) = zb(ib + Nbx) IF ((zi(1,1)<-9990.0).OR.(zi(1,2)<-9990.0).OR.(zi(2,1)<-9990.0).OR.(zi(2,2)<-9990.0)) THEN IF ((zi(1,1)<-9990.0).AND.(zi(1,2)<-9990.0).AND.(zi(2,1)<-9990.0).AND.(zi(2,2)<-9990.0)) THEN ! Find the nearest point avalable Rmin = 9999.0 DO i=1, Nb IF (zb(i)>0.0) THEN R = SQRT((x-xb(i))**2.0+(y-yb(i))**2.0) IF (R 0.0) THEN zbed = zbed + zi(i,j) Npt = Npt + 1 END IF END DO END DO zbed = zbed / Npt END IF ELSE zbed = (zi(1,1)*(x2-x)*(y2-y)+zi(2,1)*(x-x1)*(y2-y)+zi(1,2)*(x2-x)*(y-y1)+zi(2,2)*(x-x1)*(y-y1))/(dbx*dby) END IF END FUNCTION InterpolateDEM !!!!!!!!!!!!!!!!!!! ! Returns the elevation of the top surface from the DEM. Pay attention that if the elevation ! of the top surface is estimated from interpolation the mesh can get tangled (the elevation ! of the top node is smaller than the elavation of the bedrock). In this situation the top node ! will have the elavation of the bedrock plus MinIceThickness !!------------------------------------------------------------------------------!! FUNCTION getSurfaceElevation ( Model, nodenumber, znode) RESULT(Zsurf) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver INTEGER :: nodenumber REAL(KIND=dp) :: znode, Zbed, Zsurf INTEGER :: Nb, Nbx, Nby, i, j INTEGER :: Ns, Nsx, Nsy REAL(KIND=dp) :: x, y, MinIceThickness REAL(KIND=dp) :: xb0, yb0, xs0, ys0 REAL(KIND=dp) :: lbx, lby, lsx, lsy REAL(KIND=dp) :: InterpolateDEM REAL(KIND=dp), ALLOCATABLE :: xb(:), yb(:), zb(:) REAL(KIND=dp), ALLOCATABLE :: xs(:), ys(:), zs(:) CHARACTER(LEN=100) :: SurfaceDEM, BedrockDEM LOGICAL :: Found LOGICAL :: FirstTime=.True. SAVE FirstTime SAVE xb, yb, zb, xs, ys, zs MinIceThickness = GetConstReal(Model % Constants, "MinIceThickness", Found) IF(.NOT.Found) CALL Fatal("getSurfaceElevation", "Unable to find MinIceThickness") SurfaceDEM = GetString(Model % Constants, "SurfaceDEM", Found) IF(.NOT.Found) CALL Fatal("getSurfaceElevation", "Unable to find SurfaceDEM") BedrockDEM = GetString(Model % Constants, "BedrockDEM", Found) IF(.NOT.Found) CALL Fatal("getSurfaceElevation", "Unable to find BedrockDEM") Nsx = 801 Nsy = 420 xs0 = 610010.31250 ys0 = 5108310.50000 lsx = 24000 lsy = 12570 Ns = Nsx * Nsy Nbx = 801 Nby = 420 xb0 = 610010.31250 yb0 = 5108310.50000 lbx = 24000 lby = 12570 Nb = Nbx * Nby IF (FirstTime) THEN FirstTime = .False. OPEN(10,file=SurfaceDEM) ALLOCATE(xs(Ns), ys(Ns), zs(Ns)) READ(10,*)(xs(i), ys(i), zs(i), i=1,Ns) CLOSE(10) OPEN(10,file=BedrockDEM) ALLOCATE(xb(Nb), yb(Nb), zb(Nb)) READ(10,*)(xb(j), yb(j), zb(j), j=1,Nb) CLOSE(10) END IF x = Model % Nodes % x (nodenumber) y = Model % Nodes % y (nodenumber) Zbed = InterpolateDEM(x,y,xb,yb,zb,Nbx,Nby,xb0,yb0,lbx,lby) Zsurf = InterpolateDEM(x,y,xs,ys,zs,Nsx,Nsy,xs0,ys0,lsx,lsy) ! Check that mesh don't get tangled IF (Zsurf-Zbed < MinIceThickness) THEN Zsurf = Zbed + MinIceThickness END IF END FUNCTION getSurfaceElevation !!!!!!!!!!!!!!!!!!! ! Returns the bedrock altitude for a given nodes !!------------------------------------------------------------------------------!! FUNCTION getBedrockElevation ( Model, nodenumber, znode) RESULT(Zbed) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver INTEGER :: nodenumber REAL(KIND=dp) :: znode, Zbed INTEGER :: imin, Npt, t INTEGER :: NMAX, i, j, Nb, Nbx, Nby REAL(KIND=dp) :: x, y, z, xb0, yb0 REAL(KIND=dp) :: R, Rmin, lbx, lby, InterpolateDEM REAL(KIND=dp), ALLOCATABLE :: xb(:), yb(:), zb(:) CHARACTER(LEN=100) :: BedrockDEM LOGICAL :: Found, FirstTime=.True. SAVE FirstTime SAVE xb, yb, zb Nbx = 801 Nby = 420 xb0 = 610010.31250 yb0 = 5108310.50000 lbx = 24000 lby = 12570 Nb = Nbx * Nby BedrockDEM = GetString(Model % Constants, "BedrockDEM", Found) IF(.NOT.Found) CALL Fatal("getBedrockElevation", "Unable to find BedrockDEM") IF (FirstTime) THEN FirstTime = .False. OPEN(10,file=BedrockDEM) ALLOCATE(xb(Nb), yb(Nb), zb(Nb)) READ(10,*)(xb(i), yb(i), zb(i), i=1,Nb) CLOSE(10) END IF ! zbed for point (x,y) x = Model % Nodes % x (nodenumber) y = Model % Nodes % y (nodenumber) zbed = InterpolateDEM(x,y,xb,yb,zb,Nbx,Nby,xb0,yb0,lbx,lby) END FUNCTION getBedrockElevation !!!!!!!!!!!!!!!!!!! ! Returns the minimum elevation of the free surface. The top node has the elavation of the bedrock ! plus MinIceThickness. !!------------------------------------------------------------------------------!! FUNCTION getMinSurfaceElevation ( Model, nodenumber, znode) RESULT(Zsurf) USE types USE CoordinateSystems USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver INTEGER :: nodenumber REAL(KIND=dp) :: znode, Zbed, Zsurf INTEGER :: imin, Npt, t INTEGER :: NMAX, i, j, Nb, Nbx, Nby REAL(KIND=dp) :: x, y, z, xb0, yb0, MinIceThickness REAL(KIND=dp) :: R, Rmin, lbx, lby, InterpolateDEM REAL(KIND=dp), ALLOCATABLE :: xb(:), yb(:), zb(:) CHARACTER(LEN=100) :: BedrockDEM LOGICAL :: Found LOGICAL :: FirstTime=.True. SAVE FirstTime SAVE xb, yb, zb MinIceThickness = GetConstReal(Model % Constants, "MinIceThickness", Found) IF(.NOT.Found) CALL Fatal("getMinSurfaceElevation", "Unable to find MinIceThickness") BedrockDEM = GetString(Model % Constants, "BedrockDEM", Found) IF(.NOT.Found) CALL Fatal("getMinSurfaceElevation", "Unable to find BedrockDEM") Nbx = 801 Nby = 420 xb0 = 610010.31250 yb0 = 5108310.50000 lbx = 24000 lby = 12570 Nb = Nbx * Nby IF (FirstTime) THEN FirstTime = .False. OPEN(10,file=BedrockDEM) ALLOCATE(xb(Nb), yb(Nb), zb(Nb)) READ(10,*)(xb(i), yb(i), zb(i), i=1,Nb) CLOSE(10) END IF ! zbed for point (x,y) x = Model % Nodes % x (nodenumber) y = Model % Nodes % y (nodenumber) zbed = InterpolateDEM(x,y,xb,yb,zb,Nbx,Nby,xb0,yb0,lbx,lby) Zsurf = zbed + MinIceThickness END FUNCTION getMinSurfaceElevation !!!!!!!!!!!!!!!!!!! ! Returns the velocity of the bedrock to be set as boundary condition. !!------------------------------------------------------------------------------!! FUNCTION getBottomVelocity ( Model, nodenumber) RESULT(bottomVelocity) USE types USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver TYPE(Variable_t), POINTER :: TimeVar TYPE(Mesh_t), POINTER :: MeshVar INTEGER :: nodenumber, timeStatus, checkTime REAL(KIND=dp) :: BasalWinterVelocity, BasalSummerVelocity, SwitchPeriod REAL(KIND=dp) :: bottomVelocity, prevVelocityValue=0, Time LOGICAL :: Found SAVE prevVelocityValue SwitchPeriod = GetConstReal(Model % Constants, "SwitchPeriod", Found) IF(.NOT.Found) CALL Fatal("getBottomVelocity", "Unable to find SwitchPeriod") BasalWinterVelocity = GetConstReal(Model % Constants, "BasalWinterVelocity", Found) IF(.NOT.Found) CALL Fatal("getBottomVelocity", "Unable to find BasalWinterVelocity") BasalSummerVelocity = GetConstReal(Model % Constants, "BasalSummerVelocity", Found) IF(.NOT.Found) CALL Fatal("getBottomVelocity", "Unable to find BasalSummerVelocity") TimeVar => VariableGet( Model % Variables, 'Time' ) Time = TimeVar % Values(1) ! check whether to apply the accumulation timeStatus = checkTime(Time, SwitchPeriod); IF (timeStatus == 1) THEN bottomVelocity = BasalWinterVelocity; ! CALL INFO("getBottomVelocity", "Switch to 1st month of the year.", Level=3) ELSE IF (timeStatus == 0) THEN bottomVelocity = BasalSummerVelocity; ! CALL INFO("getBottomVelocity", "Switch to middle of the year.", Level=3) ELSE bottomVelocity = prevVelocityValue ! CALL INFO("getBottomVelocity", "Keep bottomVelocity of previous timestep", Level=3) END IF IF(prevVelocityValue /= bottomVelocity) THEN WRITE(Message, '(A, ES14.7)') 'Set value of bottomVelocity to ', bottomVelocity CALL INFO("getBottomVelocity", Message, Level=3) END IF prevVelocityValue = bottomVelocity; END FUNCTION getBottomVelocity !!!!!!!!!!!!!!!!!!! ! Returns the arrhenius factor. !!------------------------------------------------------------------------------!! FUNCTION getArrheniusFactor ( Model, nodenumber) RESULT(Arrhenius) USE types USE SolverUtils USE ElementDescription USE DefUtils IMPLICIT NONE TYPE(Model_t) :: Model TYPE(Solver_t), TARGET :: Solver TYPE(Variable_t), POINTER :: TimeVar TYPE(Mesh_t), POINTER :: MeshVar INTEGER :: nodenumber, timeStatus, checkTime REAL(KIND=dp) :: WinterArrheniusFactor, SummerArrheniusFactor REAL(KIND=dp) :: Arrhenius, prevArrheniusValue, Time, SwitchPeriod LOGICAL :: Found SAVE prevArrheniusValue TimeVar => VariableGet( Model % Variables, 'Time' ) Time = TimeVar % Values(1) SwitchPeriod = GetConstReal(Model % Constants, "SwitchPeriod", Found) IF(.NOT.Found) CALL Fatal("getBalance", "Unable to find SwitchPeriod") WinterArrheniusFactor = GetConstReal(Model % Constants, "WinterArrheniusFactor", Found) IF(.NOT.Found) CALL Fatal("getArrheniusFactor", "Unable to find WinterArrheniusFactor") SummerArrheniusFactor = GetConstReal(Model % Constants, "SummerArrheniusFactor", Found) IF(.NOT.Found) CALL Fatal("getArrheniusFactor", "Unable to find SummerArrheniusFactor") ! check whether to apply the accumulation timeStatus = checkTime(Time, SwitchPeriod); IF (timeStatus == 1) THEN Arrhenius = WinterArrheniusFactor; ! CALL INFO("getArrheniusFactor", "Switch to 1st month of the year.", Level=3) ELSE IF (timeStatus == 0) THEN Arrhenius = SummerArrheniusFactor; ! CALL INFO("getArrheniusFactor", "Switch to middle of the year.", Level=3) ELSE Arrhenius = prevArrheniusValue; ! CALL INFO("getArrheniusFactor", "Keep Arrhenius Factor of previous timestep", Level=3) END IF IF(prevArrheniusValue /= Arrhenius) THEN WRITE(Message, '(a, ES14.7)') 'Set value of Arrhenius Factor to ', Arrhenius CALL INFO("getArrheniusFactor", Message, Level=3) END IF prevArrheniusValue = Arrhenius; END FUNCTION getArrheniusFactor