CoCalc -- USF_IceProperties.F90

GitHub Repository: ElmerCSC/elmerfem
Path: blob/devel/elmerice/UserFunctions/USF_IceProperties.F90
³¹⁹⁶ views
1
!/*****************************************************************************/
2
! *
3
! *  Elmer/Ice, a glaciological add-on to Elmer
4
! *  http://elmerice.elmerfem.org
5
! *
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
! *  Authors: Denis Cohen, Thomas Zwinger
27
! *  Email:   
28
! *  Web:     http://elmerice.elmerfem.org
29
! *
30
! Contains functions for material properties of ice:
31
! - IceConductivity_SI
32
! - IceConductivity_m_Mpa_a
33
! - IceCapacity_SI
34
! - IceCapacity_m_MPa_a
35
! - IcePressureMeltingPoint
36
!
37
MODULE IceProperties
38
  USE DefUtils
39
  USE Types
40
  
41
  IMPLICIT None
42

43
  CONTAINS
44

45
  !==============================================================================
46
  FUNCTION GetIceConductivity(temp) RESULT(cond)
47
  !==============================================================================
48

49

50
    IMPLICIT None
51
    REAL(KIND=dp) :: temp, cond
52

53
    
54
    
55
    cond = 9.828_dp*exp(-5.7d-03*temp)
56

57
  END FUNCTION GetIceConductivity
58
  
59
  !==============================================================================
60
  FUNCTION GetIceCapacity(temp) RESULT(capac)
61
  !==============================================================================
62

63
    USE DefUtils
64

65
    IMPLICIT None
66

67
    REAL(KIND=dp) :: temp, capac
68

69

70

71
    capac = 146.3_dp + (7.253_dp * temp)
72
  END FUNCTION GetIceCapacity
73
  
74
  !==============================================================================
75
  FUNCTION GetIcePressureMeltingPoint(ClausiusClapeyron, press) RESULT(Tpmp)
76
  !==============================================================================
77

78
    USE DefUtils
79
    
80
    IMPLICIT None
81

82
    REAL(KIND=dp) :: Tpmp, ClausiusClapeyron, press
83

84
    Tpmp = 273.15_dp - ClausiusClapeyron*MAX(press, 0.0_dp)
85

86
  END FUNCTION GetIcePressureMeltingPoint
87

88
END MODULE IceProperties
89

90
!==============================================================================
91
FUNCTION IceConductivity(Model, Node, temp) RESULT(cond)
92
!==============================================================================
93
  USE IceProperties
94
  
95
  TYPE(Model_t) :: Model
96
  INTEGER :: Node
97
  REAL(KIND=dp) :: temp, cond
98

99
  REAL(KIND=dp) :: scalingfactor
100
  TYPE(Element_t), POINTER :: Element
101
  TYPE(ValueList_t), POINTER  :: Material
102
  LOGICAL :: Found
103
  
104
  Element => Model % CurrentElement
105
  Material => GetMaterial(Element)
106
  scalingfactor = GetConstReal(Material,"Heat Conductivity Scaling Factor",Found)
107
  IF (.NOT.Found) scalingfactor = 1.0_dp
108
  cond = scalingfactor * GetIceConductivity(temp)
109
END FUNCTION IceConductivity
110

111
!==============================================================================
112
FUNCTION IceCapacity(Model, Node, temp) RESULT(capac)
113
!==============================================================================
114
  USE IceProperties
115

116
  IMPLICIT None
117

118
  TYPE(Model_t) :: Model
119
  INTEGER :: Node
120
  REAL(KIND=dp) :: temp, capac
121

122
  REAL(KIND=dp) :: scalingfactor
123
  TYPE(Element_t), POINTER :: Element
124
  TYPE(ValueList_t), POINTER  :: Material
125
  LOGICAL :: Found
126
  
127
  Element => Model % CurrentElement
128
  Material => GetMaterial(Element)
129
  scalingfactor = GetConstReal(Material,"Heat Capacity Scaling Factor",Found)
130
  IF (.NOT.Found) scalingfactor = 1.0_dp
131
  capac = scalingfactor * GetIceCapacity(temp)
132
END FUNCTION IceCapacity
133

134
!==============================================================================
135
FUNCTION IcePressureMeltingPoint(Model, Node, press) RESULT(Tpmp)
136
!==============================================================================
137

