| 1 | subroutine soilthprop(porosity,fill,rhocobs,tiobs,layertype, & |
|---|
| 2 | & newrhoc,newti,icefrac) |
|---|
| 3 | !*********************************************************************** |
|---|
| 4 | ! soilthprop: assign thermal properties of icy soil or dirty ice |
|---|
| 5 | ! |
|---|
| 6 | ! porositiy = void space / total volume |
|---|
| 7 | ! rhof = density of free ice in space not occupied by regolith [kg/m^3] |
|---|
| 8 | ! fill = rhof/icedensity <=1 (only relevant for layertype 1) |
|---|
| 9 | ! rhocobs = heat capacity per volume of dry regolith [J/m^3] |
|---|
| 10 | ! tiobs = thermal inertia of dry regolith [SI-units] |
|---|
| 11 | ! layertype: 1=interstitial ice, 2=pure ice or ice with dirt |
|---|
| 12 | ! 3=pure ice, 4=ice-cemented soil, 5=custom values |
|---|
| 13 | ! icefrac: fraction of ice in icelayer (only relevant for layertype 2) |
|---|
| 14 | ! output are newti and newrhoc |
|---|
| 15 | !*********************************************************************** |
|---|
| 16 | implicit none |
|---|
| 17 | integer, intent(IN) :: layertype |
|---|
| 18 | real(8), intent(IN) :: porosity, fill, rhocobs, tiobs |
|---|
| 19 | real(8), intent(OUT) :: newti, newrhoc |
|---|
| 20 | real(8), intent(IN) :: icefrac |
|---|
| 21 | real(8) kobs, cice, icedensity, kice |
|---|
| 22 | !parameter (cice=2000.d0, icedensity=926.d0, kice=2.4d0) ! unaffected by scaling |
|---|
| 23 | parameter (cice=1540.d0, icedensity=927.d0, kice=3.2d0) ! at 198 Kelvin |
|---|
| 24 | real(8) fA, ki0, ki, k |
|---|
| 25 | real(8), parameter :: kw=3. ! Mellon et al., JGR 102, 19357 (1997) |
|---|
| 26 | |
|---|
| 27 | kobs = tiobs**2/rhocobs |
|---|
| 28 | ! k, rhoc, and ti are defined in between grid points |
|---|
| 29 | ! rhof and T are defined on grid points |
|---|
| 30 | |
|---|
| 31 | newrhoc = -9999. |
|---|
| 32 | newti = -9999. |
|---|
| 33 | |
|---|
| 34 | select case (layertype) |
|---|
| 35 | case (1) ! interstitial ice |
|---|
| 36 | newrhoc = rhocobs + porosity*fill*icedensity*cice |
|---|
| 37 | if (fill>0.) then |
|---|
| 38 | !--linear addition (option A) |
|---|
| 39 | k = porosity*fill*kice + kobs |
|---|
| 40 | !--Mellon et al. 1997 (option B) |
|---|
| 41 | ki0 = porosity/(1/kobs-(1-porosity)/kw) |
|---|
| 42 | fA = sqrt(fill) |
|---|
| 43 | ki = (1-fA)*ki0 + fA*kice |
|---|
| 44 | !k = kw*ki/((1-porosity)*ki+porosity*kw) |
|---|
| 45 | else |
|---|
| 46 | k = kobs |
|---|
| 47 | endif |
|---|
| 48 | newti = sqrt(newrhoc*k) |
|---|
| 49 | |
|---|
| 50 | case (2) ! massive ice (pure or dirty ice) |
|---|
| 51 | newrhoc = rhocobs*(1.-icefrac)/(1.-porosity) + icefrac*icedensity*cice |
|---|
| 52 | k = icefrac*kice + (1.-icefrac)*kw |
|---|
| 53 | newti = sqrt(newrhoc*k) |
|---|
| 54 | |
|---|
| 55 | case (3) ! all ice, special case of layertype 2, which doesn't use porosity |
|---|
| 56 | newrhoc = icedensity*cice |
|---|
| 57 | k = kice |
|---|
| 58 | newti = sqrt(newrhoc*k) |
|---|
| 59 | |
|---|
| 60 | case (4) ! pores completely filled with ice, special case of layertype 1 |
|---|
| 61 | newrhoc = rhocobs + porosity*icedensity*cice |
|---|
| 62 | k = porosity*kice + kobs ! option A, end-member case of type 1, option A |
|---|
| 63 | !k = kw*kice/((1-porosity)*kice+porosity*kw) ! option B, harmonic average |
|---|
| 64 | newti = sqrt(newrhoc*k) |
|---|
| 65 | |
|---|
| 66 | case (5) ! custom values |
|---|
| 67 | ! values from Mellon et al. (2004) for ice-cemented soil |
|---|
| 68 | newrhoc = 2018.*1040. |
|---|
| 69 | k = 2.5 |
|---|
| 70 | newti = sqrt(newrhoc*k) |
|---|
| 71 | |
|---|
| 72 | case default |
|---|
| 73 | error stop 'invalid layer type' |
|---|
| 74 | |
|---|
| 75 | end select |
|---|
| 76 | |
|---|
| 77 | end subroutine soilthprop |
|---|
