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) asny 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! * Authors: Thomas Zwinger, Denis Cohen, Juha Hartikainen 26! * Email: thomas Zwinger [at] csc.fi 27! * Web: http://elmerice.elmerfem.org 28! * Address: CSC - Scientific Computing Ltd. 29! * Keilaranta 14 30! * 02101 Espoo, Finland 31! * 32! * Original Date: January 2017 - 33! * 34! ***************************************************************************** 35!> Module containing material declarations and material functions 36!> for enhanced permafrost problem 37! ***************************************************************************** 38MODULE PermafrostMaterials 39 40 USE Types 41 USE DefUtils 42 USE SolverUtils 43 IMPLICIT NONE 44 !--------------------------------- 45 ! type for solvent (water and ice) 46 !--------------------------------- 47 TYPE SolventMaterial_t 48 REAL(KIND=dp) :: & 49 Mw,rhow0,rhoi0,hw0,hi0,vi0,bccw0,& 50 Ei0, nui0, betai, & 51 kw0th,ki0th,bw,bi, & 52 cw0,acw(0:5),bcw(0:5), & 53 ci0,aci(0:5),& 54 aw0,kw0,zw0,aaw(0:5),bzw(0:5),ckw(0:5), & 55 ai0,ki0,aai(0:5),cki(0:5),& 56 muw0,nu10,anw(0:5),bnw(0:5) 57 INTEGER :: & 58 acwl,bcwl,aawl,bzwl,ckwl,& 59 acil,aail,ckil,anwl,bnwl 60 END type SolventMaterial_t 61 62 !--------------------------------- 63 ! type for solute (ions) 64 !--------------------------------- 65 TYPE SoluteMaterial_t 66 REAL(KIND=dp) :: Mc,vc0,kc0th,& 67 d1,d2,bc,& 68 cc0,acc(0:5),bcc(0:5),& 69 rhoc0,ac0,kc0,zc0,aac(0:5),ckc(0:5),bzc(0:5),& 70 nu20,anc(0:5),bnc(0:5) 71 INTEGER :: accl, bccl,aacl,ckcl,bzcl,ancl,bncl 72 CHARACTER(LEN=MAX_NAME_LEN) :: SoluteName 73 END type SoluteMaterial_t 74 75 !--------------------------------- 76 ! type for rock material 77 !--------------------------------- 78 TYPE RockMaterial_t 79 INTEGER :: NumerOfRockRecords 80 REAL(KIND=dp), ALLOCATABLE :: ks0th(:),e1(:),bs(:),rhos0(:),& 81 Xi0(:),eta0(:),etak(:),hs0(:),Kgwh0(:,:,:),qexp(:),alphaL(:),alphaT(:),RadGen(:),& 82 cs0(:),acs(:,:),as0(:),aas(:,:),ks0(:),cks(:,:),Es0(:),nus0(:),betas(:) 83 INTEGER, ALLOCATABLE :: acsl(:),aasl(:),cksl(:) 84 CHARACTER(LEN=MAX_NAME_LEN), ALLOCATABLE :: VariableBaseName(:) 85 END TYPE RockMaterial_t 86 87 TYPE(SolventMaterial_t), TARGET :: GlobalSolventMaterial 88 TYPE(SoluteMaterial_t), TARGET :: GlobalSoluteMaterial 89 TYPE(RockMaterial_t) :: GlobalRockmaterial 90 91 92CONTAINS 93 !--------------------------------------------------------------------------------------------- 94 !--------------------------------------------------------------------------------------------- 95 !--------------------------------------------------------------------------------------------- 96 97 !------------------------------------------------- 98 ! I/O related functions 99 !------------------------------------------------- 100 101 SUBROUTINE SetPermafrostSolventMaterial( CurrentSolventMaterial) 102 IMPLICIT NONE 103 TYPE(ValueList_t), POINTER :: Params, Constants 104 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 105 ! ----------- local 106 TYPE(SolventMaterial_t), TARGET :: LocalSolventMaterial 107 LOGICAL :: FirstTime=.TRUE. 108 CHARACTER(LEN=MAX_NAME_LEN) :: SubroutineName='SetPermafrostSolventMaterial' 109 SAVE LocalSolventMaterial, FirstTime 110 111 IF (FirstTime) THEN 112 !------------------------------------------------------------------------------ 113 ! set constants for water and ice 114 ! Mw,rhow0,rhoi0,hw0,hi0,vi0,cw0,ci0,acw(3),bcw(0:2),aci(0:1),kw0th,ki0th, bi, bw 115 !------------------------------------------------------------------------------ 116 LocalSolventMaterial % Mw = 1.8015d-2 117 LocalSolventMaterial % hw0 = 0.0_dp 118 LocalSolventMaterial % hi0 = -333360.0_dp !! 119 120 ! --------------------- polynomials 121 122 !!! water !!! 123 124 ! heat capacity water 125 LocalSolventMaterial % cw0 = 4207.7_dp 126 LocalSolventMaterial % acw(0:5) = & 127 RESHAPE([1.0_dp,-0.0887_dp,0.2859_dp,0.0_dp,0.0_dp,0.0_dp], & 128 SHAPE(LocalSolventMaterial % acw)) 129 LocalSolventMaterial % acwl=2 130 LocalSolventMaterial % bcw(0:5) = & 131 RESHAPE([1.0_dp,1.5852_dp,8.0686_dp,0.0_dp,0.0_dp,0.0_dp],& 132 SHAPE(LocalSolventMaterial % bcw)) 133 LocalSolventMaterial % bcwl=2 134 135 !heat conductivity of water 136 LocalSolventMaterial % kw0th = 0.56_dp 137 LocalSolventMaterial % bw = 0.0_dp 138 139 ! density water 140 LocalSolventMaterial % rhow0 = 999.9_dp ! density at reference temperature 141 LocalSolventMaterial % kw0 = 4.4534d-10 ! Isothermal compressibility 142 143 LocalSolventMaterial % ckw(0:5) = & 144 RESHAPE([1.0_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp],& 145 SHAPE(LocalSolventMaterial % ckw)) 146 LocalSolventMaterial % ckwl=0 147 LocalSolventMaterial % aw0 = -5.3358d-05 ! Isobaric thermal expansion 148 LocalSolventMaterial % aaw(0:5) = & 149 RESHAPE([1.0_dp,-79.1305_dp,207.4836_dp,-403.8270_dp,395.5347_dp,-166.1466_dp],& 150 SHAPE(LocalSolventMaterial % aaw)) 151 LocalSolventMaterial % aawl=5 152 LocalSolventMaterial % zw0 = -2.0217d-01 ! Isothermal chemical compaction 153 LocalSolventMaterial % bzw(0:5) = & 154 RESHAPE([1.0_dp,12.8298_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp],& 155 SHAPE(LocalSolventMaterial % bzw)) 156 LocalSolventMaterial % bzwl=1 157 158 !viscosity water 159 LocalSolventMaterial % muw0 = 1.7914d-03 ! viscosity at reference temperature 160 LocalSolventMaterial % nu10 = -0.034376_dp ! temperature dependence of viscosity 161 LocalSolventMaterial % anw(0:5) = & 162 RESHAPE([1.0_dp,-2.3302_dp,4.0084_dp,-2.9697_dp,0.0_dp,0.0_dp],& 163 SHAPE(LocalSolventMaterial % anw)) 164 LocalSolventMaterial % anwl=3 165 LocalSolventMaterial % bnw(0:5) = & 166 RESHAPE([1.0_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp],& 167 SHAPE(LocalSolventMaterial % bnw)) 168 LocalSolventMaterial % bnwl=0 169 170 !!!!! ice !!!!!! 171 172 ! density ice 173 LocalSolventMaterial % rhoi0 = 916.8_dp !reference density 174 LocalSolventMaterial % ki0 = 1.1417d-10 ! Isothermal compressibility 175 LocalSolventMaterial % cki(0:5) = & 176 RESHAPE([1.0_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp],& 177 SHAPE(LocalSolventMaterial % cki)) 178 LocalSolventMaterial % ckil=0 179 LocalSolventMaterial % ai0 = 1.6781d-04 ! Isobaric thermal expansion 180 LocalSolventMaterial % aai(0:5) = & 181 RESHAPE([1.0_dp,1.1923_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp], & 182 SHAPE(LocalSolventMaterial % aai)) 183 LocalSolventMaterial % aail=1 184 185 ! specivic volume ice 186 LocalSolventMaterial % vi0 = 1.0_dp/(LocalSolventMaterial % rhoi0) ! reference specific volume 187 188 ! heat conductivity ice 189 LocalSolventMaterial % ki0th = 2.24_dp!! 190 LocalSolventMaterial % bi = 0.0_dp 191 192 ! heat capacity 193 LocalSolventMaterial % ci0 = 2088.8_dp ! reference value 194 LocalSolventMaterial % aci(0:5) = & 195 RESHAPE([1.0_dp,0.9557_dp,0.0_dp,0.0_dp,0.0_dp,0.0_dp], & 196 SHAPE(LocalSolventMaterial % aci)) 197 LocalSolventMaterial % acil=1 198 199 ! elastic properties 200 LocalSolventMaterial % Ei0 = 9.33d09 201 LocalSolventMaterial % nui0= 0.325_dp 202 LocalSolventMaterial % betai = 0.0_dp ! CHANGE 203 204 CALL INFO(SubroutineName,"-----------------------------------------------------------------",Level=9) 205 CALL INFO(SubroutineName,"Solvent related constants",Level=9) 206 WRITE(Message,*) "Mw",LocalSolventMaterial % Mw,"rhow0",LocalSolventMaterial % rhow0,"rhoi0",LocalSolventMaterial % rhoi0,& 207 "hw0",LocalSolventMaterial % hw0,"hi0",LocalSolventMaterial % hi0,"vi0",LocalSolventMaterial % vi0,& 208 "cw0",LocalSolventMaterial % cw0,"ci0",LocalSolventMaterial % ci0,"acw(3)",LocalSolventMaterial % acw(0:2) 209 CALL INFO(SubroutineName,Message,Level=9) 210 WRITE(Message,*) "bcw(0:2)",LocalSolventMaterial % bcw(0:2),"aci(0:1)",LocalSolventMaterial % aci(0:1),& 211 "kw0th",LocalSolventMaterial % kw0th,"ki0th",LocalSolventMaterial % ki0th," bi",LocalSolventMaterial % bi,& 212 "bw",LocalSolventMaterial % bw 213 CALL INFO(SubroutineName,Message,Level=9) 214 CALL INFO(SubroutineName,"-----------------------------------------------------------------",Level=9) 215 FirstTime = .FALSE. 216 END IF 217 CurrentSolventMaterial => LocalSolventMaterial 218 END SUBROUTINE SetPermafrostSolventMaterial 219 220 !--------------------------------------------------------------------------------------------- 221 SUBROUTINE ReadPermafrostSoluteMaterial( Params,Constants,CurrentSoluteMaterial ) 222 IMPLICIT NONE 223 TYPE(ValueList_t), POINTER :: Params, Constants 224 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 225 ! ----------- local 226 TYPE(SoluteMaterial_t), TARGET :: LocalSoluteMaterial 227 INTEGER :: i,j,k,l, n,t, active, DIM, ok,InitialNumerOfSoluteRecords, EntryNumber 228 INTEGER,parameter :: io=20 229 LOGICAL :: Found, DataRead=.FALSE. 230 CHARACTER(LEN=MAX_NAME_LEN) :: SoluteFileName, Comment 231 CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: SubroutineName='ReadPermafrostSoluteMaterial' 232 233 SAVE DataRead,SoluteFileName, LocalSoluteMaterial 234 235 IF (DataRead) THEN 236 CurrentSoluteMaterial => LocalSoluteMaterial 237 RETURN 238 ELSE 239 DIM = CoordinateSystemDimension() 240 !------------------------------------------------------------------------------ 241 ! Inquire and open file 242 !------------------------------------------------------------------------------ 243 ! give preference to a defined material database 244 SoluteFileName = GetString( Params, 'Solute Material File', Found ) 245 IF (.NOT.Found) THEN 246 CALL INFO(SubroutineName," 'Solute Material File' keyword not found - assigning default values!",Level=1) 247 DataRead=.TRUE. 248 LocalSoluteMaterial % SoluteName = TRIM('Sea Salt') 249 LocalSoluteMaterial % Mc = 0.031404_dp 250 LocalSoluteMaterial % rhoc0 = 2206.6_dp 251 LocalSoluteMaterial % vc0 = 1.0/(LocalSoluteMaterial % rhoc0) 252 LocalSoluteMaterial % kc0th = 0.56_dp 253 LocalSoluteMaterial % d1 = 0.87_dp 254 LocalSoluteMaterial % d2 = 2.00_dp 255 LocalSoluteMaterial % bc = 0.0_dp 256 ! heat capacity polynomials 257 LocalSoluteMaterial % cc0 = 1906.6_dp 258 LocalSoluteMaterial % acc(0:5) = & 259 RESHAPE([1.0_dp,-0.0887_dp,0.2859_dp,0.0_dp,0.0_dp,0.0_dp],& 260 SHAPE(LocalSoluteMaterial % acc)) 261 LocalSoluteMaterial % accl = 2 262 LocalSoluteMaterial % bcc(0:5) = & 263 RESHAPE([1.0_dp,-1.5852_dp,8.0686_dp,0.0_dp,0.0_dp,0.0_dp],& 264 SHAPE(LocalSoluteMaterial % bcc)) 265 LocalSoluteMaterial % bccl = 2 266 ! density 267 LocalSoluteMaterial % ac0 = -5.3358d-05 ! thermal expansion 268 LocalSoluteMaterial % kc0 = 4.4534d-10 ! compressibility 269 LocalSoluteMaterial % zc0 = 2.0217d-01 ! chemical compaction 270 LocalSoluteMaterial % aac = & 271 RESHAPE([1.0_dp, -79.1305_dp, 207.4835_dp, -403.827_dp, 395.5347_dp, -166.1466_dp],& 272 SHAPE(LocalSoluteMaterial % aac)) 273 LocalSoluteMaterial % aacl = 5 274 LocalSoluteMaterial % ckc = & 275 RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp],& 276 SHAPE(LocalSoluteMaterial % ckc)) 277 LocalSoluteMaterial % ckcl = 0 278 LocalSoluteMaterial % bzc = & 279 RESHAPE([1.0_dp, -12.8298_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp],& 280 SHAPE(LocalSoluteMaterial % bzc)) 281 LocalSoluteMaterial % bzcl = 1 282 ! viscosity 283 LocalSoluteMaterial % nu20 = 2.6870_dp !influence of salinity on viscosity 284 LocalSoluteMaterial % anc = & 285 RESHAPE([1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp],& 286 SHAPE(LocalSoluteMaterial % anc)) 287 LocalSoluteMaterial % ancl = 0 288 LocalSoluteMaterial % bnc= & 289 RESHAPE([1.0_dp, 6.7992_dp, -31.3293_dp, 44.7717_dp, 0.0_dp, 0.0_dp],& 290 SHAPE(LocalSoluteMaterial % bnc)) 291 LocalSoluteMaterial % bncl = 3 292 DataRead=.TRUE. 293 ELSE 294 OPEN(unit = io, file = TRIM(SoluteFileName), status = 'old',iostat = ok) 295 IF (ok /= 0) THEN 296 WRITE(Message,'(A,A)') 'Unable to open file ',TRIM(SoluteFileName) 297 CALL FATAL(Trim(SubroutineName),Trim(message)) 298 ELSE 299 !------------------------------------------------------------------------------ 300 ! Read in the number of records in file (first line integer) 301 !------------------------------------------------------------------------------ 302 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % SoluteName 303 WRITE(Message,'(A,A)') 'Reading entry', TRIM(LocalSoluteMaterial % SoluteName) 304 CALL INFO(Trim(SubroutineName),Trim(Message),Level=3) 305 !------------------------------------------------------------------------------ 306 ! Read in information for (currently fixed) solute (= salts) 307 ! 308 ! Mc, rhoc0 309 ! vc0,cc0,acc(0:2),bcc(0:2) 310 ! kc0th,muw0,d1,d2,bc 311 !------------------------------------------------------------------------------ 312 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % Mc, Comment 313 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % vc0, Comment 314 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % kc0th, Comment 315 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % d1, Comment 316 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % d2, Comment 317 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bc, Comment 318 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % cc0, Comment 319 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % acc(0:5), Comment 320 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % accl, Comment 321 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bcc(0:5), Comment 322 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bccl, Comment 323 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % rhoc0, Comment 324 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % ac0, Comment 325 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % kc0, Comment 326 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % zc0 , Comment 327 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % aac(0:5), Comment 328 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % aacl, Comment 329 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % ckc(0:5), Comment 330 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % ckcl, Comment 331 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bzc(0:5), Comment 332 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bzcl , Comment 333 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % nu20, Comment 334 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % anc(0:5), Comment 335 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % ancl , Comment 336 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bnc(0:5), Comment 337 READ (io, *, END=10, IOSTAT=OK, ERR=20) LocalSoluteMaterial % bncl, Comment 338 DataRead=.TRUE. 33910 CLOSE(io) 340 IF (.NOT.DataRead) THEN 341 WRITE(Message,'(A,A,A)') 'Not all entries in "Solute material File" ',TRIM(SoluteFileName),' found.' 342 CALL FATAL(Trim(SubroutineName),Trim(message)) 343 END IF 344 END IF 345 END IF 346 CurrentSoluteMaterial => LocalSoluteMaterial 347 END IF 348 CALL INFO(SubroutineName,"-----------------------------------------------------------------",Level=9) 349 CALL INFO(SubroutineName,"Solute related constants",Level=9) 350 WRITE(Message,*) "Mc",CurrentSoluteMaterial % Mc,"vc0",CurrentSoluteMaterial % vc0,& 351 "kc0th", CurrentSoluteMaterial %kc0th 352 CALL INFO(SubroutineName,Message,Level=9) 353 WRITE(Message,*) "d1",CurrentSoluteMaterial % d1,"d2",& 354 CurrentSoluteMaterial % d2,"bc",CurrentSoluteMaterial % bc 355 CALL INFO(SubroutineName,Message,Level=9) 356 WRITE(Message,*) "cc0",CurrentSoluteMaterial % cc0,"acc(0:5)",CurrentSoluteMaterial % acc(0:5),& 357 "bcc(0:5)",CurrentSoluteMaterial % bcc(0:5) 358 CALL INFO(SubroutineName,Message,Level=9) 359 WRITE(Message,*) "rhoc0",CurrentSoluteMaterial % rhoc0,"ac0",CurrentSoluteMaterial % ac0,"kc0",& 360 CurrentSoluteMaterial % kc0,"zc0",CurrentSoluteMaterial % zc0 361 WRITE(Message,*) "aac(0:5)",CurrentSoluteMaterial % aac(0:5),"ckc(0:5)",CurrentSoluteMaterial % ckc(0:5),& 362 "bzc(0:5)",CurrentSoluteMaterial % bzc(0:5) 363 CALL INFO(SubroutineName,Message,Level=9) 364 WRITE(Message,*) "bzc(0:5)",CurrentSoluteMaterial % bzc(0:5),"bnc(0:5)",CurrentSoluteMaterial % bnc(0:5),& 365 "nu20",CurrentSoluteMaterial % nu20 366 CALL INFO(SubroutineName,Message,Level=9) 367 WRITE(Message,*) "aacl",LocalSoluteMaterial % aacl,"ckcl", LocalSoluteMaterial % ckcl 368 CALL INFO(SubroutineName,Message,Level=9) 369 WRITE(Message,*) "bzcl",LocalSoluteMaterial % bzcl,"ancl", LocalSoluteMaterial % ancl,& 370 "bncl",LocalSoluteMaterial % bncl 371 CALL INFO(SubroutineName,Message,Level=9) 372 CALL INFO(SubroutineName,"-----------------------------------------------------------------",Level=9) 373 RETURN 37420 WRITE(Message,'(A,A,A)') 'Not all entries in "Solute material File" ',TRIM(SoluteFileName),' found.' 375 CLOSE(io) 376 CALL FATAL(Trim(SubroutineName),Trim(message)) 377 END SUBROUTINE ReadPermafrostSoluteMaterial 378 379 !--------------------------------------------------------------------------------------------- 380 FUNCTION ReadPermafrostRockMaterial( Params ) RESULT(NumerOfRockRecords) 381 IMPLICIT NONE 382 TYPE(ValueList_t), POINTER :: Params 383 Integer :: NumerOfRockRecords 384 !-------------- 385 INTEGER :: i,j,k,l, n,t, active, DIM, ok,InitialNumerOfRockRecords, EntryNumber 386 INTEGER,parameter :: io=21 387 LOGICAL :: Found, fexist, FirstTime=.TRUE., AllocationsDone=.FALSE., DataRead=.FALSE. 388 CHARACTER(LEN=MAX_NAME_LEN) :: MaterialFileName, NewMaterialFileName, str, Comment 389 CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: FunctionName='ReadPermafrostRockMaterial' 390 391 SAVE AllocationsDone,DataRead,InitialNumerOfRockRecords,MaterialFileName 392 393 IF (DataRead) THEN 394 NumerOfRockRecords = InitialNumerOfRockRecords 395 RETURN 396 ELSE ! we read Data from file database 397 DIM = CoordinateSystemDimension() 398 !------------------------------------------------------------------------------ 399 ! Inquire and open file 400 !------------------------------------------------------------------------------ 401 ! give preference to a defined material database 402 MaterialFileName = GetString( Params, 'Rock Material File', Found ) 403 IF (.NOT.Found) THEN 404 CALL INFO(FunctionName," 'Rock Material File' keyword not found - looking for default DB!") 405 fexist = .FALSE. 406#ifdef USE_ISO_C_BINDINGS 407 str = 'ELMER_LIB' 408#else 409 str = 'ELMER_LIB'//CHAR(0) 410#endif 411 CALL envir( str,MaterialFileName,k ) 412 IF ( k > 0 ) THEN 413 MaterialFileName = MaterialFileName(1:k) // '/permafrostmaterialdb.dat' 414 INQUIRE(FILE=TRIM(MaterialFileName), EXIST=fexist) 415 END IF 416 IF (.NOT. fexist) THEN 417#ifdef USE_ISO_C_BINDINGS 418 str = 'ELMER_HOME' 419#else 420 str = 'ELMER_HOME'//CHAR(0) 421#endif 422 CALL envir( str,MaterialFileName,k ) 423 IF ( k > 0 ) THEN 424 MaterialFileName = MaterialFileName(1:k) // '/share/elmersolver/lib/' // 'permafrostmaterialdb.dat' 425 INQUIRE(FILE=TRIM(MaterialFileName), EXIST=fexist) 426 END IF 427 IF ((.NOT. fexist) .AND. k>0) THEN 428 MaterialFileName = MaterialFileName(1:k) // '/permafrostmaterialdb.dat' 429 INQUIRE(FILE=TRIM(MaterialFileName), EXIST=fexist) 430 END IF 431 END IF 432 IF (.NOT. fexist) THEN 433 CALL Fatal('CheckKeyWord', 'permafrostmaterialdb.dat not found') 434 END IF 435 END IF 436 437 ! if we are still here, we open the file (what ever it may be) 438 OPEN(unit = io, file = TRIM(MaterialFileName), status = 'old',iostat = ok) 439 IF (ok /= 0) THEN 440 WRITE(Message,'(A,A)') 'Unable to open file ',TRIM(MaterialFileName) 441 CALL FATAL(Trim(FunctionName),Trim(message)) 442 ELSE 443 !------------------------------------------------------------------------------ 444 ! Read in the number of records in file (first line integer) 445 !------------------------------------------------------------------------------ 446 READ (io, *, END=30, IOSTAT=OK, ERR=40) NumerOfRockRecords, Comment 447 WRITE (Message,*) "Attempting to read ",NumerOfRockRecords," ",& 448 TRIM(Comment)," records from data file ",TRIM(MaterialFileName) 449 CALL INFO(FunctionName,Message,level=3) 450 InitialNumerOfRockRecords = NumerOfRockRecords 451 END IF 452 !------------------------------------------------------------------------------ 453 ! Allocate and read stuff 454 !------------------------------------------------------------------------------ 455 !M = Model % Mesh % NumberOfNodes 456 IF (AllocationsDone) THEN 457 DEALLOCATE(& 458 GlobalRockMaterial % ks0th,& 459 GlobalRockMaterial % e1,& 460 GlobalRockMaterial % bs,& 461 GlobalRockMaterial % rhos0,& 462 GlobalRockMaterial % cs0,& 463 GlobalRockMaterial % Xi0,& 464 GlobalRockMaterial % eta0,& 465 GlobalRockMaterial % etak,& 466 GlobalRockMaterial % hs0,& 467 GlobalRockMaterial % Kgwh0, & 468 GlobalRockMaterial % qexp, & 469 GlobalRockMaterial % alphaL, & 470 GlobalRockMaterial % alphaT, & 471 GlobalRockMaterial % RadGen, & 472 GlobalRockMaterial % acs, & 473 GlobalRockMaterial % as0, & 474 GlobalRockMaterial % aas, & 475 GlobalRockMaterial % ks0, & 476 GlobalRockMaterial % cks, & 477 GlobalRockMaterial % Es0, & 478 GlobalRockMaterial % nuS0, & 479 GlobalRockMaterial % acsl, & 480 GlobalRockMaterial % aasl, & 481 GlobalRockMaterial % cksl, & 482 GlobalRockMaterial % VariableBaseName) 483 END IF 484 ALLOCATE(& 485 GlobalRockMaterial % ks0th(NumerOfRockRecords),& 486 GlobalRockMaterial % e1(NumerOfRockRecords),& 487 GlobalRockMaterial % bs(NumerOfRockRecords),& 488 GlobalRockMaterial % rhos0(NumerOfRockRecords),& 489 GlobalRockMaterial % cs0(NumerOfRockRecords),& 490 GlobalRockMaterial % Xi0(NumerOfRockRecords),& 491 GlobalRockMaterial % eta0(NumerOfRockRecords),& 492 GlobalRockMaterial % etak(NumerOfRockRecords),& 493 GlobalRockMaterial % hs0(NumerOfRockRecords),& 494 GlobalRockMaterial % Kgwh0(3,3,NumerOfRockRecords),& 495 GlobalRockMaterial % qexp(NumerOfRockRecords), & 496 GlobalRockMaterial % alphaL(NumerOfRockRecords), & 497 GlobalRockMaterial % alphaT(NumerOfRockRecords), & 498 GlobalRockMaterial % RadGen(NumerOfRockRecords), & 499 GlobalRockMaterial % acs(0:5,NumerOfRockRecords), & 500 GlobalRockMaterial % as0(NumerOfRockRecords), & 501 GlobalRockMaterial % aas(0:5,NumerOfRockRecords), & 502 GlobalRockMaterial % ks0(NumerOfRockRecords), & 503 GlobalRockMaterial % cks(0:5,NumerOfRockRecords), & 504 GlobalRockMaterial % Es0(NumerOfRockRecords), & 505 GlobalRockMaterial % nuS0(NumerOfRockRecords), & 506 GlobalRockMaterial % acsl(NumerOfRockRecords), & 507 GlobalRockMaterial % aasl(NumerOfRockRecords), & 508 GlobalRockMaterial % cksl(NumerOfRockRecords), & 509 GlobalRockMaterial % VariableBaseName(NumerOfRockRecords),& 510 STAT=OK) 511 AllocationsDone = .TRUE. 512 DataRead = .TRUE. 513 IF (OK /= 0) THEN 514 CLOSE(io) 515 CALL FATAL(FunctionName, 'Allocation Error of input data array') 516 END IF 517 518 DO I=1,NumerOfRockRecords 519 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % VariableBaseName(I), EntryNumber 520 IF (EntryNumber /= I) THEN 521 WRITE(Message,'(A,I3,A,I3)') & 522 "Entry number", EntryNumber, "does not match expected number ",I 523 CLOSE(io) 524 CALL FATAL(FunctionName,Message) 525 ELSE 526 WRITE(Message,'(A,A,A,I3,A)')& 527 "Material ", TRIM(GlobalRockMaterial % VariableBaseName(I)),& 528 " entry number ", EntryNumber, " will be read in" 529 CALL INFO(FunctionName,Message,Level=3) 530 END IF 531 WRITE(Message,'(A,I2,A,A)') "Input for Variable No.",I,": ", GlobalRockMaterial % VariableBaseName(I) 532 CALL INFO(FunctionName,Message,Level=9) 533 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % Xi0(I), Comment 534 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % eta0(I), Comment 535 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % etak(I), Comment 536 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % ks0th(I), Comment 537 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % e1(I), Comment 538 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % bs(I), Comment 539 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % rhos0(I), Comment 540 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % cs0(I), Comment 541 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % hs0(I), Comment 542 DO J=1,3 543 DO K=1,3 544 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % Kgwh0(J,K,I), Comment 545 END DO 546 END DO 547 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % qexp(I), Comment 548 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % alphaL(I), Comment 549 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % alphaT(I), Comment 550 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % RadGen(I), Comment 551 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % acs(0:5,I), Comment 552 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % as0(I), Comment 553 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % aas(0:5,I), Comment 554 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % ks0(I), Comment 555 !-------------------- 556 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % cks(0:5,I), Comment 557 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % Es0(I), Comment 558 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % nuS0(I), Comment 559 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % acsl(I), Comment 560 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % aasl(I), Comment 561 READ (io, *, END=30, IOSTAT=OK, ERR=40) GlobalRockMaterial % cksl(I), Comment 562 END DO 563 WRITE(Message,'(A,I2,A,A)') "Read ",NumerOfRockRecords," rock material records from file ", TRIM(MaterialFileName) 564 CALL INFO(FunctionName,Message,Level=1) 56530 CLOSE(io) 566 IF (I < NumerOfRockRecords) THEN 567 WRITE(Message,'(I3,A,I3)') I,"records read, which is smaller than given number ", NumerOfRockRecords 568 CALL FATAL(FunctionName,Message) 569 ELSE 570 WRITE(Message,'(A,I2,A,A)') "Read ",NumerOfRockRecords," rock material records from file ", TRIM(MaterialFileName) 571 CALL INFO(FunctionName,Message,Level=1) 572 END IF 573 RETURN 574 END IF 575 57640 CALL WARN(FunctionName,"I/O error! Last successfully read variable:") 577 CALL WARN(FunctionName,Comment) 578 CALL FATAL(FunctionName,"Stopping simulation") 579 END FUNCTION ReadPermafrostRockMaterial 580 581 !--------------------------------------------------------------------------------------------- 582 FUNCTION ReadPermafrostElementRockMaterial(MaterialFileName,Solver,DIM,SkipInit) RESULT(NumberOfRockRecords) 583 IMPLICIT NONE 584 CHARACTER(LEN=MAX_NAME_LEN), INTENT(IN) :: MaterialFileName 585 TYPE(Solver_t) :: Solver 586 INTEGER :: NumberOfRockRecords,DIM 587 LOGICAL, OPTIONAL :: SkipInit 588 !----------------------------------------------------------- 589 CHARACTER(LEN=MAX_NAME_LEN) :: SubroutineName="ReadPermafrostElementRockMaterial" 590 LOGICAL :: FirstTime=.TRUE., Parallel=.FALSE. 591 TYPE(RockMaterial_t), TARGET :: LocalRockMaterial 592 TYPE(Element_t), POINTER :: CurrentElement 593 INTEGER, ALLOCATABLE :: GlobalToLocalPerm(:) 594 INTEGER :: OK, CurrentNo, I, J, NoElements, LocalNoElements, & 595 minglobalelementnumber, maxglobalelementnumber, mmaxglobalelementnumber, ierr 596 INTEGER, PARAMETER :: io=21 597 REAL(KIND=dp) :: ReceivingArray(50) 598 599 SAVE LocalRockMaterial, FirstTime, Parallel, minglobalelementnumber, maxglobalelementnumber,& 600 GlobalToLocalPerm, NoElements, mmaxglobalelementnumber 601 602 IF (PRESENT(SkipInit) .AND. FirstTime) CALL FATAL(SubroutineName,'Initialization error') 603 604 IF (FirstTime) THEN 605 NoElements = Solver % NumberOfActiveElements ! active elements in partition/serial mesh 606 Parallel = ( ParEnv % PEs > 1 ) 607 IF ( Parallel ) THEN 608 DO I=1,NoElements 609 CurrentElement => Solver % Mesh % Elements(I) 610 IF (FirstTime) THEN 611 minglobalelementnumber = CurrentElement % GElementIndex 612 maxglobalelementnumber = minglobalelementnumber 613 FirstTime = .FALSE. 614 ELSE 615 minglobalelementnumber = MIN((CurrentElement % GElementIndex),minglobalelementnumber) 616 maxglobalelementnumber = MAX((CurrentElement % GElementIndex),maxglobalelementnumber) 617 END IF 618 END DO 619 CALL MPI_ALLREDUCE(maxglobalelementnumber,mmaxglobalelementnumber,1,& 620 MPI_INTEGER,MPI_MAX,ELMER_COMM_WORLD,ierr) 621 622 !IF (ParEnv % myPE == 0) & 623 ! PRINT *,"ReadPermafrostElementRockMaterial:",Parenv % myPE, "ming/maxg",& 624 ! minglobalelementnumber,maxglobalelementnumber 625 IF ((maxglobalelementnumber - minglobalelementnumber) < 1) & 626 CALL FATAL("ReadPermafrostElementRockMaterial","Failed to create global to local permutation") 627 ALLOCATE(GlobalToLocalPerm(maxglobalelementnumber - minglobalelementnumber + 1), STAT=OK) 628 IF (OK /= 0) CALL FATAL("ReadPermafrostElementRockMaterial","Allocation error of GlobalToLocalPerm") 629 GlobalToLocalPerm = 0 630 DO I=1,NoElements 631 CurrentElement => Solver % Mesh % Elements(I) 632 GlobalToLocalPerm((CurrentElement % GElementIndex) - minglobalelementnumber + 1) = I 633 !IF (ParEnv % myPE == 0) & 634 ! PRINT *,"ReadPermafrostElementRockMaterial:",Parenv % myPE, & 635 ! "GlobalToLocalPerm(",(CurrentElement % GElementIndex)," - ",minglobalelementnumber," + 1)=",I 636 END DO 637 ELSE 638 minglobalelementnumber = 1 639 maxglobalelementnumber = NoElements 640 mmaxglobalelementnumber = NoElements 641 END IF 642 ALLOCATE(& 643 GlobalRockMaterial % ks0th(NoElements),& 644 GlobalRockMaterial % e1(NoElements),& 645 GlobalRockMaterial % bs(NoElements),& 646 GlobalRockMaterial % rhos0(NoElements),& 647 GlobalRockMaterial % Xi0(NoElements),& 648 GlobalRockMaterial % eta0(NoElements),& 649 GlobalRockMaterial % etak(NoElements),& 650 GlobalRockMaterial % hs0(NoElements),& 651 GlobalRockMaterial % Kgwh0(3,3,NoElements),& 652 GlobalRockMaterial % qexp(NoElements), & 653 GlobalRockMaterial % alphaL(NoElements), & 654 GlobalRockMaterial % alphaT(NoElements), & 655 GlobalRockMaterial % RadGen(NoElements), & 656 GlobalRockMaterial % cs0(NoElements),& 657 GlobalRockMaterial % acs(0:5,NoElements), & 658 GlobalRockMaterial % as0(NoElements), & 659 GlobalRockMaterial % aas(0:5,NoElements), & 660 GlobalRockMaterial % ks0(NoElements), & 661 GlobalRockMaterial % cks(0:5,NoElements), & 662 GlobalRockMaterial % Es0(NoElements),& 663 GlobalRockMaterial % nus0(NoElements),& 664 GlobalRockMaterial % acsl(NoElements), & 665 GlobalRockMaterial % aasl(NoElements), & 666 GlobalRockMaterial % cksl(NoElements), & 667 GlobalRockMaterial % VariableBaseName(NoElements),& 668 STAT=OK) 669 OPEN(unit = io, file = TRIM(MaterialFileName), status = 'old',iostat = ok) 670 IF (ok /= 0) THEN 671 WRITE(Message,'(A,A)') 'Unable to open file ',TRIM(MaterialFileName) 672 CALL FATAL(Trim(SubroutineName),Trim(message)) 673 ELSE 674 !------------------------------------------------------------------------------ 675 ! Read in the number of records in file (first line integer) 676 ! MIND: all receiving array numbers are shifted by -1 in index with resepect 677 ! to J. Hartikainen's instructions in input_data_forsmark_2d_example.pdf! 678 !------------------------------------------------------------------------------ 679 WRITE (Message,*) "Attempting to read ",mmaxglobalelementnumber,& 680 " records from data file ",TRIM(MaterialFileName) 681 CALL INFO(SubroutineName,Message,level=3) 682 LocalNoElements = 0 683 DO J=1,maxglobalelementnumber 684 READ (io, *, END=50, ERR=60, IOSTAT=OK) CurrentNo, ReceivingArray(1:50) 685 IF ( Parallel ) THEN 686 IF (J < minglobalelementnumber) CYCLE 687 I = GlobalToLocalPerm(J - minglobalelementnumber +1) 688 !IF (I> 0) & 689 ! PRINT *,"ReadPermafrostElementRockMaterial:", Parenv % myPE, & 690 ! "GlobalToLocalPerm(",J," -", minglobalelementnumber," +1) =", & 691 ! GlobalToLocalPerm(J - minglobalelementnumber +1) 692 IF (I == 0) CYCLE 693 LocalNoElements = LocalNoElements + 1 694 ELSE 695 I=J 696 LocalNoElements = LocalNoElements + 1 697 END IF 698 CurrentElement => Solver % Mesh % Elements(I) 699 !! IMPORTANT: Mind that all ReceivingArray numbers ar N-1 with respect to the document (input_data_forsmark_2d) 700 GlobalRockMaterial % ks0th(I) = ReceivingArray(11) ! shall be changed to tensor 701 !----------------------------- 702 GlobalRockMaterial % e1(I) = ReceivingArray(33) ! e1 (mail from Juha 11.10.) 703 !IF (GlobalRockMaterial % e1(I) > 0.01) PRINT *,"e1:", ReceivingArray(34) 704 !IF (GlobalRockMaterial % e1(I) < 0.0) PRINT *,"e1:", ReceivingArray(34) 705 GlobalRockMaterial % bs(I) = ReceivingArray(23) ! b11,1 (mail from Juha 11.10.) 706 GlobalRockMaterial % rhos0(I) = ReceivingArray(1) 707 GlobalRockMaterial % Xi0(I) = ReceivingArray(32) 708 !----------------------------- 709 GlobalRockMaterial % eta0(I) = ReceivingArray(30) ! eta_t (mail from Juha 11.10.) 710 GlobalRockMaterial % etak(I) = ReceivingArray(31) 711 GlobalRockMaterial % hs0(I) = 0.0_dp! will be removed 712 !----------------------------- Hydrol. Conductivity 713 !GlobalRockMaterial % Kgwh0 = 0.0_dp 714 715 IF(DIM==2) THEN 716 GlobalRockMaterial % Kgwh0(1,1,I) = ReceivingArray(35) 717 GlobalRockMaterial % Kgwh0(2,2,I) = ReceivingArray(37) 718 GlobalRockMaterial % Kgwh0(1,2,I) = ReceivingArray(39) 719 GlobalRockMaterial % Kgwh0(2,1,I) = ReceivingArray(39) 720 GlobalRockMaterial % Kgwh0(3,1:3,I) = 0.0_dp 721 GlobalRockMaterial % Kgwh0(1:3,3,I) = 0.0_dp 722 ELSE 723 GlobalRockMaterial % Kgwh0(1,1,I) = ReceivingArray(35) 724 GlobalRockMaterial % Kgwh0(2,2,I) = ReceivingArray(36) 725 GlobalRockMaterial % Kgwh0(3,3,I) = ReceivingArray(37) 726 GlobalRockMaterial % Kgwh0(1,2,I) = ReceivingArray(38) 727 GlobalRockMaterial % Kgwh0(1,3,I) = ReceivingArray(39) 728 GlobalRockMaterial % Kgwh0(2,3,I) = ReceivingArray(40) 729 GlobalRockMaterial % Kgwh0(2,1,I) = ReceivingArray(38) 730 GlobalRockMaterial % Kgwh0(3,1,I) = ReceivingArray(39) 731 GlobalRockMaterial % Kgwh0(3,2,I) = ReceivingArray(40) 732 END IF 733 !----------------------------- 734 GlobalRockMaterial % qexp(I) = ReceivingArray(41) !????????????????????????????????????????????? 735 GlobalRockMaterial % alphaL(I) = ReceivingArray(47) 736 GlobalRockMaterial % alphaT(I) = ReceivingArray(48) 737 GlobalRockMaterial % RadGen(I) = ReceivingArray(29) 738 !----------------------------- 739 GlobalRockMaterial % cs0(I) = ReceivingArray(8) 740 GlobalRockMaterial % acs(0,I) = ReceivingArray(9) 741 GlobalRockMaterial % acs(1,I) = ReceivingArray(10) 742 GlobalRockMaterial % acs(2:5,I) = 0.0_dp 743 GlobalRockMaterial % acsl(I)= 1 744 !----------------------------- 745 GlobalRockMaterial % as0(I)= ReceivingArray(2) 746 GlobalRockMaterial % aas(0,I) = ReceivingArray(3) 747 GlobalRockMaterial % aas(1,I) = ReceivingArray(4) 748 GlobalRockMaterial % aas(2:5,I) = 0.0_dp 749 GlobalRockMaterial % aasl(I)= 1 750 !----------------------------- 751 GlobalRockMaterial % ks0(I)= ReceivingArray(5) 752 GlobalRockMaterial % cks(0,I) = ReceivingArray(6) 753 GlobalRockMaterial % cks(1,I) = ReceivingArray(7) 754 GlobalRockMaterial % cks(2:5,I)= 0.0_dp 755 GlobalRockMaterial % cksl(I)= 1 756 !----------------------------- 757 GlobalRockMaterial % Es0(I) = ReceivingArray(49) 758 GlobalRockMaterial % nus0(I) = ReceivingArray(50) 759 WRITE(Message,*) 'Element',I 760 GlobalRockMaterial % VariableBaseName(I) = TRIM(Message) 761 END DO 762 GlobalRockMaterial % NumerOfRockRecords = NoElements 763 NumberOfRockRecords = NoElements 76450 CLOSE(io) 765 IF (LocalNoElements < NoElements) THEN 766 IF (Parallel) THEN 767 WRITE (Message,*) 'Parallel proc.', ParEnv % myPe, '. Found ONLY ',& 768 LocalNoElements,' entries in file ',TRIM(MaterialFileName),& 769 ' for ', NoElements, ' elements in mesh partition.' 770 ELSE 771 WRITE (Message,*) 'Found ONLY ',LocalNoElements,' entries in file ',& 772 TRIM(MaterialFileName),& 773 ' for ', NoElements, ' elements in mesh' 774 END IF 775 CALL FATAL(TRIM(SubroutineName),Message) 776 ELSE 777 IF (Parallel) THEN 778 PRINT *, TRIM(SubroutineName), ': Parallel proc.', ParEnv % myPe, '. Read ', NoElements, 'entries from file' 779 ELSE 780 WRITE (Message,*) 'Read ',LocalNoElements,' entries from file ',TRIM(MaterialFileName) 781 CALL INFO(TRIM(SubroutineName),Message,Level=3) 782 END IF 783 END IF 784 END IF 785 IF (Parallel) DEALLOCATE(GlobalToLocalPerm) 786 CALL MPI_BARRIER(ELMER_COMM_WORLD,ierr) 787 FirstTime = .FALSE. 788 ELSE 789 NumberOfRockRecords = NoElements 790 END IF 791 RETURN 79260 WRITE (Message,*) 'I/O error at entry ',CurrentNo,' of file ',TRIM(MaterialFileName) 793 CALL FATAL(TRIM(SubroutineName),Message) 794 END FUNCTION ReadPermafrostElementRockMaterial 795 796 !--------------------------------------------------------------------------------------------- 797 FUNCTION ReadPermafrostConstants(Model, FunctionName,& 798 DIM, GasConstant, N0, DeltaT, T0, p0, eps, Gravity) RESULT(Constantsread) 799 !------------------------------------------------------------------------------ 800 TYPE(Model_t) :: Model 801 CHARACTER(LEN=MAX_NAME_LEN) :: FunctionName 802 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 803 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 804 INTEGER :: DIM 805 REAL(KIND=dp) :: GasConstant, N0, DeltaT, T0, p0,eps, Gravity(3) ! constants read only once 806 LOGICAL :: Constantsread 807 !------------------------------------------------------------------------------ 808 REAL(KIND=dp), POINTER :: gWork(:,:) 809 LOGICAL :: Found 810 INTEGER :: I 811 !------------------------------------------------------------------------------ 812 DIM = CoordinateSystemDimension() 813 gWork => ListGetConstRealArray( Model % Constants,"Gravity",Found) 814 IF (.NOT.Found) THEN 815 Gravity = 0.0 816 CALL WARN(FunctionName,'Gravity not found in Constants section. Setting to zero') 817 ELSE 818 Gravity = gWork(1:3,1)*gWork(4,1) 819 END IF 820 !------------------------------------------------------------------------------ 821 ! Constants 822 ! GasConstant, N0, T0, p0, DeltaT, eps 823 !------------------------------------------------------------------------------ 824 !IF (.NOT.ASSOCIATED(Model % Constants)) STOP 825 GasConstant = GetConstReal(Model % Constants, "Gas Constant", Found) 826 IF (.NOT.Found) THEN 827 GasConstant = 8.3145_dp 828 CALL INFO(FunctionName, ' "Gas Constant" not found in Constants and set to default value 8.3145',Level=3) 829 END IF 830 N0 = GetConstReal(Model % Constants, "Avogadro Number", Found) 831 IF (.NOT.Found) THEN 832 N0 = 6.022140857d23 833 CALL INFO(FunctionName, ' "Avogadro Number" not found in Constants and set to default value 6.022140857E23',Level=3) 834 END IF 835 T0 = GetConstReal(Model % Constants, 'Reference Temperature', Found) 836 IF (.NOT.Found) THEN 837 T0 = 273.15_dp 838 CALL INFO(FunctionName, ' "Reference Temperature" not found in Constants and set to default value T0=273.15',Level=3) 839 END IF 840 p0 = GetConstReal(Model % Constants, 'Reference Pressure', Found) 841 IF (.NOT.Found) THEN 842 p0 = 100132.0_dp 843 CALL INFO(FunctionName, ' "Reference Pressure not found in Constants and set to default value p0=100132.0',Level=3) 844 END IF 845 DeltaT = GetConstReal(Model % Constants,"Permafrost DeltaT",Found) 846 IF (.NOT.Found) THEN 847 DeltaT = 1.0_dp 848 CALL INFO(FunctionName, ' "Permafrost DeltaT" not found in Constants and set to default value DeltaT=1.0',Level=3) 849 END IF 850 Eps = GetConstReal(Model % Constants,"Permafrost eps",Found) 851 IF (.NOT.Found) THEN 852 eps = 0.99_dp 853 CALL INFO(FunctionName, ' "Permafrost eps" not found in Constants and set to default value eps=0.99',Level=3) 854 END IF 855 ConstantsRead = .TRUE. 856 CALL INFO(FunctionName,"-----------------------------------------------------------------",Level=9) 857 CALL INFO(FunctionName,"Model Constants:", Level=9) 858 WRITE(Message,*) "GasConstant, T0, p0, DeltaT, eps:" 859 CALL INFO(FunctionName,Message, Level=9) 860 WRITE(Message,*) GasConstant, T0, p0, DeltaT, eps 861 CALL INFO(FunctionName,Message, Level=9) 862 CALL INFO(FunctionName,"-----------------------------------------------------------------",Level=9) 863 END FUNCTION ReadPermafrostConstants 864 !--------------------------------------------------------------------------------------------- 865 ! assign single nodal variable 866 !--------------------------------------------------------------------------------------------- 867 SUBROUTINE AssignSingleVar(Solver,Model,NodalVariable,VariableVar,VariablePerm,Variable,& 868 VariableName,VariableDOFS,VariableExists,PrevNodalVariable, PrevVariable) 869 IMPLICIT NONE 870 871 TYPE(Solver_t) :: Solver 872 TYPE(Model_t) :: Model 873 REAL(KIND=dp),POINTER :: NodalVariable(:), Variable(:) 874 LOGICAL :: VariableExists 875 CHARACTER(LEN=MAX_NAME_LEN) :: VariableName 876 TYPE(Variable_t), POINTER :: VariableVar 877 INTEGER, POINTER :: VariablePerm(:) 878 INTEGER :: VariableDOFS 879 REAL(KIND=dp),POINTER,OPTIONAL :: PrevNodalVariable(:), PrevVariable(:) 880 881 ! ---- 882 INTEGER :: N, istat 883 CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: SolverName='AssignSingleVar' 884 885 886 !IF (.NOT.ASSOCIATED(VariableVar)) & 887 VariableVar => VariableGet(Solver % Mesh % Variables,VariableName) 888 IF (.NOT.ASSOCIATED(VariableVar)) THEN 889 VariableExists = .FALSE. 890 WRITE (Message,*) 'Variable ',TRIM(VariableName),' not found' 891 CALL WARN(SolverName,Message) 892 RETURN 893 ELSE 894 VariableDOFS = VariableVar % DOFs 895 VariablePerm => VariableVar % Perm 896 Variable => VariableVar % Values 897 IF (.NOT.ASSOCIATED(VariablePerm) .OR. .NOT.ASSOCIATED(Variable)) & 898 CALL FATAL(SolverName, ' Error in assignments of variable pointers') 899 !PRINT *, PRESENT(PrevVariable) 900 IF (PRESENT(PrevVariable)) THEN 901 PrevVariable => VariableVar % PrevValues(:,1) 902 CALL INFO(SolverName,"Assigned previous time values variable pointer",Level=1) 903 END IF 904 END IF 905 906 IF ((.NOT.VariableExists) .OR. (Model % Mesh % Changed)) THEN 907 N = MAX( Solver % Mesh % MaxElementDOFs, Solver % Mesh % MaxElementNodes ) 908 IF (VariableExists) THEN 909 DEALLOCATE(NodalVariable) 910 IF (PRESENT(PrevNodalVariable)) DEALLOCATE(PrevNodalVariable) 911 END IF 912 ALLOCATE(NodalVariable(N*VariableDOFS),STAT=istat ) 913 IF (PRESENT(PrevNodalVariable) .AND. (istat == 0)) ALLOCATE(PrevNodalVariable(N*VariableDOFS),STAT=istat ) 914 IF ( istat /= 0 ) THEN 915 CALL FATAL(SolverName,"Allocation error") 916 ELSE 917 VariableExists = .TRUE. 918 WRITE(Message,*) "Allocations done for nodal variable of ",TRIM(VariableName) 919 CALL INFO(SolverName,Message,Level=1) 920 END IF 921 END IF 922 923 END SUBROUTINE AssignSingleVar 924 !--------------------------------------------------------------------------------------------- 925 ! assign single nodal variable time derivative 926 !--------------------------------------------------------------------------------------------- 927 SUBROUTINE AssignSingleVarTimeDer(Solver,Model,Element,NodalVariableTimeDer,& 928 VariableVar,VariableTimeDerExists,dt) 929 IMPLICIT NONE 930 931 TYPE(Solver_t) :: Solver 932 TYPE(Model_t) :: Model 933 TYPE(Element_t) :: Element 934 REAL(KIND=dp),POINTER :: NodalVariableTimeDer(:) 935 REAL(KIND=dp) :: dt 936 LOGICAL :: VariableTimeDerExists 937 TYPE(Variable_t), POINTER :: VariableVar 938 ! ---- 939 INTEGER :: VariableDOFS 940 !LOGICAL :: AllocationsDone=.FALSE. 941 INTEGER :: MaxNodes, I, J, istat,CurrentvariableNodeIndex 942 INTEGER, POINTER :: VariablePerm(:) 943 CHARACTER(LEN=MAX_NAME_LEN), PARAMETER :: SolverName='AssignSingleVarTimeDer' 944 REAL(KIND=dp), POINTER :: Variable(:),VariablePrev(:,:) 945 946 SAVE MaxNodes 947 948 IF (dt <= 0.0_dp) CALL FATAL(SolverName, "Negative or zero timestep") 949 950 IF ((.NOT.VariableTimeDerExists) .OR. (Model % Mesh % Changed)) THEN 951 MaxNodes = MAX( Solver % Mesh % MaxElementDOFs, Solver % Mesh % MaxElementNodes ) 952 IF (VariableTimeDerExists) THEN 953 CALL INFO(SolverName,"Deallocation of nodal time derivtive") 954 DEALLOCATE(NodalVariableTimeDer) 955 END IF 956 957 VariableDOFS = VariableVar % DOFs 958 ALLOCATE(NodalVariableTimeDer(MaxNodes*VariableDOFS),STAT=istat ) 959 IF ( istat /= 0 ) THEN 960 CALL FATAL(SolverName,"Allocation error") 961 ELSE 962 VariableTimeDerExists = .TRUE. 963 CALL INFO(SolverName,"Allocations Done",Level=1) 964 END IF 965 END IF 966 967 VariableDOFS = VariableVar % DOFs 968 VariablePrev => VariableVar % PrevValues 969 NodalVariableTimeDer(1:maxNodes*VariableDOFS) = 0.0_dp 970 IF ( gettimestep() == 1 ) RETURN ! use zero value in 1st timestep 971 972 IF (.NOT.ASSOCIATED(VariablePrev)) THEN 973 VariableTimeDerExists = .FALSE. 974 ELSE 975 VariablePerm => VariableVar % Perm 976 VariableDOFS = VariableVar % DOFs 977 Variable => VariableVar % Values 978 VariableTimeDerExists = .TRUE. 979 IF (MaxNodes < GetElementNOFNodes(Element)) CALL FATAL(SolverName,"Number of Nodes exceeds allocation") 980 DO I=1,GetElementNOFNodes(Element) 981 CurrentVariableNodeIndex = VariablePerm(Element % NodeIndexes(I)) 982 DO J=1,VariableDOFS 983 NodalVariableTimeDer((I-1)*VariableDOFS + J) = & 984 (Variable((CurrentVariableNodeIndex - 1) * VariableDOFS + J) & 985 - VariablePrev((CurrentVariableNodeIndex - 1) * VariableDOFS + J,1))/dt 986 END DO 987 !PRINT *,"AssignSingleVarTimeDer:", Variable((CurrentVariableNodeIndex - 1) * VariableDOFS + J),& 988 ! VariablePrev((CurrentVariableNodeIndex - 1) * VariableDOFS + J,1) 989 END DO 990 END IF 991 992 END SUBROUTINE AssignSingleVarTimeDer 993 ! assign variables 994 !--------------------------------------------------------------------------------------------- 995 SUBROUTINE AssignVars(Solver,Model,AllocationsDone,& 996 NodalTemperature,NodalPressure,NodalPorosity,NodalSalinity,NodalGWflux, & 997 NodalTemperatureDt,NodalPressureDt,NodalSalinityDt, & 998 TemperatureVar, PressureVar, PorosityVar,SalinityVar, & 999 TemperatureDtVar, PressureDtVar,SalinityDtVar, & 1000 GWFluxVar1,GWFluxVar2,GWFluxVar3, & 1001 TemperaturePerm, PressurePerm, PorosityPerm,SalinityPerm, & 1002 TemperatureDtPerm, PressureDtPerm,SalinityDtPerm, & 1003 GWfluxPerm1, GWfluxPerm2,GWfluxPerm3, & 1004 Temperature, Pressure, Porosity,Salinity,& 1005 TemperatureDt, PressureDt,SalinityDt,& 1006 GWFlux1,GWFlux2,GWFlux3, & 1007 NoPressure, NoSalinity,ConstantPorosity,GivenGWFlux, DIM, ComputeDt,CallerSolverName) 1008 1009 IMPLICIT NONE 1010 1011 TYPE(Solver_t):: Solver 1012 TYPE(Model_t) :: Model 1013 REAL(KIND=dp),POINTER :: NodalTemperature(:),NodalPressure(:),& 1014 NodalPorosity(:),NodalSalinity(:),NodalGWflux(:,:),& 1015 NodalTemperatureDt(:),NodalPressureDt(:),NodalSalinityDt(:) 1016 REAL(KIND=dp),POINTER :: Temperature(:), Pressure(:), Porosity(:), Salinity(:),& 1017 GWflux1(:),GWflux2(:),GWflux3(:),TemperatureDt(:), PressureDt(:),SalinityDt(:) 1018 INTEGER ,POINTER :: TemperaturePerm(:), PressurePerm(:), PorosityPerm(:),SalinityPerm(:),& 1019 GWfluxPerm1(:),GWfluxPerm2(:),GWfluxPerm3(:),& 1020 TemperatureDtPerm(:), PressureDtPerm(:),SalinityDtPerm(:) 1021 TYPE(Variable_t), POINTER :: TemperatureVar, PressureVar, PorosityVar,SalinityVar,& 1022 GWFluxVar1,GWFluxVar2,GWFluxVar3,& 1023 TemperatureDtVar, PressureDtVar,SalinityDtVar 1024 INTEGER :: DIM 1025 LOGICAL :: NoPressure, NoSalinity,AllocationsDone,ConstantPorosity,GivenGWFlux,ComputeDt 1026 CHARACTER(LEN=MAX_NAME_LEN) :: CallerSolverName 1027 !------------------------------ 1028 CHARACTER(LEN=MAX_NAME_LEN) :: TemperatureName,PressureName,PorosityName,SalinityName,& 1029 GWfluxName,SolverName 1030 TYPE(ValueList_t), POINTER :: Params 1031 LOGICAL :: Found 1032 INTEGER :: N, istat 1033 !------------------------------ 1034 1035 SolverName='PermaFrost(AssignVars <-'//TRIM(CallerSolverName)//')' 1036 Params => GetSolverParams() 1037 1038 IF ((.NOT.AllocationsDone) .OR. (Model % Mesh % Changed)) THEN 1039 DIM = CoordinateSystemDimension() 1040 N = MAX( Solver % Mesh % MaxElementDOFs, Solver % Mesh % MaxElementNodes ) 1041 IF (AllocationsDone) & 1042 DEALLOCATE(NodalTemperature,NodalPorosity,NodalPressure,& 1043 NodalSalinity,NodalGWflux,NodalTemperatureDt,NodalPressureDt,& 1044 NodalSalinityDt) 1045 ALLOCATE(NodalTemperature(N),NodalPorosity(N),NodalPressure(N),& 1046 NodalSalinity(N),NodalGWflux(3,N),NodalTemperatureDt(N),& 1047 NodalPressureDt(N),NodalSalinityDt(N),STAT=istat ) 1048 IF ( istat /= 0 ) THEN 1049 CALL FATAL(SolverName,"Allocation error") 1050 ELSE 1051 AllocationsDone = .TRUE. 1052 CALL INFO(SolverName,"Allocations Done",Level=1) 1053 END IF 1054 END IF 1055 1056 IF (TRIM(CallerSolverName) == "PermafrostHeatEquation") THEN 1057 TemperatureVar => Solver % Variable 1058 TemperatureName = Solver % Variable % Name 1059 ELSE 1060 TemperatureName = ListGetString(Params, & 1061 'Temperature Variable', Found ) 1062 IF (.NOT.Found) THEN 1063 CALL WARN(SolverName," 'Temperature Variable' not found. Using default 'Temperature' ") 1064 WRITE(TemperatureName,'(A)') 'Temperature' 1065 ELSE 1066 WRITE(Message,'(A,A)') "'Temperature Variable' found and set to: ", TemperatureName 1067 CALL INFO(SolverName,Message,Level=9) 1068 END IF 1069 TemperatureVar => VariableGet(Solver % Mesh % Variables,TemperatureName) 1070 END IF 1071 IF (.NOT.ASSOCIATED(TemperatureVar)) THEN 1072 WRITE(Message,'(A,A,A)') "'Temperature Variable ", TRIM(TemperatureName), " not associated" 1073 CALL FATAL(SolverName,Message) 1074 ELSE 1075 Temperature => TemperatureVar % Values 1076 TemperaturePerm => TemperatureVar % Perm 1077 WRITE(Message,'(A,A,A)') "'Temperature Variable ", TRIM(TemperatureName), " associated" 1078 CALL INFO(SolverName,Message,Level=9) 1079 IF (ComputeDt .AND. (TRIM(CallerSolverName) == "PermafrostHeatEquation")) THEN 1080 TemperatureDtVar => VariableGet(Solver % Mesh % Variables,TRIM(TemperatureName) // ' Velocity') 1081 IF(.NOT.ASSOCIATED(TemperatureDtVar)) THEN 1082 WRITE (Message,*) ' "Compute Time Derivatives" set to true, but " ', TRIM(TemperatureName), ' Velocity " not found' 1083 CALL WARN(SolverName,Message) 1084 CALL WARN(SolverName,' Switching all time derivatives in source terms off ') 1085 ComputeDt = .FALSE. 1086 ELSE 1087 TemperatureDt => TemperatureDtVar % Values 1088 TemperatureDtPerm => TemperatureDtVar % Perm 1089 END IF 1090 END IF 1091 END IF 1092 1093 1094 1095 1096 IF (TRIM(CallerSolverName) == "PermafrostGroundWaterFlow") THEN 1097 PressureVar => Solver % Variable 1098 PressureName = Solver % Variable % Name 1099 ELSE 1100 PressureName = ListGetString(Params, & 1101 'Pressure Variable', Found ) 1102 IF (.NOT.Found) THEN 1103 CALL WARN(SolverName," 'Pressure Variable' not found. Using default 'Pressure' ") 1104 WRITE(PressureName,'(A)') 'Pressure' 1105 ELSE 1106 WRITE(Message,'(A,A)') "'Pressure Variable' found and set to: ", PressureName 1107 CALL INFO(SolverName,Message,Level=9) 1108 END IF 1109 PressureVar => VariableGet(Solver % Mesh % Variables,PressureName) 1110 END IF 1111 IF (.NOT.ASSOCIATED(PressureVar)) THEN 1112 NULLIFY(Pressure) 1113 NoPressure = .TRUE. 1114 WRITE(Message,'(A,A,A)') "'Pressure Variable ", TRIM(PressureName), " not associated" 1115 CALL WARN(SolverName,Message) 1116 ELSE 1117 Pressure => PressureVar % Values 1118 PressurePerm => PressureVar % Perm 1119 NoPressure = .FALSE. 1120 WRITE(Message,'(A,A,A)') "'Pressure Variable ", TRIM(PressureName), " associated" 1121 CALL INFO(SolverName,Message,Level=9) 1122 IF (ComputeDt .AND. (TRIM(CallerSolverName) == 'PermafrostGroundWaterFlow')) THEN 1123 PressureDtVar => VariableGet(Solver % Mesh % Variables,TRIM(PressureName) // ' Velocity') 1124 IF(.NOT.ASSOCIATED(PressureDtVar)) THEN 1125 WRITE (Message,*) ' "Compute Time Derivatives" set to true, but " ', TRIM(PressureName), ' Velocity " not found' 1126 CALL WARN(SolverName,Message) 1127 CALL WARN(SolverName,' Switching all time derivatives in source terms off ') 1128 ComputeDt = .FALSE. 1129 ELSE 1130 PressureDt => PressureDtVar % Values 1131 PressureDtPerm => PressureDtVar % Perm 1132 END IF 1133 END IF 1134 END IF 1135 1136 PorosityName = ListGetString(Params, & 1137 'Porosity Variable', Found ) 1138 IF (.NOT.Found) THEN 1139 CALL WARN(SolverName," 'Porosity Variable' not found. Using default 'Porosity' ") 1140 WRITE(PorosityName,'(A)') 'Porosity' 1141 ELSE 1142 WRITE(Message,'(A,A)') "'Porosity Variable' found and set to: ", PorosityName 1143 CALL INFO(SolverName,Message,Level=9) 1144 END IF 1145 ConstantPorosity= GetLogical(Params,'Constant Porosity', Found) 1146 IF ((.NOT.Found) .OR. (.NOT.ConstantPorosity)) THEN 1147 PorosityVar => VariableGet(Solver % Mesh % Variables,PorosityName) 1148 IF (.NOT.ASSOCIATED(PorosityVar)) THEN 1149 CALL FATAL(SolverName,'Porosity Variable not found') 1150 ELSE 1151 Porosity => PorosityVar % Values 1152 PorosityPerm => PorosityVar % Perm 1153 END IF 1154 ELSE 1155 NULLIFY(PorosityVar) 1156 END IF 1157 1158 IF (TRIM(CallerSolverName) == 'PermafrostSoluteTransport') THEN 1159 SalinityVar => Solver % Variable 1160 ELSE 1161 SalinityName = ListGetString(Params, & 1162 'Salinity Variable', Found ) 1163 IF (.NOT.Found) THEN 1164 CALL WARN(SolverName," 'Salinity Variable' not found. Using default 'Salinity' ") 1165 WRITE(SalinityName,'(A)') 'Salinity' 1166 ELSE 1167 WRITE(Message,'(A,A)') "'Salinity Variable' found and set to: ", SalinityName 1168 CALL INFO(SolverName,Message,Level=9) 1169 END IF 1170 SalinityVar => VariableGet(Solver % Mesh % Variables,SalinityName) 1171 END IF 1172 IF (.NOT.ASSOCIATED(SalinityVar)) THEN 1173 CALL WARN(SolverName,'Salinity Variable not found. Switching Salinity off') 1174 NoSalinity = .TRUE. 1175 ELSE 1176 Salinity => SalinityVar % Values 1177 SalinityPerm => SalinityVar % Perm 1178 IF (ComputeDt .AND. (TRIM(CallerSolverName) == "PermafrostSoluteTransport")) THEN 1179 SalinityDtVar => VariableGet(Solver % Mesh % Variables,TRIM(SalinityName) // ' Velocity') 1180 IF(.NOT.ASSOCIATED(SalinityDtVar)) THEN 1181 WRITE (Message,*) ' "Compute Time Derivatives" set to true, but " ', TRIM(SalinityName), ' Velocity " not found' 1182 CALL WARN(SolverName,Message) 1183 CALL WARN(SolverName,' Switching all time derivatives in source terms off ') 1184 ComputeDt = .FALSE. 1185 ELSE 1186 SalinityDt => SalinityDtVar % Values 1187 SalinityDtPerm => SalinityDtVar % Perm 1188 END IF 1189 END IF 1190 NoSalinity=.FALSE. 1191 END IF 1192 1193 GWfluxName = ListGetString(Params, & 1194 'Groundwater Flux Variable', GivenGWFlux ) 1195 IF (GivenGWFlux) THEN 1196 WRITE(Message,'(A,A)') "'Groundwater flux Variable' found and set to: ", GWfluxName 1197 CALL INFO(SolverName,Message,Level=9) 1198 GWFluxVar1 => VariableGet(Solver % Mesh % Variables,TRIM(GWfluxName) // " 1") 1199 IF (.NOT.ASSOCIATED(GWFluxVar1)) THEN 1200 PRINT *, TRIM(GWfluxName) // " 1", " not found" 1201 GivenGWflux = .FALSE. 1202 END IF 1203 IF (DIM > 1) THEN 1204 GWFluxVar2 => VariableGet(Solver % Mesh % Variables,TRIM(GWfluxName) // " 2") 1205 IF (.NOT.ASSOCIATED(GWFluxVar2)) THEN 1206 PRINT *, TRIM(GWfluxName) // " 2", " not found" 1207 GivenGWflux = .FALSE. 1208 END IF 1209 IF (DIM > 2) THEN 1210 GWFluxVar3 => VariableGet(Solver % Mesh % Variables,TRIM(GWfluxName) // " 3") 1211 IF (.NOT.ASSOCIATED(GWFluxVar2)) THEN 1212 PRINT *, TRIM(GWfluxName) // " 3", " not found" 1213 GivenGWflux = .FALSE. 1214 END IF 1215 END IF 1216 END IF 1217 GWflux1 => GWFluxVar1 % Values 1218 GWfluxPerm1 => GWFluxVar1 % Perm 1219 IF (DIM > 1) THEN 1220 GWflux2 => GWFluxVar2 % Values 1221 GWfluxPerm2 => GWFluxVar2 % Perm 1222 IF (DIM > 2) THEN 1223 GWflux3 => GWFluxVar3 % Values 1224 GWfluxPerm3 => GWFluxVar3 % Perm 1225 END IF 1226 END IF 1227 CALL INFO(SolverName,'Groundwater flux Variable found. Using this as prescribed groundwater flux',Level=9) 1228 END IF 1229 END SUBROUTINE AssignVars 1230 ! compute element-wise single nodal variable 1231 SUBROUTINE ReadSingleVar(N,Element,VariablePerm,NodalVariable,Variable,VariableDOFs) 1232 IMPLICIT NONE 1233 1234 INTEGER :: N,VariableDOFs 1235 INTEGER, POINTER :: VariablePerm(:) 1236 TYPE(Element_t) :: Element 1237 REAL(KIND=dp),POINTER :: NodalVariable(:),Variable(:) 1238 !----------------------- 1239 INTEGER :: I,J 1240 1241 DO I=1,N 1242 DO J=1,VariableDOFs 1243 NodalVariable((VariableDOFs*I - 1) + J) = & 1244 Variable(VariableDOFs*(VariablePerm(Element % NodeIndexes(I))-1) + J) 1245 END DO 1246 END DO 1247 END SUBROUTINE ReadSingleVar 1248 1249 1250 SUBROUTINE ReadVarsDt(N,Element,Model,Material,& 1251 NodalTemperatureDt,NodalPressureDt,NodalSalinityDt,& 1252 TemperatureDtPerm, PressureDtPerm, SalinityDtPerm,& 1253 TemperatureDt, PressureDt, SalinityDt,& 1254 NoSalinity,NoPressure,CallerSolverName,DIM) 1255 IMPLICIT NONE 1256 1257 INTEGER :: N, DIM 1258 TYPE(Model_t) :: Model 1259 TYPE(Element_t) :: Element 1260 TYPE(ValueList_t), POINTER :: Material 1261 REAL(KIND=dp),POINTER :: NodalTemperatureDt(:),NodalPressureDt(:),NodalSalinityDt(:) 1262 REAL(KIND=dp),POINTER :: TemperatureDt(:), PressureDt(:), SalinityDt(:) 1263 INTEGER ,POINTER :: TemperatureDtPerm(:), PressureDtPerm(:),SalinityDtPerm(:) 1264 CHARACTER(LEN=MAX_NAME_LEN) :: CallerSolverName 1265 LOGICAL :: NoSalinity,NoPressure 1266 1267 IF (.NOT.NoPressure) NodalPressureDt(1:N) = PressureDt(PressureDtPerm(Element % NodeIndexes(1:N))) 1268 IF (.NOT.NoSalinity) NodalSalinityDt(1:N) = SalinityDt(SalinityDtPerm(Element % NodeIndexes(1:N))) 1269 NodalTemperatureDt(1:N) = TemperatureDt(TemperatureDtPerm(Element % NodeIndexes(1:N))) 1270 1271 END SUBROUTINE ReadVarsDt 1272 1273 ! compute element-wise nodal variables 1274 SUBROUTINE ReadVars(N,Element,Model,Material,& 1275 NodalTemperature,NodalPressure,NodalPorosity,NodalSalinity,NodalGWflux,& 1276 Temperature, Pressure, Porosity,Salinity,GWFlux1,GWFlux2,GWFlux3,& 1277 TemperaturePerm, PressurePerm, PorosityPerm,SalinityPerm,& 1278 GWfluxPerm1,GWfluxPerm2,GWfluxPerm3,& 1279 NoSalinity,NoPressure,ConstantPorosity,GivenGWFlux,& 1280 PorosityName,CallerSolverName,DIM) 1281 1282 IMPLICIT NONE 1283 1284 INTEGER :: N, DIM 1285 TYPE(Model_t) :: Model 1286 TYPE(Element_t) :: Element 1287 TYPE(ValueList_t), POINTER :: Material 1288 REAL(KIND=dp),POINTER :: NodalTemperature(:),NodalPressure(:),& 1289 NodalPorosity(:),NodalSalinity(:),NodalGWflux(:,:) 1290 REAL(KIND=dp),POINTER :: Temperature(:), Pressure(:), Porosity(:), Salinity(:),& 1291 GWflux1(:),GWflux2(:),GWflux3(:) 1292 INTEGER ,POINTER :: TemperaturePerm(:), PressurePerm(:), PorosityPerm(:),SalinityPerm(:),& 1293 GWfluxPerm1(:),GWfluxPerm2(:),GWfluxPerm3(:) 1294 LOGICAL :: NoPressure, NoSalinity,ConstantPorosity,GivenGWFlux 1295 CHARACTER(LEN=MAX_NAME_LEN) :: PorosityName, CallerSolverName 1296 !------------------------- 1297 REAL(KIND=dp) :: p0 1298 INTEGER :: I 1299 CHARACTER(LEN=MAX_NAME_LEN) ::SolverName 1300 LOGICAL :: Found 1301 !------------------------- 1302 1303 SolverName='PermaFrost(ReadVars <-'//TRIM(CallerSolverName)//')' 1304 1305 NodalPressure(1:N) = 0.0_dp 1306 NodalSalinity(1:N) = 0.0_dp 1307 NodalGWflux(1:3,1:N) = 0.0_dp 1308 NodalPorosity(1:N) = 0.0_dp 1309 ! Nodal variable dependencies 1310 NodalTemperature(1:N) = Temperature(TemperaturePerm(Element % NodeIndexes(1:N))) 1311 IF (ConstantPorosity) THEN 1312 NodalPorosity(1:N) = ListGetReal(Material,PorosityName,N,Element % NodeIndexes, Found) 1313 IF (.NOT.Found) THEN 1314 WRITE (Message,'(A,A,A)') "No '",TRIM(PorosityName) ,"'found in Material" 1315 CALL FATAL(SolverName,Message) 1316 END IF 1317 ELSE 1318 NodalPorosity(1:N) = Porosity(PorosityPerm(Element % NodeIndexes(1:N))) 1319 END IF 1320 DO I=1,N 1321 IF (NodalPorosity(I) .NE. NodalPorosity(I)) THEN 1322 PRINT *,SolverName,": Invalid value dedected in NodalPorosity" 1323 PRINT *,SolverName,":", Porosity(PorosityPerm(Element % NodeIndexes(1:N))) 1324 CALL FATAL(SolverName,"Exiting") 1325 END IF 1326 END DO 1327 IF (NoPressure) THEN 1328 CALL INFO(SolverName,'No Pressure variable found - setting to "Reference Pressure"',Level=9) 1329 p0 = GetConstReal(Model % Constants, 'Reference Pressure', Found) 1330 IF (.NOT.Found) THEN 1331 p0 = 101032.0_dp 1332 CALL INFO(SolverName, ' "Reference Pressure not found in Constants and set to default value p0=101032.0',Level=9) 1333 END IF 1334 NodalPressure(1:N) = p0 1335 ELSE 1336 NodalPressure(1:N) = Pressure(PressurePerm(Element % NodeIndexes(1:N))) 1337 END IF 1338 IF (NoSalinity) THEN 1339 NodalSalinity(1:N) = 0.0_dp 1340 ELSE 1341 NodalSalinity(1:N) = Salinity(SalinityPerm(Element % NodeIndexes(1:N))) 1342 END IF 1343 IF (GivenGWflux) THEN 1344 NodalGWflux(1,1:N) = & 1345 GWflux1(GWfluxPerm1(Element % NodeIndexes(1:N))) 1346 IF (DIM > 1) THEN 1347 NodalGWflux(2,1:N) = & 1348 GWflux2(GWfluxPerm2(Element % NodeIndexes(1:N))) 1349 IF (DIM > 2) THEN 1350 NodalGWflux(3,1:N) = & 1351 GWflux3(GWfluxPerm3(Element % NodeIndexes(1:N))) 1352 END IF 1353 END IF 1354 END IF 1355 END SUBROUTINE ReadVars 1356 !--------------------------------------------------------------------------------------------- 1357 !--------------------------------------------------------------------------------------------- 1358 !--------------------------------------------------------------------------------------------- 1359 !--------------------------- 1360 ! model parameters functions 1361 !--------------------------- 1362 !--------------------------------------------------------------------------------------------- 1363 !--------------------------------------------------------------------------------------------- 1364 ! 1365 !--------------------------------------------------------------------------------------------- 1366 ! general functions 1367 !--------------------------------------------------------------------------------------------- 1368 REAL(KIND=dp) FUNCTION GeneralPolynomial(Variable,ReferenceValue,Normation,coeff,pdeg) 1369 IMPLICIT NONE 1370 !------- 1371 REAL(KIND=dp), INTENT(IN) :: Variable,ReferenceValue,Normation,coeff(0:5) 1372 INTEGER, INTENT(IN) :: pdeg 1373 REAL(KIND=dp) outval 1374 ! ------ 1375 REAL(KIND=dp) currpot 1376 INTEGER :: i 1377 1378 outval = 0.0_dp 1379 currpot = 1.0_dp 1380 DO i=0,pdeg 1381 outval = outval + coeff(i) * currpot 1382 currpot = currpot * (Variable - ReferenceValue)/Normation 1383 END DO 1384 GeneralPolynomial = outval 1385 END FUNCTION GeneralPolynomial 1386 !--------------------------------------------------------------------------------------------- 1387 REAL(KIND=dp) FUNCTION GeneralIntegral(Variable,ReferenceValue,Normation,coeff0,coeff,pdeg) 1388 IMPLICIT NONE 1389 !------- 1390 REAL(KIND=dp), INTENT(IN) :: Variable,ReferenceValue,Normation,coeff0,coeff(0:5) 1391 INTEGER, INTENT(IN) :: pdeg 1392 REAL(KIND=dp) prefactor, summation 1393 ! ------ 1394 REAL(KIND=dp) currpot 1395 INTEGER :: currdeg 1396 prefactor = coeff0*(Variable - ReferenceValue) 1397 summation = 0.0_dp 1398 DO currdeg=0,pdeg 1399 summation = summation +& 1400 ( coeff(currdeg) * & 1401 ( (Variable - ReferenceValue)/Normation )**(DBLE(currdeg)) )/(DBLE(currdeg) + 1.0_dp) 1402 END DO 1403 !currpot = 1.0_dp 1404 !DO currdeg=0,pdeg 1405 ! outval = outval * coeff(currdeg) * currpot/(DBLE(currdeg) + 1.0_dp) 1406 ! currpot = currpot * (Variable - ReferenceValue)/Normation 1407 !END DO 1408 GeneralIntegral = prefactor * summation 1409 END FUNCTION GeneralIntegral 1410 !--------------------------------------------------------------------------------------------- 1411 ! functions specific to heat transfer and phase change 1412 !--------------------------------------------------------------------------------------------- 1413 FUNCTION GetXiAnderson(A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity) RESULT(XiAnderson) 1414 REAL(KIND=dp), INTENT(IN) :: A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity 1415 REAL(KIND=dp) :: Tstar, XiAnderson 1416 IF (Porosity <= 0.0) & 1417 CALL FATAL("Permafrost(GetXiAnderson)","Zero or negative porosity detected") 1418 Tstar = T0 - Beta * Pressure - Temperature 1419 XiAnderson = MAX(MIN((rhos0/rhow)*(A*(Tstar**B)/Porosity),1.0_dp),0.0_dp) 1420 END FUNCTION GetXiAnderson 1421 !--------------------------------------------------------------------------------------------- 1422 REAL (KIND=dp) FUNCTION XiAndersonT(Xi,A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity) 1423 REAL(KIND=dp), INTENT(IN) :: Xi,A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity 1424 REAL(KIND=dp) :: Tstar 1425 IF (Porosity <= 0.0) & 1426 CALL FATAL("Permafrost(GetXiAndersonT)","Zero or negative porosity detected") 1427 Tstar = T0 - Beta * Pressure - Temperature 1428 IF (Xi == 1.0_dp .OR. Xi == 0.0_dp) THEN 1429 XiAndersonT = 0.0_dp 1430 ELSE 1431 XiAndersonT = -(rhos0/rhow)*(A*B*(Tstar**(B - 1.0_dp)))/Porosity 1432 END IF 1433 END FUNCTION XiAndersonT 1434 !--------------------------------------------------------------------------------------------- 1435 REAL (KIND=dp) FUNCTION XiAndersonP(Xi,A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity) 1436 REAL(KIND=dp), INTENT(IN) :: Xi,A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity 1437 REAL(KIND=dp) :: Tstar 1438 IF (Porosity <= 0.0_dp) & 1439 CALL FATAL("Permafrost(GetXiAndersonT)","Zero or negative porosity detected") 1440 Tstar = T0 - Beta * Pressure - Temperature 1441 IF (Xi == 1_dp .OR. Xi == 0.0_dp) THEN 1442 XiAndersonP = 0.0_dp 1443 ELSE 1444 XiAndersonP = -Beta*(rhos0/rhow)*(A*B*(Tstar**(B - 1.0_dp)))/Porosity 1445 END IF 1446 END FUNCTION XiAndersonP 1447 !--------------------------------------------------------------------------------------------- 1448 REAL (KIND=dp) FUNCTION XiAndersonEta(Xi,A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity) 1449 REAL(KIND=dp), INTENT(IN) :: Xi,A,B,Beta,rhow,rhos0,T0,Temperature,Pressure,Porosity 1450 REAL(KIND=dp) :: Tstar 1451 IF (Porosity <= 0.0) & 1452 CALL FATAL("Permafrost(GetXiAndersonEta)","Zero or negative porosity detected") 1453 Tstar = T0 - Beta * Pressure - Temperature 1454 IF (Xi == 1.0_dp .OR. Xi == 0.0_dp) THEN 1455 XiAndersonEta = 0.0_dp 1456 ELSE 1457 XiAndersonEta = -(rhos0/rhow)*(A*(Tstar**B))/(Porosity*Porosity) 1458 END IF 1459 END FUNCTION XiAndersonEta 1460 !--------------------------------------------------------------------------------------------- 1461 REAL (KIND=dp) FUNCTION delta(CurrentSolventMaterial,& 1462 eps,DeltaT,T0,GasConstant) 1463 IMPLICIT NONE 1464 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1465 REAL(KIND=dp), INTENT(IN) :: eps,DeltaT,T0,GasConstant 1466 REAL(KIND=dp) :: aux,Mw,hi0,cw0,ci0 1467 LOGICAL :: FirstTime=.TRUE. 1468 SAVE FirstTime,Mw,hi0,ci0,cw0 1469 IF(FirstTime) THEN 1470 !assign needed variables 1471 Mw = CurrentSolventMaterial % Mw 1472 hi0 = CurrentSolventMaterial % hi0 1473 ci0 = CurrentSolventMaterial % ci0 1474 cw0 = CurrentSolventMaterial % cw0 1475 FirstTime = .FALSE. 1476 END IF 1477 aux = -0.5_dp*hi0*DeltaT/T0 & 1478 + (cw0 - ci0)*((T0 + 0.5_dp*DeltaT)*LOG(1.0_dp + 0.5_dp*DeltaT/T0) - 0.5_dp*DeltaT) 1479 delta = aux*(eps*(1.0_dp - eps)/(2.0_dp*eps - 1.0_dp))* Mw/(GasConstant*(T0 + 0.5_dp*DeltaT)) 1480 !IF (delta < 1.0d-10) PRINT *, "delta=", delta, "(aux,Mw,T0,DeltaT,eps)",hi0, Mw, T0,DeltaT,eps 1481 END FUNCTION delta 1482 !--------------------------------------------------------------------------------------------- 1483 FUNCTION GetAcAlphatilde(CurrentSolventMaterial,ComputeIce) RESULT(acAlphatilde) 1484 IMPLICIT NONE 1485 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1486 LOGICAL, INTENT(IN) :: ComputeIce 1487 !--------------------------------- 1488 REAL(KIND=dp) :: acAlphaTilde(0:5) 1489 REAL(KIND=dp) :: acAlpha(0:5), sumation 1490 INTEGER :: acAlphal,I 1491 !assign needed properties 1492 acAlphaTilde(0:5) = 0.0_dp 1493 IF (ComputeIce) THEN 1494 acAlphal = CurrentSolventMaterial % acil 1495 acAlpha(0:5) = CurrentSolventMaterial % aci(0:5) 1496 ELSE 1497 acAlphal = CurrentSolventMaterial % acwl 1498 acAlpha(0:5) = CurrentSolventMaterial % acw(0:5) 1499 END IF 1500 ! acAlphal-entries 1501 sumation = 0.0_dp 1502 DO I=acAlphal,1,-1 1503 sumation = acAlpha(I)- sumation 1504 acAlphaTilde(I) = ( (1.0_dp/(DBLE(I) + 1.0_dp)) - 1.0/DBLE(I) ) * sumation 1505 END DO 1506 ! zero-entry only for acAlphaTilde(0) 1507 acAlphaTilde(0)= acAlpha(0)- sumation 1508 END FUNCTION GetAcAlphatilde 1509 !--------------------------------------------------------------------------------------------- 1510 REAL (KIND=dp) FUNCTION gwa(CurrentSolventMaterial,& 1511 p0,T0,rhow,Temperature,Pressure) 1512 IMPLICIT NONE 1513 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1514 REAL(KIND=dp), INTENT(IN) :: p0,T0,rhow,Temperature,Pressure 1515 REAL(KIND=dp) :: cw0,kw0,bcw(0:5) 1516 INTEGER :: I,bcwl 1517 REAL(KIND=dp) :: acwtilde(0:5),aux 1518 LOGICAL :: FirstTime=.TRUE. 1519 1520 SAVE FirstTime,acwtilde 1521 1522 IF (FirstTime) THEN 1523 acwtilde = GetAcAlphatilde(CurrentSolventMaterial,.FALSE.) 1524 FirstTime = .FALSE. 1525 END IF 1526 1527 aux = -(CurrentSolventMaterial % cw0)*(acwtilde(0) * Temperature * LOG(Temperature/T0) & 1528 - (Temperature - T0) & 1529 * GeneralPolynomial(Temperature,T0,T0,acwtilde,CurrentSolventMaterial % acwl)) 1530 gwa = aux + ((Pressure - p0)*(1.0_dp + 0.5_dp*(CurrentSolventMaterial % kw0)*(Pressure - p0))/rhow) 1531 1532 IF (gwa .NE. gwa) THEN 1533 PRINT *, "gwa:", gwa 1534 PRINT *, GeneralPolynomial(Temperature,T0,T0,acwtilde,CurrentSolventMaterial % acwl) 1535 PRINT *, acwtilde(0), Temperature, LOG(Temperature/T0) 1536 PRINT *, ((Pressure - p0)*(1.0_dp + 0.5_dp*(CurrentSolventMaterial % kw0)*(Pressure - p0))/rhow) 1537 PRINT *,rhow, Pressure 1538 CALL FATAL("PermafrostMaterials(gwa)","Error in gwa") 1539 END IF 1540 END FUNCTION gwa 1541 !--------------------------------------------------------------------------------------------- 1542 REAL (KIND=dp) FUNCTION gia(CurrentSolventMaterial,& 1543 p0,T0,rhoi,Temperature,Pressure) 1544 IMPLICIT NONE 1545 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1546 REAL(KIND=dp), INTENT(IN) :: p0,T0,rhoi,Temperature,Pressure 1547 !--------------------------------- 1548 INTEGER :: acil,I 1549 REAL(KIND=dp) ::acitilde(0:5),aux,aux1 1550 LOGICAL :: FirstTime=.TRUE. 1551 !----------------------------------- 1552 SAVE FirstTime,acitilde 1553 1554 IF (FirstTime) THEN 1555 acitilde = GetAcAlphatilde(CurrentSolventMaterial,.TRUE.) 1556 FirstTime = .FALSE. 1557 END IF 1558 aux = -(CurrentSolventMaterial % hi0)*((Temperature - T0)/T0)& 1559 - (CurrentSolventMaterial % ci0) *(acitilde(0) * Temperature * LOG(Temperature/T0) & 1560 - (Temperature - T0) * GeneralPolynomial(Temperature,T0,T0,acitilde,CurrentSolventMaterial % acil)) 1561 gia = aux + ((Pressure - p0)*(1.0_dp + 0.5_dp*(CurrentSolventMaterial % ki0)*(Pressure - p0))/rhoi) 1562 END FUNCTION gia 1563 !--------------------------------------------------------------------------------------------- 1564 REAL (KIND=dp) FUNCTION gwaT(CurrentSolventMaterial,& 1565 p0,T0,rhow,Temperature) 1566 IMPLICIT NONE 1567 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1568 REAL(KIND=dp), INTENT(IN) :: p0,T0,rhow,Temperature 1569 !---------------------------- 1570 INTEGER :: I 1571 REAL(KIND=dp) :: acwtilde(0:5),aux, currpot 1572 LOGICAL :: FirstTime=.TRUE. 1573 1574 SAVE FirstTime, acwtilde 1575 1576 IF (FirstTime) THEN 1577 acwtilde = GetAcAlphatilde(CurrentSolventMaterial,.FALSE.) 1578 FirstTime = .FALSE. 1579 END IF 1580 1581 ! get the derivative 1582 aux = 0.0_dp 1583 currpot = 1.0_dp 1584 DO i=0,CurrentSolventMaterial % acwl 1585 !IF (pdeg == 0) PRINT*,i,"of",pdeg, currpot, Variable, ReferenceValue, Normation 1586 aux = aux + (i + 1.0_dp)*acwtilde(i) * currpot 1587 currpot = currpot * (Temperature - T0)/T0 1588 END DO 1589 gwaT = -(CurrentSolventMaterial % cw0)*(acwtilde(0) * (1.0_dp + LOG(Temperature/T0)) - aux) 1590 1591 ! neglected term 1592 !gwaT = aux & 1593 !* GeneralPolynomial(watercont,1.0_dp,1.0_dp,& 1594 ! CurrentSolventMaterial % bcw(0:5),& 1595 ! CurrentSolventMaterial % bcwl) 1596 1597 IF (gwaT .NE. gwaT) THEN 1598 PRINT *, "gwaT" 1599 CALL FATAL("PermafrostMaterials(gwaT)","Error in gwaT") 1600 END IF 1601 END FUNCTION gwaT 1602 !--------------------------------------------------------------------------------------------- 1603 REAL (KIND=dp) FUNCTION giaT(CurrentSolventMaterial,& 1604 p0,T0,rhoi,Temperature) 1605 IMPLICIT NONE 1606 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1607 REAL(KIND=dp), INTENT(IN) :: p0,T0,rhoi,Temperature 1608 INTEGER :: I 1609 REAL(KIND=dp) :: acitilde(0:5),aux, currpot 1610 LOGICAL :: FirstTime=.TRUE. 1611 1612 SAVE FirstTime, acitilde 1613 1614 IF (FirstTime) THEN 1615 acitilde = GetAcAlphatilde(CurrentSolventMaterial,.TRUE.) 1616 !FirstTime = .FALSE. 1617 END IF 1618 1619 ! get the derivative 1620 aux = 0.0_dp 1621 currpot = 1.0_dp 1622 DO i=0,CurrentSolventMaterial % acil 1623 aux = aux + (i + 1.0_dp)*acitilde(i) * currpot 1624 currpot = currpot * (Temperature - T0)/T0 1625 END DO 1626 aux = (CurrentSolventMaterial % ci0)*(acitilde(0) * (1.0_dp + LOG(Temperature/T0)) - aux) 1627 giaT = -(CurrentSolventMaterial % hi0)/T0 - aux 1628 IF (giaT .NE. giaT) THEN 1629 PRINT *, "giaT" 1630 CALL FATAL("PermafrostMaterials(gwa)","Error in giaT") 1631 END IF 1632 FirstTime = .FALSE. 1633 END FUNCTION giaT 1634 !--------------------------------------------------------------------------------------------- 1635 REAL (KIND=dp) FUNCTION deltaG(gwa,gia) 1636 IMPLICIT NONE 1637 REAL(KIND=dp), INTENT(IN) :: gwa,gia 1638 deltaG = gwa - gia 1639 END FUNCTION deltaG 1640 !--------------------------------------------------------------------------------------------- 1641 FUNCTION GetBi(CurrentSoluteMaterial,RockMaterialID,& 1642 Xi0Tilde,Salinity,Update) RESULT(bi) 1643 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 1644 REAL(KIND=dp), INTENT(IN) :: Xi0Tilde,Salinity 1645 INTEGER, INTENT(IN) :: RockMaterialID 1646 REAL(KIND=dp):: bi(4) 1647 LOGICAL :: Update 1648 !---------- 1649 REAL(KIND=dp):: aux,d1,d2,e1 1650 1651 IF (Update) THEN 1652 e1 = GlobalRockMaterial % e1(RockMaterialID) 1653 bi(3) = (1.0_dp - Xi0Tilde)*e1 1654 bi(4) = Xi0Tilde*e1 1655 ELSE 1656 d1 = CurrentSoluteMaterial % d1 1657 d2 = CurrentSoluteMaterial % d2 1658 aux = Salinity/(1.0_dp - Salinity) 1659 bi(1) = aux*(d1 + 0.5_dp*d2*aux) 1660 bi(2) = aux*(d1 + d2*aux)/(1.0_dp - Salinity) 1661 bi(3) = 0.0_dp 1662 bi(4) = 0.0_dp 1663 END IF 1664 END FUNCTION GetBi 1665 !--------------------------------------------------------------------------------------------- 1666 FUNCTION GetBiYc(CurrentSoluteMaterial,Salinity) RESULT(biYc) 1667 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 1668 REAL(KIND=dp), INTENT(IN) :: Salinity 1669 REAL(KIND=dp):: biYc(2) 1670 !---------- 1671 REAL(KIND=dp):: aux,d1,d2 1672 1673 d1 = CurrentSoluteMaterial % d1 1674 d2 = CurrentSoluteMaterial % d2 1675 1676 aux = 1.0_dp/(1.0_dp - Salinity) 1677 biYc(1) = (d1 + d2*Salinity*aux)*aux*aux 1678 biYc(2) = (d1*(1.0_dp + Salinity) + d2*Salinity*(2.0_dp + Salinity))*aux**3.0_dp 1679 END FUNCTION GetBiYc 1680 !--------------------------------------------------------------------------------------------- 1681 FUNCTION GetB(RockMaterialID,CurrentSolventMaterial,& 1682 Xi0tilde,delta,deltaG,GasConstant,bi,Temperature) RESULT(B) 1683 IMPLICIT NONE 1684 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1685 REAL(KIND=dp), INTENT(IN) :: Xi0tilde,delta,deltaG,GasConstant,bi(4),Temperature 1686 INTEGER, INTENT(IN) :: RockMaterialID 1687 REAL(KIND=dp) :: B 1688 REAL(KIND=dp) :: e1,Mw 1689 Mw = CurrentSolventMaterial % Mw 1690 e1 = GlobalRockMaterial % e1(RockMaterialID) 1691 1692 B =(Mw*deltaG/(GasConstant*Temperature) - bi(1) + bi(3))/(delta + bi(2) + bi(4)) 1693 1694 IF (B .NE. B) THEN 1695 PRINT *, "B:", Mw, deltaG,Temperature,bi(1),e1,delta,bi(2),bi(4) 1696 CALL FATAL("PermafrostMaterials(getB)","Error in getB") 1697 END IF 1698 END FUNCTION GetB 1699 !--------------------------------------------------------------------------------------------- 1700 REAL (KIND=dp) FUNCTION D(RockMaterialID,delta,bi) 1701 IMPLICIT NONE 1702 INTEGER, INTENT(IN) :: RockMaterialID 1703 REAL(KIND=dp), INTENT(IN) :: delta,bi(4) 1704 ! local 1705 D = delta/(delta + bi(2) + bi(4)) 1706 !IF (D .NE. D) THEN 1707 ! PRINT *, "D" 1708 ! CALL FATAL("PermafrostMaterials(getD)","Error in getD") 1709 !END IF 1710 END FUNCTION D 1711 !--------------------------------------------------------------------------------------------- 1712 FUNCTION GetXi0Tilde(RockMaterialID,Porosity) RESULT(Xi0tilde) 1713 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1714 INTEGER, INTENT(IN) :: RockMaterialID 1715 REAL(KIND=dp), INTENT(IN) :: Porosity 1716 REAL(KIND=dp) :: Xi0tilde 1717 REAL(KIND=dp) :: Xi0,eta0 1718 LOGICAL :: FirstTime = .TRUE. 1719 SAVE FirstTime 1720 1721 Xi0 = GlobalRockMaterial % Xi0(RockMaterialID) 1722 eta0 = GlobalRockMaterial % eta0(RockMaterialID) 1723 IF (Porosity <= 0.0_dp) THEN 1724 IF (Xi0 == 0.0_dp) THEN 1725 Xi0tilde = 1.0_dp 1726 ELSE 1727 CALL FATAL("Permafrost(GetXi)","Zero or negative porosity detected") 1728 END IF 1729 ELSE 1730 !Xi0tilde = MIN(Xi0 * (eta0/Porosity) * (1.0_dp - Porosity)/(1.0_dp - eta0),1.0_dp) 1731 Xi0tilde = Xi0 * (eta0/Porosity) * (1.0_dp - Porosity)/(1.0_dp - eta0) 1732 END IF 1733 END FUNCTION GetXi0Tilde 1734 !--------------------------------------------------------------------------------------------- 1735 REAL(KIND=dp) FUNCTION fw(RockMaterialID,CurrentSolventMaterial,& 1736 Xi0tilde,rhow,Xi,GasConstant,Temperature) 1737 IMPLICIT NONE 1738 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1739 INTEGER, INTENT(IN) :: RockMaterialID 1740 REAL(KIND=dp), INTENT(IN) :: Xi0tilde,rhow,Xi,GasConstant,Temperature 1741 !----------------- 1742 REAL(KIND=dp) :: Mw, e1 1743 !---------------- 1744 ! IF (Xi > Xi0tilde) THEN 1745 ! fw = 0.0_dp 1746 ! ELSE 1747 ! e1 = GlobalRockMaterial % e1(RockMaterialID) 1748 ! Mw = CurrentSolventMaterial % Mw 1749 ! fw = rhow*GasConstant*Temperature*e1*Xi0tilde/(Mw*Xi) 1750 ! END IF 1751 fw = 0.0_dp !! CHANGE BACK WHEN JUHA TELLS US TO DO SO !! 1752 END FUNCTION fw 1753 !--------------------------------------------------------------------------------------------- 1754 FUNCTION GetXi(B,D) RESULT(Xi) 1755 REAL(KIND=dp), INTENT(IN) :: B,D 1756 REAL(KIND=dp) :: Xi 1757 Xi= 1.0_dp/(1.0_dp + 0.5_dp*B + SQRT(0.25_dp*B*B + D)) 1758 END FUNCTION GetXi 1759 !--------------------------------------------------------------------------------------------- 1760 REAL (KIND=dp) FUNCTION XiT(CurrentSolventMaterial,& 1761 B,D,Xi,bi,p0,delta,deltaG,T0,gwa,gia,gwaT,giaT,GasConstant,Temperature) 1762 IMPLICIT NONE 1763 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1764 REAL(KIND=dp), INTENT(IN) :: B,D,Xi,bi(4),p0,delta,deltaG,& 1765 T0,gwa,gia,gwaT,giaT,GasConstant,Temperature 1766 !local 1767 REAL(KIND=dp) :: aux1, aux2, aux3, Mw,e1, hi0,hw0,rhow0,rhoi0,cw0,ci0 1768 LOGICAL :: FirstTime=.TRUE. 1769 1770 SAVE Mw,hi0,hw0,rhow0,rhoi0,cw0,ci0, FirstTime 1771 1772 IF (FirstTime) THEN 1773 Mw = CurrentSolventMaterial % Mw 1774 hi0 = CurrentSolventMaterial % hi0 1775 hw0 = CurrentSolventMaterial % hw0 1776 rhow0 = CurrentSolventMaterial % rhow0 1777 rhoi0 = CurrentSolventMaterial % rhoi0 1778 cw0 = CurrentSolventMaterial % cw0 1779 ci0 = CurrentSolventMaterial % ci0 1780 FirstTime=.FALSE. 1781 END IF 1782 aux1 = 1.0_dp/(delta + bi(2) + bi(4)) 1783 aux2 = (1.0_dp + B/SQRT(B*B + 4.0_dp*D)) 1784 aux3 = ((gwa - gia)/Temperature - (gwaT - giaT)) 1785 XiT = 0.5_dp*(Mw/(GasConstant*Temperature))*aux1*aux2*aux3*Xi*Xi 1786 END FUNCTION XiT 1787 !--------------------------------------------------------------------------------------------- 1788 REAL (KIND=dp) FUNCTION XiP(CurrentSolventMaterial,& 1789 B,D,bi,Xi,gwap,giap,delta,GasConstant,Temperature) 1790 IMPLICIT NONE 1791 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1792 REAL(KIND=dp), INTENT(IN) :: B,D,bi(4),Xi,gwap,giap,delta,GasConstant,Temperature 1793 !local 1794 REAL(KIND=dp) :: aux1, aux2, rhow0,rhoi0, Mw 1795 LOGICAL :: FirstTime=.TRUE. 1796 1797 SAVE Mw,rhow0,rhoi0,FirstTime 1798 IF (FirstTime) THEN 1799 Mw = CurrentSolventMaterial % Mw 1800 rhow0 = CurrentSolventMaterial % rhow0 1801 rhoi0 = CurrentSolventMaterial % rhoi0 1802 FirstTime=.FALSE. 1803 END IF 1804 IF (Temperature <= 0.0_dp) CALL FATAL("Permafrost (XiP)","(sub-)Zero Temperature detected") 1805 aux1 = 1.0_dp/(delta + bi(2) + bi(4)) 1806 aux2 = (1.0_dp + B/SQRT(B*B + 4.0_dp*D)) 1807 XiP = 0.5_dp * aux1 * aux2 *(giap - gwap)* Mw/(GasConstant*Temperature)*Xi*Xi 1808 END FUNCTION XiP 1809 !--------------------------------------------------------------------------------------------- 1810 REAL (KIND=dp) FUNCTION XiYc(B,D,bi,biYc,Xi,delta) 1811 IMPLICIT NONE 1812 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 1813 REAL(KIND=dp), INTENT(IN) :: B,D,bi(4),biYc(2),Xi,delta 1814 !local 1815 REAL(KIND=dp) :: aux1, aux2, aux3, aux_sqrt 1816 1817 aux_sqrt = B*B + 4.0_dp*D 1818 aux1 = 1.0_dp/(delta + bi(2) + bi(4)) 1819 aux2 = ( 1.0_dp + B/SQRT(aux_sqrt) )*(biYc(1) + B*biYc(2)) & 1820 + biYc(2)/(SQRT(aux_sqrt)) 1821 aux3 = 2.0_dp*D*biYc(2)/(SQRT(aux_sqrt)) 1822 XiYc = 0.5_dp*aux1*(aux2 + aux3)*Xi*Xi 1823 !IF (XiYc .NE. XiYc) THEN 1824 ! PRINT *, "XiYc:", aux1, aux2 1825 ! PRINT *, B, D, biYc(1), biYc(2),delta 1826 ! STOP 1827 !END IF 1828 END FUNCTION XiYc 1829 !--------------------------------------------------------------------------------------------- 1830 REAL (KIND=dp) FUNCTION XiEta(RockMaterialID,& 1831 B,D,bi,biYc,Xi,delta,Porosity) 1832 IMPLICIT NONE 1833 INTEGER, INTENT(IN) :: RockMaterialID 1834 REAL(KIND=dp), INTENT(IN) :: B,D,bi(4),biYc(2),Xi,delta,Porosity 1835 !local 1836 REAL(KIND=dp) :: aux1, aux2, aux3, aux_sqrt,Xi0,eta0 1837 1838 Xi0 = GlobalRockMaterial % Xi0(RockMaterialID) 1839 eta0 = GlobalRockMaterial % eta0(RockMaterialID) 1840 1841 aux_sqrt = B*B + 4.0_dp*D 1842 1843 IF (Porosity >= 0.0_dp) THEN 1844 aux1 = 1.0_dp/(delta + bi(2) + bi(4)) 1845 aux2 = ( 1.0_dp + B/SQRT(aux_sqrt) )*(1.0_dp + B) 1846 aux3 = 2.0_dp*D*biYc(2)/SQRT(aux_sqrt) 1847 XiEta = 0.5_dp*aux1*(aux2 + aux3) * (Xi0*eta0/(1.0_dp - eta0))& 1848 *(1.0_dp/(Porosity**2.0_dp))*Xi*Xi 1849 ELSE 1850 CALL WARN("Permafrost(XiEta)","Porosity out of physical range - returning zero") 1851 XiEta = 0.0_dp 1852 END IF 1853 END FUNCTION XiEta 1854 !---------------------------------------------------------------------- 1855 SUBROUTINE GetXiHartikainen (RockMaterialID,& 1856 CurrentSoluteMaterial,CurrentSolventMaterial,& 1857 TemperatureAtIP,PressureAtIP,SalinityAtIP,PorosityAtIP,& 1858 Xi0tilde,deltaInElement,rhowAtIP,rhoiAtIP,& 1859 GasConstant,p0,T0,& 1860 XiAtIP,XiTAtIP,XiYcAtIP,XiPAtIP,XiEtaAtIP,& 1861 ComputeXi,ComputeXiT, ComputeXiYc, ComputeXiP, ComputeXiEta) 1862 1863 IMPLICIT NONE 1864 1865 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 1866 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1867 INTEGER :: RockMaterialID 1868 REAL(KIND=dp), INTENT(IN) :: Xi0tilde,deltaInElement,rhowAtIP,rhoiAtIP 1869 REAL(KIND=dp), INTENT(IN) :: GasConstant,p0,T0 1870 REAL(KIND=dp), INTENT(IN) :: TemperatureAtIP,PressureAtIP,SalinityAtIP,PorosityAtIP 1871 REAL(KIND=dp), INTENT(OUT) :: XiAtIP,XiTAtIP,XiYcAtIP,XiPAtIP,XiEtaAtIP 1872 LOGICAL, INTENT(IN) :: ComputeXi,ComputeXiT, ComputeXiYc, ComputeXiP, ComputeXiEta 1873 !--------------------------- 1874 REAL(KIND=dp) :: biAtIP(4),biYcAtIP(2),gwaAtIP,gwaTAtIP,gwapAtIP,& 1875 giaAtIP,giaTAtIP,giapAtIP,deltaGAtIP,DAtIP,BAtIP 1876 !--------------------------- 1877 IF (ComputeXi .OR. (ComputeXiT .OR. ComputeXiYC .OR. ComputeXiP)) THEN 1878 biAtIP = GetBi(CurrentSoluteMaterial,RockMaterialID,& 1879 Xi0Tilde,SalinityAtIP,.FALSE.) 1880 gwaAtIP = gwa(CurrentSolventMaterial,& 1881 p0,T0,rhowAtIP,TemperatureAtIP,PressureAtIP) 1882 giaAtIP = gia(CurrentSolventMaterial,& 1883 p0,T0,rhoiAtIP,TemperatureAtIP,PressureAtIP) 1884 deltaGAtIP = deltaG(gwaAtIP,giaAtIP) 1885 DAtIP= D(RockMaterialID,deltaInElement,biAtIP) 1886 BAtIP = GetB(RockMaterialID,CurrentSolventMaterial,& 1887 Xi0tilde,deltaInElement,deltaGAtIP,GasConstant,biAtIP,TemperatureAtIP) 1888 ELSE 1889 CALL WARN("GetXiHartikainen","Nothing to be done - why did you call this routine?") 1890 END IF 1891 IF (ComputeXi) THEN 1892 XiAtIP = GetXi(BAtIP,DAtIP) 1893 END IF 1894 1895 ! updates of derivatives 1896 IF (XiAtIP < Xi0tilde) THEN 1897 biAtIP = GetBi(CurrentSoluteMaterial,RockMaterialID,& 1898 Xi0Tilde,SalinityAtIP,.TRUE.) 1899 XiAtIP = GetXi(BAtIP,DAtIP) 1900 END IF 1901 !---------------------------------------------------- 1902 XiTAtIP = 0.0_dp 1903 XiYcAtIP = 0.0_dp 1904 XiPAtIP = 0.0_dp 1905 IF (ComputeXiT) THEN 1906 giaTAtIP = giaT(CurrentSolventMaterial,& 1907 p0,T0,rhoiAtIP,TemperatureAtIP) 1908 gwaTAtIP = gwaT(CurrentSolventMaterial,& 1909 p0,T0,rhowAtIP,TemperatureAtIP)! 1910 XiTAtIP= XiT(CurrentSolventMaterial,& 1911 BAtIP,DAtIP,XiAtIP,biAtIP,p0,& 1912 deltaInElement,deltaGAtIP,T0,gwaAtIP,giaAtIP,gwaTAtIP,giaTAtIP,GasConstant,TemperatureAtIP) 1913 END IF 1914 IF (ComputeXiYC) THEN 1915 biYcAtIP = GetBiYc(CurrentSoluteMaterial,SalinityAtIP) 1916 XiYcAtIP = XiYc(BAtIP,DAtIP,biAtIP,biYcAtIP,XiAtIP,deltaInElement) 1917 END IF 1918 IF (ComputeXiP) THEN 1919 giapAtIP = 1.0_dp/rhoiAtIP 1920 gwapAtIP = 1.0_dp/rhowAtIP 1921 XiPAtIP = XiP(CurrentSolventMaterial,& 1922 BAtIP,DAtIP,biAtIP,gwapAtIP,giapAtIP,XiAtIP,& 1923 deltaInElement,GasConstant,TemperatureAtIP) 1924 END IF 1925 END SUBROUTINE GetXiHartikainen 1926 !--------------------------------------------------------------------------------------------- 1927 ! Densities and their derivatives, thermal expansion, isothermal chemical compaction and 1928 ! compressibility coefficients 1929 !--------------------------------------------------------------------------------------------- 1930 REAL (KIND=dp) FUNCTION rhos(RockMaterialID,& 1931 T0,p0,Temperature,Pressure,ConstVal) 1932 IMPLICIT NONE 1933 INTEGER, INTENT(IN) :: RockMaterialID 1934 REAL(KIND=dp), INTENT(IN) :: T0,p0,Temperature,Pressure 1935 LOGICAL :: ConstVal 1936 !---------------------- 1937 REAL(KIND=dp) :: aux1,aux2 1938 !---------------------- 1939 rhos = GlobalRockMaterial % rhos0(RockMaterialID) 1940 IF (.NOT.ConstVal) THEN 1941 !aux1 = GeneralIntegral(Pressure,p0,p0,& 1942 ! GlobalRockMaterial % ks0(RockMaterialID),& 1943 ! GlobalRockMaterial % cks(0:5,RockMaterialID),& 1944 ! GlobalRockMaterial % cksl(RockMaterialID)) 1945 ! a shortcut, as only cks(0) = 1 1946 aux1 = (GlobalRockMaterial % ks0(RockMaterialID)) * (Pressure - p0) 1947 aux2 = GeneralIntegral(Temperature,T0,T0,& 1948 GlobalRockMaterial % as0(RockMaterialID),& 1949 GlobalRockMaterial % aas(0:5,RockMaterialID),& 1950 GlobalRockMaterial % aasl(RockMaterialID)) 1951 rhos = rhos * EXP(aux1 - aux2) 1952 END IF 1953 END FUNCTION rhos 1954 !--------------------------------------------------------------------------------------------- 1955 REAL(KIND=dp) FUNCTION rhosT(RockMaterialID,rhos,T0,Temperature) 1956 IMPLICIT NONE 1957 REAL(KIND=dp), INTENT(IN) :: rhos,T0,Temperature 1958 INTEGER, INTENT(IN) :: RockMaterialID 1959 REAL(KIND=dp) :: alphaS 1960 1961 alphaS = GlobalRockMaterial % as0(RockMaterialID) *& 1962 GeneralPolynomial(Temperature,T0,T0,& 1963 GlobalRockMaterial % aas(0:5,RockMaterialID),& 1964 GlobalRockMaterial % aasl(RockMaterialID)) 1965 rhosT = rhos * alphaS 1966 END FUNCTION rhosT 1967!--------------------------------------------------------------------------------------------- 1968 REAL(KIND=dp) FUNCTION rhosp(RockMaterialID,rhos,p0,Pressure) 1969 IMPLICIT NONE 1970 INTEGER, INTENT(IN) :: RockMaterialID 1971 REAL(KIND=dp), INTENT(IN) :: rhos,p0,Pressure 1972 !-------------------- 1973 REAL(KIND=dp) :: kappas 1974 !-------------------- 1975 kappas = ( GlobalRockMaterial % ks0(RockMaterialID)) 1976 rhosP = rhos * kappas 1977 END FUNCTION rhosp 1978 !--------------------------------------------------------------------------------------------- 1979 REAL (KIND=dp) FUNCTION rhow(CurrentSolventMaterial,T0,p0,Temperature,Pressure,ConstVal) 1980 IMPLICIT NONE 1981 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 1982 REAL(KIND=dp), INTENT(IN) :: T0,p0,Temperature,Pressure 1983 LOGICAL :: ConstVal 1984 !---------------------- 1985 REAL(KIND=dp) :: aux1, aux11, aux2, aux22, aux3, aux33, watercont 1986 !---------------------- 1987 rhow = CurrentSolventMaterial % rhow0 1988 IF (.NOT.ConstVal) THEN 1989 !aux1 = GeneralIntegral(Pressure,p0,p0,& 1990 ! CurrentSolventMaterial % kw0,& 1991 ! CurrentSolventMaterial % ckw(0:5),& 1992 ! CurrentSolventMaterial % ckwl) 1993 ! a shortcut, as only ckw(0) = 1 1994 aux1 = (CurrentSolventMaterial % kw0) * (Pressure - p0) 1995 aux2 = GeneralIntegral(Temperature,T0,T0,& 1996 CurrentSolventMaterial % aw0,& 1997 CurrentSolventMaterial % aaw(0:5),& 1998 CurrentSolventMaterial % aawl) 1999 rhow = rhow * EXP(aux1 - aux2) 2000 IF (rhow < 800.0) THEN 2001 PRINT *, "rhow:", rhow,CurrentSolventMaterial % rhow0,aux1, aux2,Pressure,Temperature 2002 END IF 2003 IF (rhow .NE. rhow) THEN 2004 PRINT *, "rhow:", rhow,CurrentSolventMaterial % rhow0,aux1, aux2,Pressure,Temperature 2005 CALL FATAL("PermafrostMaterials(rhow)","Error in rhow") 2006 END IF 2007 END IF 2008 END FUNCTION rhow 2009 !--------------------------------------------------------------------------------------------- 2010 REAL (KIND=dp) FUNCTION rhowupdate(CurrentSolventMaterial,& 2011 previousrhow,Xi,Salinity,ConstVal) 2012 IMPLICIT NONE 2013 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2014 REAL(KIND=dp), INTENT(IN) :: previousrhow,Xi,Salinity 2015 LOGICAL :: ConstVal 2016 !---------------------- 2017 REAL(KIND=dp) :: aux3, aux33, watercont 2018 !---------------------- 2019 rhowupdate = previousrhow 2020 IF (.NOT.ConstVal) THEN 2021 watercont = MAX(0.0_dp, 1.0_dp - Salinity/Xi) 2022 aux3 = GeneralIntegral(watercont,1.0_dp,1.0_dp,& 2023 CurrentSolventMaterial % zw0,& 2024 CurrentSolventMaterial % bzw(0:5),& 2025 CurrentSolventMaterial % bzwl) 2026 !aux33 = (watercont - 1.0_dp)* (CurrentSolventMaterial % zw0) *& 2027 ! ( (CurrentSolventMaterial % bzw(0))& 2028 ! + 0.5_dp*((CurrentSolventMaterial % bzw(1)) * (watercont - 1.0_dp))) 2029 rhowupdate = previousrhow * EXP(aux3) 2030 !IF (aux3 .NE. aux33) THEN 2031 ! PRINT *, "rhowupdate:", previousrhow, EXP(aux3), aux3, aux33, watercont,Salinity,Xi 2032 ! PRINT *, "zw0",CurrentSolventMaterial % zw0, "bzw",CurrentSolventMaterial % bzw(0:CurrentSolventMaterial % bzwl) 2033 ! PRINT *, CurrentSolventMaterial % bzwl, CurrentSolventMaterial % zw0, "---", CurrentSolventMaterial % bzw(0:5) 2034 !END IF 2035 !IF (rhowupdate .NE. rhowupdate) THEN 2036 ! PRINT *, "rhowupdate:" 2037 ! STOP 2038 !END IF 2039 END IF 2040 END FUNCTION rhowupdate 2041 !--------------------------------------------------------------------------------------------- 2042 REAL(KIND=dp) FUNCTION rhowT(CurrentSolventMaterial,rhow,T0,Temperature) 2043 IMPLICIT NONE 2044 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2045 REAL(KIND=dp), INTENT(IN) :: rhow,T0,Temperature 2046 !-------------------- 2047 REAL(KIND=dp) :: alphaW 2048 !-------------------- 2049 alphaW = (CurrentSolventMaterial % aw0) *& 2050 GeneralPolynomial(Temperature,T0,T0,& 2051 CurrentSolventMaterial % aaw(0:5),& 2052 CurrentSolventMaterial % aawl) 2053 rhowT = rhow * alphaW 2054 END FUNCTION rhowT 2055 !--------------------------------------------------------------------------------------------- 2056 REAL(KIND=dp) FUNCTION rhowP(CurrentSolventMaterial,rhow,p0,Pressure) 2057 IMPLICIT NONE 2058 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2059 REAL(KIND=dp), INTENT(IN) :: rhow,p0,Pressure 2060 !-------------------- 2061 REAL(KIND=dp) :: kappaW 2062 !-------------------- 2063 !kappaW = (CurrentSolventMaterial % kw0) *& 2064 ! GeneralPolynomial(Pressure,p0,p0,& 2065 ! CurrentSolventMaterial % ckw(0:5),& 2066 ! CurrentSolventMaterial % ckwl) 2067 ! a shortcut, as only ckw(0) = 1.0 2068 kappaW = (CurrentSolventMaterial % kw0) 2069 rhowP = rhow * kappaW 2070 END FUNCTION rhowP 2071 !--------------------------------------------------------------------------------------------- 2072 REAL(KIND=dp) FUNCTION rhowYc(CurrentSolventMaterial,rhow,Xi,Salinity) 2073 IMPLICIT NONE 2074 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2075 REAL(KIND=dp), INTENT(IN) :: rhow,Xi,Salinity 2076 !-------------------- 2077 REAL(KIND=dp) :: zetaW, xw 2078 !-------------------- 2079 xw = 1.0_dp - (Salinity/Xi) 2080 zetaW = (CurrentSolventMaterial % zw0) *& 2081 GeneralPolynomial(xw,1.0_dp,1.0_dp,& 2082 CurrentSolventMaterial % bzw(0:5),& 2083 CurrentSolventMaterial % bzwl) 2084 rhowYc = -zetaW *rhow/Xi 2085 END FUNCTION rhowYc 2086 !--------------------------------------------------------------------------------------------- 2087 REAL (KIND=dp) FUNCTION rhoi(CurrentSolventMaterial,T0,p0,Temperature,Pressure,ConstVal) 2088 IMPLICIT NONE 2089 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2090 REAL(KIND=dp), INTENT(IN) :: T0,p0,Temperature,Pressure 2091 LOGICAL :: ConstVal 2092 !---------------------- 2093 REAL(KIND=dp) :: aux1, aux2 2094 !---------------------- 2095 rhoi = CurrentSolventMaterial % rhoi0 2096 IF (.NOT.ConstVal) THEN 2097 !aux1 = GeneralIntegral(Pressure,p0,p0,& 2098 ! CurrentSolventMaterial % ki0 ,& 2099 ! CurrentSolventMaterial % cki(0:5),& 2100 ! CurrentSolventMaterial % ckil) 2101 ! a shortcut, as only cki(0) = 1 2102 aux1 = (CurrentSolventMaterial % ki0 ) * (Pressure - p0) 2103 aux2 = GeneralIntegral(Temperature,T0,T0,& 2104 CurrentSolventMaterial % ai0,& 2105 CurrentSolventMaterial % aai(0:5),& 2106 CurrentSolventMaterial % aail) 2107 rhoi = rhoi * EXP(aux1 - aux2) 2108 END IF 2109 END FUNCTION rhoi 2110 !--------------------------------------------------------------------------------------------- 2111 REAL(KIND=dp) FUNCTION rhoiT(CurrentSolventMaterial,rhoi,T0,Temperature) 2112 IMPLICIT NONE 2113 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2114 REAL(KIND=dp), INTENT(IN) :: rhoi,T0,Temperature 2115 !------------------------ 2116 REAL(KIND=dp) :: alphaI 2117 !------------------------ 2118 alphaI = (CurrentSolventMaterial % ai0) *& 2119 GeneralPolynomial(Temperature,T0,T0,& 2120 CurrentSolventMaterial % aai(0:5),& 2121 CurrentSolventMaterial % aail) 2122 rhoiT = rhoi * alphaI 2123 END FUNCTION rhoiT 2124 !--------------------------------------------------------------------------------------------- 2125 REAL(KIND=dp) FUNCTION rhoiP(CurrentSolventMaterial,rhoi,p0,Pressure) 2126 IMPLICIT NONE 2127 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2128 REAL(KIND=dp), INTENT(IN) :: rhoi,p0,Pressure 2129 !------------------------ 2130 REAL(KIND=dp):: kappaI 2131 !------------------------ 2132 !kappaI = (CurrentSolventMaterial % ki0) *& 2133 ! GeneralPolynomial(Pressure,p0,p0,& 2134 ! CurrentSolventMaterial % cki(0:5),& 2135 ! CurrentSolventMaterial % ckil) 2136 ! a shortcut, as only cki(0) = 1 2137 kappaI = CurrentSolventMaterial % ki0 2138 rhoiP = rhoi * kappaI 2139 END FUNCTION rhoiP 2140 !--------------------------------------------------------------------------------------------- 2141 REAL (KIND=dp) FUNCTION rhoc(CurrentSoluteMaterial,T0,p0,Xi,Temperature,Pressure,Salinity,ConstVal) 2142 IMPLICIT NONE 2143 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2144 REAL(KIND=dp), INTENT(IN) :: T0,p0,Xi,Temperature,Pressure,Salinity 2145 LOGICAL :: ConstVal 2146 !---------------------- 2147 REAL(KIND=dp) :: aux1, aux2, aux3, xc 2148 !---------------------- 2149 rhoc = CurrentSoluteMaterial % rhoc0 2150 IF (.NOT.ConstVal) THEN 2151 xc = Salinity/Xi 2152 aux1 = GeneralIntegral(Pressure,p0,p0,& 2153 CurrentSoluteMaterial % kc0,& 2154 CurrentSoluteMaterial % ckc(0:5),& 2155 CurrentSoluteMaterial % ckcl) 2156 aux2 = GeneralIntegral(Temperature,T0,T0,& 2157 CurrentSoluteMaterial % ac0,& 2158 CurrentSoluteMaterial % aac(0:5),& 2159 CurrentSoluteMaterial % aacl) 2160 aux3 = GeneralIntegral(xc,0.0_dp,1.0_dp,& 2161 CurrentSoluteMaterial % zc0,& 2162 CurrentSoluteMaterial % bzc(0:5),& 2163 CurrentSoluteMaterial % bzcl) 2164 rhoc = rhoc * EXP(aux1 - aux2 + aux3) 2165 END IF 2166 END FUNCTION rhoc 2167 !--------------------------------------------------------------------------------------------- 2168 REAL(KIND=dp) FUNCTION rhocT(CurrentSoluteMaterial,rhoc,T0,Temperature,ConstVal) 2169 IMPLICIT NONE 2170 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2171 REAL(KIND=dp), INTENT(IN) :: rhoc,T0,Temperature 2172 LOGICAL :: ConstVal 2173 !------------------------- 2174 REAL(KIND=dp):: alphaC 2175 !------------------------- 2176 IF (ConstVal) THEN 2177 rhocT = 0.0_dp 2178 ELSE 2179 alphaC = (CurrentSoluteMaterial % ac0) * & 2180 GeneralPolynomial(Temperature,T0,T0,& 2181 CurrentSoluteMaterial % aac(0:5),& 2182 CurrentSoluteMaterial % aacl) 2183 rhocT = rhoc * alphaC 2184 END IF 2185 END FUNCTION rhocT 2186 !--------------------------------------------------------------------------------------------- 2187 REAL(KIND=dp) FUNCTION rhocP(CurrentSoluteMaterial,rhoc,ConstVal) 2188 IMPLICIT NONE 2189 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2190 REAL(KIND=dp), INTENT(IN) :: rhoc 2191 LOGICAL :: ConstVal 2192 !--------------- 2193 IF (ConstVal) THEN 2194 rhocP = 0.0_dp 2195 ELSE 2196 rhocP = rhoc * (CurrentSoluteMaterial % kc0) 2197 END IF 2198 END FUNCTION rhocP 2199 !--------------------------------------------------------------------------------------------- 2200 REAL(KIND=dp) FUNCTION rhocYc(CurrentSoluteMaterial,rhoc,Xi,Salinity,ConstVal) 2201 IMPLICIT NONE 2202 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2203 REAL(KIND=dp), INTENT(IN) :: rhoc, Xi, Salinity 2204 LOGICAL :: ConstVal 2205 !--------------- 2206 REAL(KIND=dp):: xc, zc 2207 !--------------- 2208 IF (ConstVal) THEN 2209 rhocYc = 0.0_dp 2210 ELSE 2211 xc = Salinity/Xi 2212 zc = (CurrentSoluteMaterial % zc0) * & 2213 GeneralPolynomial(xc,0.0_dp,1.0_dp,& 2214 CurrentSoluteMaterial % bzc(0:5),& 2215 CurrentSoluteMaterial % bzcl) 2216 rhocYc = rhoc * zc /Xi 2217 END IF 2218 END FUNCTION rhocYc 2219 !--------------------------------------------------------------------------------------------- 2220 REAL (KIND=dp) FUNCTION rhogw(rhow,rhoc,Xi,Salinity) 2221 IMPLICIT NONE 2222 REAL(KIND=dp), INTENT(IN) :: rhow,rhoc,Xi,Salinity 2223 !------------ 2224 REAL(KIND=dp) :: xc, LSalinity 2225 !------------ 2226!!$ IF (Salinity < 0.0_dp) THEN 2227!!$ CALL WARN("rhogw","Salinity smaller than 0") 2228!!$ LSalinity = 0.0_dp 2229!!$ ELSE IF (Salinity > 0.3_dp) THEN 2230!!$ CALL WARN("rhogw","Salinity larger than 0.3") 2231!!$ LSalinity = 0.3_dp 2232!!$ ELSE 2233!!$ LSalinity =Salinity 2234!!$ END IF 2235 xc = Salinity/Xi 2236 rhogw = rhow + xc*(rhoc - rhow) 2237 END FUNCTION rhogw 2238 !--------------------------------------------------------------------------------------------- 2239 REAL (KIND=dp) FUNCTION rhogwP(rhowp,rhocp,Xi,Salinity) 2240 IMPLICIT NONE 2241 REAL(KIND=dp), INTENT(IN) :: rhowp,rhocp,Xi,Salinity 2242 !------------ 2243 REAL(KIND=dp) :: xc 2244 !------------ 2245 xc = Salinity/Xi 2246 rhogwP = (1.0_dp - xc)*rhowP + xc*rhocP 2247 END FUNCTION rhogwP 2248 !--------------------------------------------------------------------------------------------- 2249 REAL (KIND=dp) FUNCTION rhogwT(rhowT,rhocT,Xi,Salinity) 2250 IMPLICIT NONE 2251 REAL(KIND=dp), INTENT(IN) :: rhowT,rhocT,Xi,Salinity 2252 !------------ 2253 REAL(KIND=dp) :: xc 2254 !------------ 2255 xc = Salinity/Xi 2256 rhogwT = (1.0_dp - xc)*rhowT + xc*rhocT 2257 END FUNCTION rhogwT 2258 !--------------------------------------------------------------------------------------------- 2259 REAL (KIND=dp) FUNCTION rhogwYc(rhow, rhoc, rhowYc,rhocYc,Xi,Salinity) 2260 IMPLICIT NONE 2261 REAL(KIND=dp), INTENT(IN) :: rhow, rhoc, rhowYc,rhocYc,Xi,Salinity 2262 !------------ 2263 REAL(KIND=dp) :: xc 2264 !------------ 2265 xc = Salinity/Xi 2266 rhogwYc = ((1.0_dp - xc)*rhowYc + xc*rhocYc + rhow + rhoc)/Xi 2267 END FUNCTION rhogwYc 2268 !--------------------------------------------------------------------------------------------- 2269 REAL (KIND=dp) FUNCTION cs(RockMaterialID,T0,Temperature,ConstVal) 2270 IMPLICIT NONE 2271 INTEGER, INTENT(IN) :: RockMaterialID 2272 REAL(KIND=dp), INTENT(IN) :: T0,Temperature 2273 LOGICAL :: ConstVal 2274 !---------------------- 2275 REAL(KIND=dp) :: aux 2276 !---------------------- 2277 cs = GlobalRockMaterial % cs0(RockMaterialID) 2278 !PRINT *,"cs:", Temperature,T0,GlobalRockMaterial % cs0(RockMaterialID),& 2279 ! GlobalRockMaterial % acs(0:5,RockMaterialID),GlobalRockMaterial % acsl(RockMaterialID) 2280 IF (.NOT.ConstVal) & 2281 cs = cs * GeneralPolynomial(Temperature,T0,T0,& 2282 GlobalRockMaterial % acs(0:5,RockMaterialID),& 2283 GlobalRockMaterial % acsl(RockMaterialID)) 2284 END FUNCTION cs 2285 !--------------------------------------------------------------------------------------------- 2286 REAL (KIND=dp) FUNCTION cw(CurrentSolventMaterial,T0,Xi,Temperature,Salinity,ConstVal) 2287 IMPLICIT NONE 2288 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2289 REAL(KIND=dp), INTENT(IN) :: T0,Xi,Temperature,Salinity 2290 LOGICAL :: ConstVal 2291 !---------------------- 2292 REAL(KIND=dp) :: aux1, aux2, watercont 2293 !---------------------- 2294 cw = CurrentSolventMaterial % cw0 2295 IF (.NOT.ConstVal) THEN 2296 watercont = MAX(1.0_dp - Salinity/Xi,0.0_dp) 2297 aux1 = GeneralPolynomial(Temperature,T0,T0,& 2298 CurrentSolventMaterial % acw(0:5),& 2299 CurrentSolventMaterial % acwl) 2300 aux2 = GeneralPolynomial(watercont,1.0_dp,1.0_dp,& 2301 CurrentSolventMaterial % bcw(0:5),& 2302 CurrentSolventMaterial % bcwl) 2303 cw = cw * aux1 * aux2 2304 END IF 2305 END FUNCTION cw 2306 !--------------------------------------------------------------------------------------------- 2307 REAL (KIND=dp) FUNCTION ci(CurrentSolventMaterial,& 2308 T0,Temperature,ConstVal) 2309 IMPLICIT NONE 2310 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2311 REAL(KIND=dp), INTENT(IN) :: T0,Temperature 2312 REAL(KIND=dp) :: ci0 2313 REAL(KIND=dp), DIMENSION(0:5) :: aci 2314 INTEGER :: acil 2315 LOGICAL :: ConstVal 2316 !---------------------- 2317 ci = CurrentSolventMaterial % ci0 2318 IF (.NOT.ConstVal) THEN 2319 !PRINT *, "ci:",ci0,aci(0:5),acil 2320 ci = ci *& 2321 GeneralPolynomial(Temperature,T0,T0,& 2322 CurrentSolventMaterial % aci(0:5),& 2323 CurrentSolventMaterial % acil) 2324 END IF 2325 END FUNCTION ci 2326 !--------------------------------------------------------------------------------------------- 2327 REAL (KIND=dp) FUNCTION cc(CurrentSoluteMaterial,& 2328 T0,Temperature,Salinity,ConstVal) 2329 IMPLICIT NONE 2330 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2331 REAL(KIND=dp), INTENT(IN) :: T0,Temperature,Salinity 2332 LOGICAL :: ConstVal 2333 !---------------------- 2334 REAL(KIND=dp) :: aux1, aux2 2335 !---------------------- 2336 cc = CurrentSoluteMaterial % cc0 2337 IF (.NOT.ConstVal) THEN 2338 aux1 = GeneralPolynomial(Temperature,T0,T0,& 2339 CurrentSoluteMaterial % acc(0:5),& 2340 CurrentSoluteMaterial % accl) 2341 aux2 = GeneralPolynomial(Salinity,0.0_dp,1.0_dp,& 2342 CurrentSoluteMaterial % bcc(0:5),& 2343 CurrentSoluteMaterial % bccl) 2344 cc = cc*aux1*aux2 2345 END IF 2346 END FUNCTION cc 2347 !--------------------------------------------------------------------------------------------- 2348 ! latent heat of water 2349 REAL (KIND=dp) FUNCTION hw(CurrentSolventMaterial,& 2350 T0,Xi,Temperature,Salinity,ConstVal) 2351 IMPLICIT NONE 2352 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2353 REAL(KIND=dp), INTENT(IN) :: T0,Xi,Temperature,Salinity 2354 LOGICAL :: ConstVal 2355 !---------------------- 2356 REAL(KIND=dp) :: aux1, aux2, watercont 2357 !---------------------- 2358 hw = CurrentSolventMaterial % hw0 2359 IF (.NOT.ConstVal) THEN 2360 watercont = MAX(1.0_dp - Salinity/Xi,0.0_dp) 2361 aux1 = GeneralPolynomial(watercont,1.0_dp,1.0_dp,& 2362 CurrentSolventMaterial % bcw(0:5),& 2363 CurrentSolventMaterial % bcwl) 2364 aux2 = GeneralIntegral(Temperature,T0,T0,& 2365 CurrentSolventMaterial % cw0,& 2366 CurrentSolventMaterial % acw(0:5),& 2367 CurrentSolventMaterial % acwl) 2368 hw = hw + aux1 * aux2 2369 END IF 2370 END FUNCTION hw 2371 !--------------------------------------------------------------------------------------------- 2372 ! latent heat of ice 2373 REAL (KIND=dp) FUNCTION hi(CurrentSolventMaterial,& 2374 T0,Temperature,ConstVal) 2375 IMPLICIT NONE 2376 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2377 REAL(KIND=dp), INTENT(IN) :: T0,Temperature 2378 LOGICAL :: ConstVal 2379 !---------------------- 2380 REAL(KIND=DP) :: hi0,ci0 2381 REAL(KIND=DP), DIMENSION(0:5) :: aci 2382 INTEGER :: acil 2383 LOGICAL :: FirstTime = .TRUE. 2384 SAVE FirstTime,hi0,ci0,aci,acil 2385 !---------------------- 2386 hi0 = CurrentSolventMaterial % hi0 2387 IF (ConstVal) THEN 2388 hi = hi0 2389 ELSE 2390 IF (FirstTime) THEN 2391 ci0 = CurrentSolventMaterial % ci0 2392 aci(0:5) = CurrentSolventMaterial % aci(0:5) 2393 acil = CurrentSolventMaterial % acil 2394 FirstTime = .FALSE. 2395 END IF 2396 hi = hi0 + GeneralIntegral(Temperature,T0,T0,ci0,aci,acil) 2397 END IF 2398 END FUNCTION hi 2399 !--------------------------------------------------------------------------------------------- 2400 ! General constituent thermal conductivity: kalpha0th and balpha have to be directly transferred 2401 FUNCTION GetKAlphaTh(kalpha0th,balpha,T0,Temperature)RESULT(kalphath) 2402 IMPLICIT NONE 2403 REAL(KIND=dp), INTENT(IN) :: kalpha0th,balpha,T0,Temperature 2404 REAL(KIND=dp) :: kalphath 2405 !------------------------- 2406 kalphath = kalpha0th/( 1.0_dp + balpha*(Temperature - T0)/T0) 2407 !kalphath = kalpha0th 2408 END FUNCTION GetKAlphaTh 2409 !--------------------------------------------------------------------------------------------- 2410 FUNCTION GetCGTT(Xi,XiT,rhos,rhow,rhoi,rhoc,cw,ci,cs,cc,hi,hw,& 2411 Porosity,Salinity)RESULT(CGTT)! All state variables or derived values 2412 IMPLICIT NONE 2413 REAL(KIND=dp), INTENT(IN) :: Xi,XiT,rhos,rhow,rhoi,rhoc,cw,ci,cs,cc,& 2414 hi,hw,Porosity,Salinity 2415 REAL(KIND=dp) :: CGTT 2416 !------------------------- 2417 REAL(KIND=dp) :: xc 2418 !------------------------- 2419 xc = Salinity/Xi 2420 CGTT = (1.0_dp - Porosity)*rhos*cs & 2421 + (Xi - Salinity) * Porosity * rhow * cw & ! mind xc * Xi = Salinity 2422 + Salinity * Porosity * rhoc * cc & ! mind xc * Xi = Salinity 2423 + (1.0_dp - Xi)*Porosity*rhoi*ci & 2424 + rhoi*(hw - hi)*Porosity*XiT 2425 END FUNCTION GetCGTT 2426 !--------------------------------------------------------------------------------------------- 2427 FUNCTION GetCGTp(rhoi,hi,hw,XiP,Porosity)RESULT(CGTp)! All state variables or derived values 2428 IMPLICIT NONE 2429 REAL(KIND=dp), INTENT(IN) :: rhoi,hi,hw,XiP,Porosity 2430 REAL(KIND=dp) :: CGTp 2431 !------------------------- 2432 CGTp = Porosity*rhoi*(hw - hi)*XiP 2433 END FUNCTION GetCGTp 2434 !--------------------------------------------------------------------------------------------- 2435 FUNCTION GetCGTyc(rhoi,hi,hw,XiYc,Porosity)RESULT(CGTyc)! All state variables or derived values 2436 IMPLICIT NONE 2437 REAL(KIND=dp), INTENT(IN) :: rhoi,hi,hw,XiYc,Porosity 2438 REAL(KIND=dp) :: CGTyc 2439 !------------------------- 2440 CGTyc = Porosity*rhoi*(hw - hi)*XiYc 2441 END FUNCTION GetCGTyc 2442 !--------------------------------------------------------------------------------------------- 2443 ! functions specific to groundwater flow 2444 !--------------------------------------------------------------------------------------------- 2445 FUNCTION GetJgwD(Kgwpp,KgwpT,Kgw,gradp,gradT,Gravity,rhogw,DIM,CryogenicSuction) RESULT(JgwD) 2446 IMPLICIT NONE 2447 REAL (KIND=dp), INTENT(IN) :: Kgwpp(3,3),KgwpT(3,3),Kgw(3,3),gradp(3),gradT(3),Gravity(3),& 2448 rhogw 2449 REAL (KIND=dp) :: JgwD(3) 2450 LOGICAL, INTENT(IN):: CryogenicSuction 2451 !------------------------- 2452 INTEGER, INTENT(IN) :: DIM 2453 INTEGER :: i 2454 REAL (KIND=dp) :: fluxp(3),fluxT(3),fluxg(3) 2455 !------------------------- 2456 fluxT(1:DIM) = 0.0_dp 2457 JgwD = 0.0_dp 2458 DO i=1,DIM 2459 fluxp(i) = -1.0_dp * SUM(Kgwpp(i,1:DIM)*gradp(1:DIM)) 2460 IF (CryogenicSuction) & 2461 fluxT(i) = -1.0_dp * SUM(KgwpT(i,1:DIM)*gradT(1:DIM)) 2462 fluxg(i) = rhogw * SUM(Kgw(i,1:DIM)*Gravity(1:DIM)) 2463 END DO 2464 JgwD(1:DIM) = fluxp(1:DIM) + fluxT(1:DIM) + fluxg(1:DIM) 2465 END FUNCTION GetJgwD 2466 !--------------------------------------------------------------------------------------------- 2467 FUNCTION GetKGTT(ksth,kwth,kith,kcth,Xi,& 2468 Salinity,Porosity,meanfactor)RESULT(KGTT) ! All state variables or derived values 2469 IMPLICIT NONE 2470 REAL(KIND=dp), INTENT(IN) :: ksth,kwth,kith,kcth,Xi,& 2471 Salinity,Porosity,meanfactor 2472 REAL(KIND=dp) :: KGTT(3,3) 2473 !------------------------- 2474 REAL(KIND=dp) :: KGaTT, KghTT, unittensor(3,3),xc 2475 !------------------------- 2476 xc = Salinity/Xi 2477 unittensor=RESHAPE([1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0], SHAPE(unittensor)) 2478 KGhTT = 1.0_dp/((1.0_dp - Porosity)/ksth + (1.0_dp - xc)*Xi*Porosity/kwth & 2479 + xc*Porosity/kcth + (1.0_dp - Xi)*Porosity/kith) 2480 KGaTT = (1.0_dp - Porosity)*ksth + (1.0_dp - xc)*Xi*Porosity*kwth & 2481 + xc*Porosity*kcth + (1.0_dp - Xi)*Porosity*kith 2482 KGTT = unittensor*((1.0_dp - meanfactor)*KGhTT + meanfactor * KGaTT) 2483 END FUNCTION GetKGTT 2484 !--------------------------------------------------------------------------------------------- 2485 FUNCTION GetDtd(RockMaterialID,Xi,Porosity,JgwD)RESULT(Dtd) 2486 IMPLICIT NONE 2487 REAL(KIND=dp), INTENT(IN) :: Xi,Porosity,JgwD(3) 2488 REAL(KIND=dp) :: Dtd(3,3) 2489 INTEGER, INTENT(IN) :: RockMaterialID 2490 !------------------------- 2491 REAL(KIND=dp) :: unittensor(3,3),absJgwD,alphaL,alphaT 2492 INTEGER :: I,J 2493 !------------------------- 2494 unittensor=RESHAPE([1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0], SHAPE(unittensor)) 2495 absJgwD = SQRT(SUM(JgwD(1:3)*JgwD(1:3))) 2496 IF(absJgwD > 0.0_dp) THEN 2497 alphaL = GlobalRockMaterial % alphaL(RockMaterialID) 2498 alphaT = GlobalRockMaterial % alphaT(RockMaterialID) 2499 DO I=1,3 2500 DO J=1,3 2501 Dtd(I,J) = alphaT*absJgwD*unittensor(I,J) + (alphaL - alphaT)*JgwD(I)*JgwD(J)/absJgwD 2502 END DO 2503 END DO 2504 ELSE 2505 Dtd = 0.0_dp 2506 END IF 2507 END FUNCTION GetDtd 2508 !--------------------------------------------------------------------------------------------- 2509 FUNCTION GetCgwTT(rhow,rhoc,cw,cc,Xi,Salinity)RESULT(CgwTT)! All state variables or derived values 2510 IMPLICIT NONE 2511 REAL(KIND=dp), INTENT(IN) :: rhow,rhoc,cw,cc,Xi,Salinity 2512 REAL(KIND=dp) :: CgwTT 2513 !------------------------- 2514 REAL(KIND=dp) :: xc 2515 !------------------------- 2516 xc = Salinity/Xi 2517 CgwTT = (1.0_dp - xc)*rhow*cw + xc*rhoc*cc 2518 END FUNCTION GetCgwTT 2519 !--------------------------------------------------------------------------------------------- 2520 FUNCTION GetCgwpp(rhogw,rhoi,rhogwp,rhoip,rhosp,& 2521 kappaG,Xi,Xip,& 2522 RockMaterialID,Porosity)RESULT(Cgwpp) 2523 IMPLICIT NONE 2524 REAL(KIND=dp), INTENT(IN) :: rhogw,rhoi,rhogwp,rhoip,rhosp,kappaG,Xi,Xip,Porosity 2525 INTEGER, INTENT(IN) :: RockMaterialID 2526 REAL(KIND=dp) :: Cgwpp 2527 !------------------------- 2528 Cgwpp = Porosity * ( (rhogw - rhoi) * Xip + Xi * rhogwp + (1.0_dp - Xi)*rhoip ) & 2529 + (Xi * rhogw + (1.0_dp - Xi)*rhoi) * ( (1.0_dp - Porosity) * rhosp + kappaG ) 2530 END FUNCTION GetCgwpp 2531 !--------------------------------------------------------------------------------------------- 2532 FUNCTION GetCgwpT(rhogw,rhoi,rhogwT,rhoiT,rhosT,Xi,XiT,Porosity)RESULT(CgwpT) 2533 IMPLICIT NONE 2534 REAL(KIND=dp), INTENT(IN) :: rhogw,rhoi,rhogwT,rhoiT,rhosT,Xi,XiT,Porosity 2535 REAL(KIND=dp) :: CgwpT 2536 !------------------------- 2537 CgwpT = Porosity * ( (rhogw - rhoi) * XiT + Xi * rhogwT + (1.0_dp - Xi)*rhoiT ) & 2538 + (Xi * rhogw + (1.0_dp - Xi)*rhoi) *(1.0_dp - Porosity) * rhosT 2539 END FUNCTION GetCgwpT 2540 !--------------------------------------------------------------------------------------------- 2541 FUNCTION GetCgwpYc(rhogw,rhoi,rhogwYc,Xi,XiYc,Porosity)RESULT(CgwpYc) 2542 IMPLICIT NONE 2543 REAL(KIND=dp), INTENT(IN) :: rhogw,rhoi,rhogwYc,Xi,XiYc,Porosity 2544 REAL(KIND=dp) :: CgwpYc 2545 !------------------------- 2546 CgwpYc = Porosity * ( (rhogw - rhoi) * XiYc + Xi * rhogwYc ) 2547 END FUNCTION GetCgwpYc 2548 !--------------------------------------------------------------------------------------------- 2549 FUNCTION GetCgwpI1(rhogw,rhoi,Xi,kappaG,RockMaterialID)RESULT(CgwpI1) 2550 IMPLICIT NONE 2551 REAL(KIND=dp), INTENT(IN) :: rhogw,rhoi,Xi,kappaG 2552 INTEGER, INTENT(IN) :: RockMaterialID 2553 REAL(KIND=dp) :: CgwpI1 2554 !------------------------- 2555 REAL(KIND=dp) :: kappas 2556 !------------------------- 2557 kappas = GlobalRockMaterial % ks0(RockMaterialID) 2558 CgwpI1 = (Xi * rhogw + (1.0_dp - Xi) * rhoi)*(kappaG)/3.0_dp 2559 END FUNCTION GetCgwpI1 2560 !--------------------------------------------------------------------------------------------- 2561 REAL (KIND=dp) FUNCTION mugw(CurrentSolventMaterial,CurrentSoluteMaterial,& 2562 Xi,T0,Salinity,Temperature,ConstVal) 2563 IMPLICIT NONE 2564 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2565 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2566 REAL(KIND=dp), INTENT(IN) :: Xi,T0,Salinity,Temperature 2567 LOGICAL :: ConstVal 2568 !------------------------- 2569 REAL(KIND=dp) :: nu1, nu2, xc 2570 !------------------------- 2571 IF (.NOT.ConstVal) THEN 2572 xc = MAX(Salinity/Xi,0.0_dp) 2573 xc = MIN(xc,0.2_dp) 2574 nu1 = (CurrentSolventMaterial % nu10) *& 2575 GeneralPolynomial(Temperature,T0,T0,& 2576 CurrentSolventMaterial % anw(0:5),& 2577 CurrentSolventMaterial % anwl) 2578 ! PRINT *,"mugw:anw,anwl,nu1", anw(0:5),anwl, nu1 2579 ! PRINT *,"mugw:bnc", bnc(0:5),bncl 2580 ! PRINT *,"bnc=",bnc(0:5),"bncl=",bncl,"xc=",xc 2581 ! STOP 2582 nu2 = (CurrentSoluteMaterial % nu20) *& 2583 GeneralPolynomial(xc,0.0_dp,1.0_dp,& 2584 CurrentSoluteMaterial % bnc(0:5),& 2585 CurrentSoluteMaterial % bncl) 2586 mugw = CurrentSolventMaterial % muw0 * & 2587 EXP(nu1 * (Temperature - T0) + nu2 * (xc - 0.0)) 2588 IF ((mugw .NE. mugw) .OR. (mugw > HUGE(mugw))) THEN 2589 WRITE (Message,*) 'invalid value:' , mugw,& 2590 'Input values: muw0/nu1/nu2/T/T0/xc/Salinity/Xi: ',& 2591 CurrentSolventMaterial % muw0, nu1,nu2,Temperature, T0, xc, Salinity, Xi 2592 ! PRINT *, nu1 * (Temperature - T0), nu2 * (xc - 0.0) 2593 CALL FATAL("PermafrostMaterials(mugw)",Message) 2594 END IF 2595 ELSE 2596 mugw = CurrentSolventMaterial % muw0 2597 END IF 2598 END FUNCTION mugw 2599 !--------------------------------------------------------------------------------------------- 2600 FUNCTION GetKgw(RockMaterialID,CurrentSolventMaterial,& 2601 mugw,Xi,MinKgw)RESULT(Kgw) 2602 IMPLICIT NONE 2603 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2604 INTEGER, INTENT(IN) :: RockMaterialID 2605 REAL(KIND=dp), INTENT(IN) :: Xi,MinKgw,mugw 2606 REAL(KIND=dp) :: Kgw(3,3) 2607 !-------------------------- 2608 REAL(KIND=dp) :: muw0,rhow0,qexp,Kgwh0(3,3),factor 2609 REAL(KIND=dp), PARAMETER :: gval=9.81_dp !hard coded, so match Kgwh0 with this value 2610 INTEGER :: I, J 2611 !------------------------- 2612 IF (mugw <= 0.0_dp) & 2613 CALL FATAL("Permafrost(GetKgw)","Unphysical viscosity detected") 2614 muw0 = CurrentSolventMaterial % muw0 2615 rhow0 = CurrentSolventMaterial % rhow0 2616 qexp = GlobalRockMaterial % qexp(RockMaterialID) 2617 Kgwh0(1:3,1:3) = GlobalRockMaterial % Kgwh0(1:3,1:3,RockMaterialID) ! hydro-conductivity 2618 ! transformation factor from hydr. conductivity to permeability hydr. conductivity tensor 2619 factor = (muw0/mugw)*(Xi**qexp)/(rhow0*gval) 2620 ! this is OK 2621 !PRINT *,"Kgw:",muw0,mugw,rhow0,Kgwh0,Xi,factor 2622 Kgw = 0.0_dp 2623 DO I=1,3 2624 DO J=1,3 2625 Kgw(i,j) = Kgwh0(i,j)*factor 2626 END DO 2627 END DO 2628 DO I=1,3 2629 Kgw(i,i) = MAX(Kgw(i,i),MinKgw) 2630 END DO 2631 END FUNCTION GetKgw 2632 !--------------------------------------------------------------------------------------------- 2633 FUNCTION GetKgwpT(fw,XiT,Kgw)RESULT(KgwpT) ! All state variables or derived values 2634 IMPLICIT NONE 2635 REAL(KIND=dp), INTENT(IN) :: fw,XiT,Kgw(3,3) 2636 REAL(KIND=dp) :: KgwpT(3,3) 2637 !------------------------- 2638 KgwpT(1:3,1:3) = fw*XiT*Kgw(1:3,1:3) 2639 END FUNCTION GetKgwpT 2640 !--------------------------------------------------------------------------------------------- 2641 FUNCTION GetKgwpp(fw,XiP,Kgw)RESULT(Kgwpp)! All state variables or derived values 2642 IMPLICIT NONE 2643 REAL(KIND=dp), INTENT(IN) :: fw,XiP,Kgw(3,3) 2644 REAL(KIND=dp) :: Kgwpp(3,3) 2645 !------------------------- 2646 Kgwpp(1:3,1:3) = (1.0_dp + fw*XiP)*Kgw(1:3,1:3) 2647 END FUNCTION GetKgwpp 2648 !--------------------------------------------------------------------------------------------- 2649 ! functions specific to solute transport 2650 !--------------------------------------------------------------------------------------------- 2651 FUNCTION GetKc(RockMaterialID,Dm,Xi,JgwD,Porosity)RESULT(Kc) 2652 IMPLICIT NONE 2653 REAL(KIND=dp), INTENT(IN) :: Dm,Xi,JgwD(3),Porosity 2654 INTEGER, INTENT(IN) :: RockMaterialID 2655 REAL(KIND=dp) :: alphaL,alphaT,Kc(3,3), unittensor(3,3), aux, eL(3),absJgwD 2656 INTEGER :: I,J 2657 !------------------------- 2658 unittensor=RESHAPE([1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0], SHAPE(unittensor)) 2659 IF (Porosity <= 0.0_dp) & 2660 CALL FATAL("GetKc","Negative/Zero Porosity detected") 2661 IF (Xi <= 0.0_dp) & 2662 CALL FATAL("GetKc","Negative/Zero water content detected") 2663 Kc = Dm * unittensor 2664 absJgwD = SQRT(SUM(JgwD(1:3) * JgwD(1:3))) 2665 IF (absJgwD > 0.0_dp) THEN 2666 alphaL = GlobalRockMaterial % alphaL(RockMaterialID) 2667 alphaT = GlobalRockMaterial % alphaT(RockMaterialID) 2668 eL = JgwD/absJgwD 2669 aux = absJgwD/(Porosity * Xi) 2670 DO I=1,3 2671 DO J=1,3 2672 Kc(I,J) = Kc(I,J) & 2673 + aux*((alphaL - alphaT)*eL(I)*eL(J) + alphaT * unittensor(I,J)) 2674 END DO 2675 END DO 2676 END IF 2677 END FUNCTION GetKc 2678 !--------------------------------------------------------------------------------------------- 2679 FUNCTION GetConstKc(DispersionCoefficient)RESULT(Kc) 2680 IMPLICIT NONE 2681 REAL(KIND=dp), INTENT(IN) :: DispersionCoefficient 2682 REAL(KIND=dp) :: Kc(3,3) 2683 !------------------------- 2684 REAL(KIND=dp) :: unittensor(3,3) 2685 unittensor=RESHAPE([1.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,1.0], SHAPE(unittensor)) 2686 Kc = DispersionCoefficient * unittensor 2687 END FUNCTION GetConstKc 2688 !--------------------------------------------------------------------------------------------- 2689 REAL(KIND=dp) FUNCTION Dm(CurrentSoluteMaterial,N0,GasConstant,rhoc,mugw,Temperature) 2690 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2691 REAL(KIND=dp), INTENT(IN) :: N0,GasConstant,rhoc,mugw,Temperature 2692 !------------------------- 2693 REAL(KIND=dp), PARAMETER :: bconst = 3.0 * PI 2694 REAL(KIND=dp) :: Mc, lcbar 2695 !------------------------- 2696 Mc = CurrentSoluteMaterial % Mc 2697 lcbar = (Mc/(rhoc * N0))**(1.0_dp/3.0_dp) 2698 Dm = GasConstant * Temperature / (bconst * mugw * lcbar * N0) 2699 END FUNCTION Dm 2700 !--------------------------------------------------------------------------------------------- 2701 FUNCTION GetR(CurrentSoluteMaterial,CurrentSolventMaterial,GasConstant,rhow,rhoc,Xi,Temperature,Salinity) RESULT(r12) 2702 TYPE(SoluteMaterial_t), POINTER :: CurrentSoluteMaterial 2703 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2704 REAL(KIND=dp), INTENT(IN) :: GasConstant,rhow,rhoc,Xi,Salinity,Temperature 2705 REAL(KIND=dp) :: r12(2) 2706 REAL(KIND=dp) :: d1, d2, Mc, Mw, aux, epsilonc 2707 !------------------------- 2708 d1 = CurrentSoluteMaterial % d1 2709 d2 = CurrentSoluteMaterial % d2 2710 Mc = CurrentSoluteMaterial % Mc 2711 Mw = CurrentSolventMaterial % Mw 2712 epsilonc = (Mc/Mw)*(rhow/rhoc) 2713 aux = Salinity/(Xi - Salinity) 2714 r12(1) = (1/epsilonc)*Mc*(1.0_dp - Salinity/Xi)/(rhoc * GasConstant * Temperature) 2715 r12(2) = epsilonc * ( d1 + (d1 + d2)*aux + d2*aux*aux ) 2716 END FUNCTION GetR 2717 !--------------------------------------------------------------------------------------------- 2718 FUNCTION GetKcYcYc(Kc,r12) RESULT(KcYcYc)! All state variables or derived values 2719 IMPLICIT NONE 2720 REAL(KIND=dp), INTENT(IN) :: Kc(3,3),r12(2) 2721 REAL(KIND=dp) :: KcYcYc(3,3) 2722 !------------------------- 2723 KcYcYc(1:3,1:3) = r12(2) * Kc(1:3,1:3) 2724 END FUNCTION GetKcYcYc 2725 !--------------------------------------------------------------------------------------------- 2726 FUNCTION GetFc(rhoc,rhow,Gravity,r12,XiT,XiP,Xi,gradP,gradT) RESULT(fc)! All state variables or derived values 2727 IMPLICIT NONE 2728 REAL(KIND=dp), INTENT(IN) :: rhoc,rhow,Gravity(3),r12(2),XiT,XiP,Xi,gradP(3),gradT(3) 2729 REAL(KIND=dp) :: fc(3) 2730 !------------------------- 2731 fc(1:3) = r12(1)*(rhoc - rhow)*Gravity(1:3) + r12(2)*(XiT*gradT(1:3) + XiP*gradP(1:3))/Xi 2732 END FUNCTION GetFc 2733 !--------------------------------------------------------------------------------------------- 2734 FUNCTION GetJcF(KcYcYc,Kc,fc,GradSalinity,Salinity) RESULT(JcF) 2735 IMPLICIT NONE 2736 REAL(KIND=dp), INTENT(IN) :: KcYcYc(3,3),Kc(3,3),fc(3), GradSalinity(3),Salinity 2737 REAL(KIND=dp) :: JcF(3) 2738 !----------- 2739 !REAL(KIND=dp) :: diffFlux(3), extForce(3) 2740 INTEGER :: i 2741 DO i=1,3 2742 !diffFlux(i) = SUM(KcYcYc(i,1:3) * GradSalinity(1:3)) 2743 !extForce(i) = SUM(Kc(i,1:3) * fc(1:3)) * Salinity 2744 JcF(i) = -SUM(KcYcYc(i,1:3)*GradSalinity(1:3)) + SUM(Kc(i,1:3) * fc(1:3)) * Salinity 2745 !SUM(KgwAtIP(i,1:DIM)*Gravity(1:DIM)) 2746 END DO 2747 END FUNCTION GetJcF 2748 !--------------------------------------------------------------------------------------------- 2749 REAL(KIND=dp) FUNCTION CcYcT(rhocT,Porosity,Salinity)! All state variables or derived values 2750 IMPLICIT NONE 2751 REAL(KIND=dp), INTENT(IN) :: rhocT,Porosity, Salinity 2752 CcYcT = Porosity*Salinity*rhocT 2753 END FUNCTION CcYcT 2754 !--------------------------------------------------------------------------------------------- 2755 REAL(KIND=dp) FUNCTION CcYcP(rhocP,Porosity, Salinity)! All state variables or derived values 2756 IMPLICIT NONE 2757 REAL(KIND=dp), INTENT(IN) :: rhocP,Porosity, Salinity 2758 !------------------------- 2759 CcYcP = Porosity*Salinity*rhocp 2760 END FUNCTION CcYcP 2761 !--------------------------------------------------------------------------------------------- 2762 REAL(KIND=dp) FUNCTION CcYcYc(rhoc,rhocYc,Porosity, Salinity)! All state variables or derived values 2763 IMPLICIT NONE 2764 REAL(KIND=dp), INTENT(IN) :: rhoc,rhocYc,Porosity, Salinity 2765 !------------------------- 2766 CcYcYc = Porosity*(rhoc + Salinity*rhocYc) 2767 END FUNCTION CcYcYc 2768 !--------------------------------------------------------------------------------------------- 2769 REAL(Kind=dp) FUNCTION RadiogenicHeatProduction(RockMaterialID,Depth,RefDepth) 2770 IMPLICIT NONE 2771 REAL(KIND=dp), INTENT(IN) :: Depth,RefDepth 2772 INTEGER, INTENT(IN) :: RockMaterialID 2773 !--------- 2774 RadiogenicHeatProduction = GlobalRockMaterial % RadGen(RockMaterialID) & 2775 * EXP(-Depth/RefDepth) 2776 END FUNCTION RadiogenicHeatProduction 2777 !--------------------------------------------------------------------------------------------- 2778 ! functions specific to ground deformation 2779 !--------------------------------------------------------------------------------------------- 2780 REAL(Kind=dp) FUNCTION EG(CurrentSolventMaterial,RockMaterialID,Xi,Porosity) 2781 IMPLICIT NONE 2782 REAL(KIND=dp), INTENT(IN) :: Xi,Porosity 2783 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2784 INTEGER, INTENT(IN) :: RockMaterialID 2785 !PRINT *,"EG", Porosity, Xi, RockMaterialID, GlobalRockMaterial % Es0(RockMaterialID) 2786 EG = (1.0_dp - Porosity)*(GlobalRockMaterial % Es0(RockMaterialID))& 2787 /(1.0_dp - (GlobalRockMaterial % eta0(RockMaterialID))) & 2788 + Porosity * (1.0_dp - Xi) * (CurrentSolventMaterial % Ei0) 2789 END FUNCTION EG 2790 !--------------------------------------------------------------------------------------------- 2791 REAL(Kind=dp) FUNCTION nuG(CurrentSolventMaterial,RockMaterialID,Xi,Porosity) 2792 IMPLICIT NONE 2793 REAL(KIND=dp), INTENT(IN) :: Xi,Porosity 2794 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2795 INTEGER, INTENT(IN) :: RockMaterialID 2796 !--------- 2797 nuG = (1.0_dp - Porosity)*(GlobalRockMaterial % nuS0(RockMaterialID))& 2798 + Porosity * (1.0_dp - Xi) * (CurrentSolventMaterial % nui0) 2799 END FUNCTION nuG 2800 !--------------------------------------------------------------------------------------------- 2801 REAL(Kind=dp) FUNCTION betaG(CurrentSolventMaterial,RockMaterialID,Xi,Porosity) 2802 IMPLICIT NONE 2803 REAL(KIND=dp), INTENT(IN) :: Xi,Porosity 2804 TYPE(SolventMaterial_t), POINTER :: CurrentSolventMaterial 2805 INTEGER, INTENT(IN) :: RockMaterialID 2806 !--------- 2807 betaG = (1.0_dp - Porosity)*(GlobalRockMaterial % betas(RockMaterialID)& 2808 + Porosity * (1.0_dp - Xi) * (CurrentSolventMaterial % betai)) 2809 END FUNCTION BetaG 2810 !--------------------------------------------------------------------------------------------- 2811 REAL(Kind=dp) FUNCTION rhoG(rhos,rhogw,rhoi,Porosity,Salinity,Xi) 2812 IMPLICIT NONE 2813 REAL(KIND=dp), INTENT(IN) :: rhos,rhogw,rhoi,Porosity,Salinity,Xi 2814 !--------- 2815 rhoG = (1.0_dp - Porosity)*rhos + Porosity*Xi*(1.0_dp - Salinity)*rhogw & 2816 + Porosity*(1.0_dp - Xi)*rhoi 2817 END FUNCTION rhoG 2818 2819 !--------------------------------------------------------------------------------------------- 2820 FUNCTION KGuu(EG,nuG,DIM) RESULT(OutKGuu) 2821 REAL(KIND=dp), INTENT(IN) :: EG,nuG 2822 REAL(KIND=dp) OutKGuu(6,6) 2823 INTEGER, INTENT(IN) :: DIM 2824 !---------- 2825 INTEGER :: I,J 2826 !--------- 2827 OutKGuu = 0.0_dp 2828 DO I=1,DIM 2829 OutKGuu(I,I) = 1.0_dp - nuG 2830 OutKGuu(DIM+I,DIM+I) = 0.5_dp - nuG 2831 DO J=1,DIM 2832 IF (J /= I) OutKGuu(I,J) = nuG 2833 END DO 2834 END DO 2835 END FUNCTION KGuu 2836 !--------------------------------------------------------------------------------------------- 2837 REAL(Kind=dp) FUNCTION kappaG(EG,nuG) ! needed directly in Darcy Model 2838 IMPLICIT NONE 2839 REAL(KIND=dp), INTENT(IN) :: EG,nuG 2840 !--------- 2841 kappaG = (3.0_dp*(1.0_dp - 2.0_dp * nuG))/EG 2842 END FUNCTION KappaG 2843 !--------------------------------------------------------------------------------------------- 2844END MODULE PermafrostMaterials 2845