138
  USE IceProperties
139

140
  IMPLICIT None
141

142
  TYPE(Model_t) :: Model
143
  INTEGER :: Node
144
  REAL(KIND=dp) :: Tpmp, press
145

146
  INTEGER :: N
147
  REAL(KIND=dp) :: ClausiusClapeyron, tmpoffset
148
  TYPE(ValueList_t), POINTER :: Constants
149
  LOGICAL :: FirstTime = .TRUE., GotIt, InCelsius
150
  REAL(KIND=dp) :: scalingfactor=1.0_dp
151
  TYPE(Element_t), POINTER :: Element
152
  TYPE(ValueList_t), POINTER  :: Material, BC
153
  
154
  SAVE FirstTime, ClausiusClapeyron
155

156
  Element => Model % CurrentElement
157
  Material => GetMaterial(Element)
158
  IF (.NOT.ASSOCIATED(Material)) THEN
159
    BC => GetBC(Element)
160
    IF (.NOT.ASSOCIATED(BC)) THEN
161
      scalingfactor=1.0_dp
162
      CALL INFO("IcePressureMeltingPoint", 'No "Pressure Scaling Factor" found - setting to 1.0',Level=9)
163
    ELSE
164
      scalingfactor = GetConstReal(BC,"Pressure Scaling Factor",GotIt)
165
    END IF
166
  ELSE
167
    scalingfactor = GetConstReal(Material,"Pressure Scaling Factor",GotIt)
168
  END IF
169
  IF (.NOT.GotIt) scalingfactor = 1.0_dp
170
  InCelsius = GetLogical(Material, "Pressure Melting Point Celsius", GotIt)
171
  IF (InCelsius) THEN
172
    tmpoffset = 273.15_dp
173
  ELSE
174
    tmpoffset = 0.0_dp
175
  END IF
176
  IF (FirstTime) THEN
177
    FirstTime = .FALSE.
178
    Constants => GetConstants()
179
    IF (ASSOCIATED(Constants)) THEN
180
      ClausiusClapeyron = GetConstReal( Constants, 'Clausius Clapeyron Constant', GotIt)
181
    ELSE
182
      GotIt=.FALSE.
183
    ENDIF
184
    IF (.NOT.GotIt) THEN
185
      ClausiusClapeyron = 9.8d-08
186
      CALL INFO("IcePressureMeltingPoint","No entry found for >Clausius Clapeyron Constant<.",Level=9)
187
      CALL INFO("IcePressureMeltingPoint","Setting to 9.8d-08 (SI units)",Level=9)
188
    END IF
189
  END IF
190
  Tpmp = GetIcePressureMeltingPoint(ClausiusClapeyron,scalingfactor*press) + tmpoffset
191
END FUNCTION IcePressureMeltingPoint
192
!==============================================================================
193
FUNCTION RelativeTemperature(Model, Node, InputArray) RESULT(Trel)
194
!==============================================================================
195
  USE IceProperties
196
  USE DefUtils
197
  IMPLICIT None
198

199
  TYPE(Model_t) :: Model
200
  INTEGER :: Node
201
  REAL(KIND=dp) :: Trel, InputArray(2)
202
  ! ----
203
  REAL(KIND=dp) :: press, Tabs, Tpmp,  ClausiusClapeyron, tmpoffset, scalingfactor=1.0_dp
204
  TYPE(Element_t), POINTER :: Element
205
  TYPE(ValueList_t), POINTER  :: Material, BC, Constants
206
  LOGICAL :: CorrectRelTemp = .FALSE., GotIt, FirstTime = .TRUE., InCelsius = .FALSE.
207

208
  SAVE FirstTime, ClausiusClapeyron
209
  
210
  Tabs = InputArray(1)
211
  press = InputArray(2)
212
  Element => Model % CurrentElement
213
  Material => GetMaterial(Element)
214
  IF (.NOT.ASSOCIATED(Material)) THEN
215
    BC => GetBC(Element)
216
    IF (.NOT.ASSOCIATED(BC)) THEN
217
      scalingfactor=1.0_dp
218
      CALL INFO("IcePressureMeltingPoint", 'No "Pressure Scaling Factor" found - setting to 1.0',Level=9)
219
    ELSE
220
      scalingfactor = GetConstReal(BC,"Pressure Scaling Factor",GotIt)
221
    END IF
222
  ELSE
223
    scalingfactor = GetConstReal(Material,"Pressure Scaling Factor",GotIt)
224
  END IF
225
  IF (.NOT.GotIt) scalingfactor = 1.0_dp
226
  IF (FirstTime) THEN
227
    FirstTime = .FALSE.
228
    Constants => GetConstants()
229
    CorrectRelTemp=ListGetLogical(Material,"Correct Relative Tempeature", GotIt)
230
    IF (ASSOCIATED(Constants)) THEN
231
      ClausiusClapeyron = GetConstReal( Constants, 'Clausius Clapeyron Constant', GotIt)
232
    ELSE
233
      GotIt=.FALSE.
234
    ENDIF
235
    IF (.NOT.GotIt) THEN
236
      ClausiusClapeyron = 9.8d-08
237
      CALL INFO("IcePressureMeltingPoint","No entry found for >Clausius Clapeyron Constant<.",Level=9)
238
      CALL INFO("IcePressureMeltingPoint","Setting to 9.8d-08 (SI units)",Level=9)
239
    END IF
240
  END IF
241

242
  InCelsius = GetLogical(Material, "Pressure Melting Point Celsius", GotIt)
243
  IF (InCelsius) THEN
244
    tmpoffset = 273.15_dp
245
  ELSE
246
    tmpoffset = 0.0_dp
247
  END IF
248
  Tpmp = GetIcePressureMeltingPoint(ClausiusClapeyron,scalingfactor*press) + tmpoffset  
249
  Trel = Tabs - Tpmp
250
  Element => Model % CurrentElement
251
  Material => GetMaterial(Element)
252
  IF (ASSOCIATED(Material)) THEN
253
    CorrectRelTemp=ListGetLogical(Material,"Correct Relative Tempeature", GotIt)
254
    IF (.NOT.GotIt) THEN
255
      CorrectRelTemp=.FALSE.
256
    ELSE IF (CorrectRelTemp) THEN
257
      IF (FirstTime) THEN
258
        CALL INFO("USF_IceProperties(IcePressureMeltingPoint)","Limiting Relative Temperature",Level=3)
259
        FirstTime = .FALSE.
260
      ENDIF
261
      Trel = MIN(Trel, 0.0_dp)
262
    ENDIF
263
  ENDIF
264
END FUNCTION RelativeTemperature
265
!==============================================================================
266
FUNCTION ArrheniusFactor(Model, Node, InputArray) RESULT(ArrhF)
267
!==============================================================================
268

269
  USE IceProperties
270

271
  IMPLICIT None
272

273
  TYPE(Model_t) :: Model
274
  INTEGER :: Node
275
  REAL(KIND=dp) :: InputArray(2), Arrhf
276

277
  INTEGER :: N
278
  REAL(KIND=dp) ::  Temp, Pressure
279
  REAL(KIND=dp) :: Tempoffset, GlenExponent, GlenEnhancementFactor, &
280
       ClausiusClapeyron, GasConstant, LimitTemperature,&
281
        Tpmp, Trel, scalingfactor, RateFactor(2), ActivationEnergy(2)
282
  TYPE(ValueList_t), POINTER :: Constants
283
  LOGICAL ::  GotIt, InCelsius, ScaleRateFactor
284
  TYPE(Element_t), POINTER :: Element
285
  TYPE(ValueList_t), POINTER  :: Material
286
  
287

288

289
  
290
  Element => Model % CurrentElement
291
  Material => GetMaterial(Element)
292
  Constants => GetConstants()
293

294
  GasConstant = GetConstReal(Constants, "Gas Constant",GotIt)
295
  IF (.NOT.GotIt) GasConstant = 8.314_dp
296
  ClausiusClapeyron = GetConstReal( Constants, 'Clausius Clapeyron Constant', GotIt)
297
  IF (.NOT.GotIt) THEN
298
    ClausiusClapeyron = 9.8d-08
299
    CALL INFO("IcePressureMeltingPoint","No entry found for >Clausius Clapeyron Constant<.",Level=5)
300
    CALL INFO("IcePressureMeltingPoint","Setting to 9.8d-08 (SI units)",Level=5)
301
  END IF
302
  scalingfactor = GetConstReal(Material,"Pressure Scaling Factor",GotIt)
303
  IF (.NOT.GotIt) scalingfactor = 1.0_dp
304
  InCelsius = GetLogical(Material, "Temperature in Celsius", GotIt)
305
  IF (InCelsius) THEN
306
    Tempoffset = 0.0_dp
307
  ELSE
308
    Tempoffset = 273.15_dp
309
  END IF
310

311
  Temp = InputArray(1) - Tempoffset
312
  Pressure = InputArray(2)
313
  
314
  GlenExponent = GetConstReal(Material,"Glen Exponent",GotIt)
315
   IF (.NOT.GotIt) THEN
316
     CALL INFO("IceProperties(ArrheniusFactor)",&
317
          '"Glen Exponent" not found. Setting to 3.0',Level=5)
318
    GlenExponent = 3.0_dp
319
  END IF  
320
  LimitTemperature = GetConstReal(Material,"Limit Temperature",GotIt)
321
  IF (.NOT.GotIt) THEN
322
    CALL INFO("IceProperties(ArrheniusFactor)",&
323
         '"Limit Temperature" not found. Setting to -10.0',Level=5)
324
    LimitTemperature = -10.0_dp
325
  END IF
326
  RateFactor(1) = GetConstReal(Material,"Rate Factor 1", GotIt)
327
  IF (.NOT.GotIt) THEN
328
    CALL INFO("IceProperties(ArrheniusFactor)",&
329
         '"Rate Factor 1" not found. Setting to (SI value) 2.89165e-13',Level=5)
330
    RateFactor(1) = 2.89165d-13
331
  END IF
332
  RateFactor(2) = GetConstReal(Material,"Rate Factor 2", GotIt)
333
  IF (.NOT.GotIt) THEN
334
    CALL INFO("IceProperties(ArrheniusFactor)",&
335
         '"Rate Factor 2" not found. Setting to (SI value) 2.42736e-02',Level=5)
336
    RateFactor(2) = 2.42736d-02
337
  END IF
338
  ScaleRateFactor = GetLogical(Material,"Scale Rate Factors",GotIt)
339
  IF (ScaleRateFactor) RateFactor = RateFactor*(scalingfactor**GlenExponent)
340
  
341
  ActivationEnergy(1) = GetConstReal(Material,"Activation Energy 1", GotIt)
342
  IF (.NOT.GotIt) THEN
343
    CALL INFO("IceProperties(ArrheniusEnergy)",&
344
         '"Activation Energy 1" not found. Setting to (SI value) 60.0e03',Level=5)
345
    ActivationEnergy(1) = 60.0d03
346
  END IF
347
  ActivationEnergy(2) = GetConstReal(Material,"Activation Energy 2", GotIt)
348
  IF (.NOT.GotIt) THEN
349
    CALL INFO("IceProperties(ArrheniusEnergy)",&
350
         '"Activation Energy 2" not found. Setting to (SI value) 115.0e3',Level=5)
351
    ActivationEnergy(2) = 115.d03
352
  END IF
353
 
354
  GlenEnhancementFactor = GetConstReal(Material,"Glen Enhancement Factor",GotIt)
355
  IF (.NOT.GotIt) THEN
356
    CALL INFO("IceProperties(ArrheniusEnergy)",&
357
         '"Glen Enhancement Factor" not found. Setting to 1.0', Level=5)
358
    GlenEnhancementFactor = 1.0_dp
359
  END IF
360

361
  ! need to scale pressure 
362
  Tpmp = GetIcePressureMeltingPoint(ClausiusClapeyron, scalingfactor * Pressure)
363
   
364
  Trel = MAX(MIN(Temp - Tpmp,0.0_dp),-60.0_dp)
365
  IF (Trel<-10) THEN
366
    Arrhf = RateFactor(1) * EXP( -ActivationEnergy(1)/( GasConstant* (273.15 + Trel)))
367
  ELSE
368
    Arrhf= RateFactor(2) * EXP( -ActivationEnergy(2)/( GasConstant* (273.15 + Trel)))
369
  END IF
370
END FUNCTION ArrheniusFactor
371

372
Product

Resources

Company
