[3286] | 1 | MODULE layering_mod |
---|
| 2 | |
---|
| 3 | !======================================================================= |
---|
[3778] | 4 | ! Purpose: Manage layered layerings_map in the PEM with a linked list data structure |
---|
[3286] | 5 | ! |
---|
| 6 | ! Author: JBC |
---|
| 7 | ! Date: 08/03/2024 |
---|
| 8 | !======================================================================= |
---|
| 9 | |
---|
[3553] | 10 | use glaciers_mod, only: rho_co2ice, rho_h2oice |
---|
[3286] | 11 | |
---|
| 12 | implicit none |
---|
| 13 | |
---|
| 14 | !---- Declarations |
---|
| 15 | ! Numerical parameter |
---|
| 16 | real, parameter :: eps = epsilon(1.), tol = 1.e2*eps |
---|
[3297] | 17 | integer :: nb_str_max |
---|
[3319] | 18 | logical :: layering_algo |
---|
[3286] | 19 | |
---|
| 20 | ! Physical parameters |
---|
[3821] | 21 | logical :: impose_dust_ratio = .false. ! Impose dust-to-ice ratio instead of dust tendency (see 'dust2ice_ratio') |
---|
| 22 | real :: dust2ice_ratio = 0.01 ! Dust-to-ice ratio when ice condenses (1% by default) |
---|
| 23 | real :: d_dust = 1.e-3 ! Tendency of dust [kg.m-2.y-1] |
---|
| 24 | real, parameter :: dry_lag_porosity = 0.4 ! Porosity of dust lag |
---|
| 25 | real, parameter :: wet_lag_porosity = 0.1 ! Porosity of water ice lag |
---|
| 26 | real, parameter :: regolith_porosity = 0.4 ! Porosity of regolith |
---|
| 27 | real, parameter :: breccia_porosity = 0.1 ! Porosity of breccia |
---|
| 28 | real, parameter :: co2ice_porosity = 0. ! Porosity of CO2 ice |
---|
| 29 | real, parameter :: h2oice_porosity = 0. ! Porosity of H2O ice |
---|
| 30 | real, parameter :: rho_dust = 2500. ! Density of dust [kg.m-3] |
---|
| 31 | real, parameter :: rho_rock = 3200. ! Density of rock [kg.m-3] |
---|
| 32 | real, parameter :: h_patchy_dust = 5.e-4 ! Height under which a dust lag is considered patchy [m] |
---|
| 33 | real, parameter :: h_patchy_ice = 5.e-4 ! Height under which a H2O ice lag is considered patchy [m] |
---|
[3286] | 34 | |
---|
[3782] | 35 | ! Parameters for CO2 flux correction (inspired by Levrard et al. 2007) |
---|
[3809] | 36 | real, parameter :: hmin_lag_co2 = 0.5 ! Minimal height of the lag deposit to reduce the sublimation rate [m] |
---|
| 37 | real, parameter :: fred_sublco2 = 0.1 ! Reduction factor of sublimation rate |
---|
[3424] | 38 | |
---|
[3286] | 39 | ! Stratum = node of the linked list |
---|
| 40 | type :: stratum |
---|
[3673] | 41 | real :: top_elevation ! Layer top_elevation (top height from the surface) [m] |
---|
[3770] | 42 | real :: h_co2ice ! Height of CO2 ice [m] |
---|
| 43 | real :: h_h2oice ! Height of H2O ice [m] |
---|
| 44 | real :: h_dust ! Height of dust [m] |
---|
| 45 | real :: h_pore ! Height of pores [m] |
---|
[3673] | 46 | real :: poreice_volfrac ! Pore-filled ice volumetric fraction |
---|
| 47 | type(stratum), pointer :: up => null() ! Upper stratum (next node) |
---|
| 48 | type(stratum), pointer :: down => null() ! Lower stratum (previous node) |
---|
[3286] | 49 | end type stratum |
---|
| 50 | |
---|
| 51 | ! Layering = linked list |
---|
| 52 | type layering |
---|
| 53 | integer :: nb_str = 0 ! Number of strata |
---|
| 54 | type(stratum), pointer :: bottom => null() ! First stratum (head of the list) |
---|
| 55 | type(stratum), pointer :: top => null() ! Last stratum (tail of the list) |
---|
[3778] | 56 | !real :: h_toplag ! Height of dust layer at the top [m] |
---|
| 57 | !type(stratum), pointer :: subice => null() ! Shallower stratum containing ice |
---|
| 58 | !type(stratum), pointer :: surf => null() ! Surface stratum |
---|
[3286] | 59 | end type layering |
---|
| 60 | |
---|
[3297] | 61 | ! Array of pointers to "stratum" |
---|
| 62 | type ptrarray |
---|
| 63 | type(stratum), pointer :: p |
---|
| 64 | end type ptrarray |
---|
| 65 | |
---|
[3286] | 66 | !======================================================================= |
---|
| 67 | contains |
---|
[3297] | 68 | ! Procedures to manipulate the data structure: |
---|
| 69 | ! > print_stratum |
---|
| 70 | ! > add_stratum |
---|
| 71 | ! > insert_stratum |
---|
[3782] | 72 | ! > del_stratum |
---|
[3297] | 73 | ! > get_stratum |
---|
| 74 | ! > ini_layering |
---|
| 75 | ! > del_layering |
---|
| 76 | ! > print_layering |
---|
| 77 | ! > get_nb_str_max |
---|
| 78 | ! > stratif2array |
---|
| 79 | ! > array2stratif |
---|
[3778] | 80 | ! > print_layerings_map |
---|
[3770] | 81 | ! Procedures to get information about a stratum: |
---|
| 82 | ! > thickness |
---|
| 83 | ! > is_co2ice_str |
---|
| 84 | ! > is_h2oice_str |
---|
| 85 | ! > is_dust_str |
---|
| 86 | ! > is_dust_lag |
---|
[3778] | 87 | ! > subsurface_ice_layering |
---|
[3297] | 88 | ! Procedures for the algorithm to evolve the layering: |
---|
| 89 | ! > make_layering |
---|
[3785] | 90 | ! > lift_dust |
---|
[3770] | 91 | ! > subl_ice_str |
---|
[3286] | 92 | !======================================================================= |
---|
| 93 | ! To display a stratum |
---|
| 94 | SUBROUTINE print_stratum(this) |
---|
| 95 | |
---|
| 96 | implicit none |
---|
| 97 | |
---|
| 98 | !---- Arguments |
---|
| 99 | type(stratum), intent(in) :: this |
---|
| 100 | |
---|
| 101 | !---- Code |
---|
[3770] | 102 | write(*,'(a,es13.6)') 'Top elevation [m]: ', this%top_elevation |
---|
| 103 | write(*,'(a,es13.6)') 'Height of CO2 ice [m]: ', this%h_co2ice |
---|
| 104 | write(*,'(a,es13.6)') 'Height of H2O ice [m]: ', this%h_h2oice |
---|
| 105 | write(*,'(a,es13.6)') 'Height of dust [m]: ', this%h_dust |
---|
| 106 | write(*,'(a,es13.6)') 'Height of pores [m]: ', this%h_pore |
---|
| 107 | write(*,'(a,es13.6)') 'Pore-filled ice [m3/m3]: ', this%poreice_volfrac |
---|
[3286] | 108 | |
---|
| 109 | END SUBROUTINE print_stratum |
---|
| 110 | |
---|
| 111 | !======================================================================= |
---|
| 112 | ! To add a stratum to the top of a layering |
---|
[3770] | 113 | SUBROUTINE add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice) |
---|
[3286] | 114 | |
---|
| 115 | implicit none |
---|
| 116 | |
---|
| 117 | !---- Arguments |
---|
| 118 | type(layering), intent(inout) :: this |
---|
[3770] | 119 | real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice |
---|
[3286] | 120 | |
---|
| 121 | !---- Local variables |
---|
| 122 | type(stratum), pointer :: str |
---|
| 123 | |
---|
| 124 | !---- Code |
---|
| 125 | ! Creation of the stratum |
---|
| 126 | allocate(str) |
---|
| 127 | nullify(str%up,str%down) |
---|
[3778] | 128 | str%top_elevation = 1. |
---|
| 129 | str%h_co2ice = 0. |
---|
| 130 | str%h_h2oice = 0. |
---|
| 131 | str%h_dust = 1. |
---|
| 132 | str%h_pore = 0. |
---|
[3673] | 133 | str%poreice_volfrac = 0. |
---|
[3778] | 134 | if (present(top_elevation)) str%top_elevation = top_elevation |
---|
| 135 | if (present(co2ice)) str%h_co2ice = co2ice |
---|
| 136 | if (present(h2oice)) str%h_h2oice = h2oice |
---|
| 137 | if (present(dust)) str%h_dust = dust |
---|
| 138 | if (present(pore)) str%h_pore = pore |
---|
| 139 | if (present(poreice)) str%poreice_volfrac = poreice |
---|
[3286] | 140 | |
---|
| 141 | ! Increment the number of strata |
---|
| 142 | this%nb_str = this%nb_str + 1 |
---|
| 143 | |
---|
| 144 | ! Add the stratum to the layering |
---|
| 145 | if (.not. associated(this%bottom)) then ! If it is the very first one |
---|
| 146 | this%bottom => str |
---|
| 147 | this%top => str |
---|
| 148 | else |
---|
| 149 | str%down => this%top |
---|
| 150 | this%top%up => str |
---|
| 151 | this%top => str |
---|
| 152 | endif |
---|
| 153 | |
---|
| 154 | END SUBROUTINE add_stratum |
---|
| 155 | |
---|
| 156 | !======================================================================= |
---|
| 157 | ! To insert a stratum after another one in a layering |
---|
[3770] | 158 | SUBROUTINE insert_stratum(this,str_prev,top_elevation,co2ice,h2oice,dust,pore,poreice) |
---|
[3286] | 159 | |
---|
| 160 | implicit none |
---|
| 161 | |
---|
| 162 | !---- Arguments |
---|
| 163 | type(layering), intent(inout) :: this |
---|
| 164 | type(stratum), pointer, intent(inout) :: str_prev |
---|
[3770] | 165 | real, intent(in), optional :: top_elevation, co2ice, h2oice, dust, pore, poreice |
---|
[3286] | 166 | |
---|
| 167 | !---- Local variables |
---|
| 168 | type(stratum), pointer :: str, current |
---|
| 169 | |
---|
| 170 | !---- Code |
---|
| 171 | if (.not. associated(str_prev%up)) then ! If str_prev is the top stratum of the layering |
---|
[3770] | 172 | call add_stratum(this,top_elevation,co2ice,h2oice,dust,pore,poreice) |
---|
[3286] | 173 | else |
---|
| 174 | ! Creation of the stratum |
---|
| 175 | allocate(str) |
---|
| 176 | nullify(str%up,str%down) |
---|
[3778] | 177 | str%top_elevation = 1. |
---|
| 178 | str%h_co2ice = 0. |
---|
| 179 | str%h_h2oice = 0. |
---|
| 180 | str%h_dust = 1. |
---|
| 181 | str%h_pore = 0. |
---|
[3673] | 182 | str%poreice_volfrac = 0. |
---|
[3778] | 183 | if (present(top_elevation)) str%top_elevation = top_elevation |
---|
| 184 | if (present(co2ice)) str%h_co2ice = co2ice |
---|
| 185 | if (present(h2oice)) str%h_h2oice = h2oice |
---|
| 186 | if (present(dust)) str%h_dust = dust |
---|
| 187 | if (present(pore)) str%h_pore = pore |
---|
| 188 | if (present(poreice)) str%poreice_volfrac = poreice |
---|
[3286] | 189 | |
---|
| 190 | ! Increment the number of strata |
---|
| 191 | this%nb_str = this%nb_str + 1 |
---|
| 192 | |
---|
| 193 | ! Insert the stratum into the layering at the right place |
---|
| 194 | str%down => str_prev |
---|
| 195 | str%up => str_prev%up |
---|
| 196 | str_prev%up => str |
---|
| 197 | str%up%down => str |
---|
| 198 | |
---|
[3297] | 199 | ! Correct the value of top_elevation for the upper strata |
---|
[3286] | 200 | current => str%up |
---|
| 201 | do while (associated(current)) |
---|
[3770] | 202 | current%top_elevation = current%down%top_elevation + thickness(current) |
---|
[3286] | 203 | current => current%up |
---|
| 204 | enddo |
---|
| 205 | endif |
---|
| 206 | |
---|
| 207 | END SUBROUTINE insert_stratum |
---|
| 208 | |
---|
| 209 | !======================================================================= |
---|
| 210 | ! To remove a stratum in a layering |
---|
[3782] | 211 | SUBROUTINE del_stratum(this,str) |
---|
[3286] | 212 | |
---|
| 213 | implicit none |
---|
| 214 | |
---|
| 215 | !---- Arguments |
---|
| 216 | type(layering), intent(inout) :: this |
---|
| 217 | type(stratum), pointer, intent(inout) :: str |
---|
| 218 | |
---|
| 219 | !---- Local variables |
---|
[3809] | 220 | type(stratum), pointer :: new_str ! To make the argument point towards something |
---|
[3286] | 221 | |
---|
| 222 | !---- Code |
---|
[3782] | 223 | ! Protect the sub-surface strata from deletion |
---|
[3297] | 224 | if (str%top_elevation <= 0.) return |
---|
[3286] | 225 | |
---|
| 226 | ! Decrement the number of strata |
---|
| 227 | this%nb_str = this%nb_str - 1 |
---|
| 228 | |
---|
[3809] | 229 | ! Making the new links |
---|
| 230 | if (associated(str%down)) then ! If there is a lower stratum |
---|
| 231 | str%down%up => str%up |
---|
| 232 | else ! If it was the first stratum (bottom of layering) |
---|
| 233 | this%bottom => str%up |
---|
| 234 | endif |
---|
| 235 | if (associated(str%up)) then ! If there is an upper stratum |
---|
[3286] | 236 | str%up%down => str%down |
---|
[3809] | 237 | else ! If it was the last stratum (top of layering) |
---|
[3286] | 238 | this%top => str%down |
---|
| 239 | endif |
---|
[3809] | 240 | |
---|
| 241 | ! Remove the stratum from the layering |
---|
[3286] | 242 | new_str => str%down |
---|
| 243 | nullify(str%up,str%down) |
---|
| 244 | deallocate(str) |
---|
| 245 | nullify(str) |
---|
| 246 | str => new_str |
---|
| 247 | |
---|
[3782] | 248 | END SUBROUTINE del_stratum |
---|
[3286] | 249 | |
---|
| 250 | !======================================================================= |
---|
[3297] | 251 | ! To get a specific stratum in a layering |
---|
| 252 | FUNCTION get_stratum(lay,i) RESULT(str) |
---|
| 253 | |
---|
| 254 | implicit none |
---|
| 255 | |
---|
| 256 | !---- Arguments |
---|
| 257 | type(layering), intent(in) :: lay |
---|
| 258 | integer, intent(in) :: i |
---|
| 259 | type(stratum), pointer :: str |
---|
| 260 | |
---|
| 261 | !---- Local variables |
---|
| 262 | integer :: k |
---|
| 263 | |
---|
| 264 | !---- Code |
---|
| 265 | if (i < 1 .or. i > lay%nb_str) error stop 'get_stratum: requested number is beyond the number of strata of the layering!' |
---|
| 266 | k = 1 |
---|
| 267 | str => lay%bottom |
---|
| 268 | do while (k /= i) |
---|
| 269 | str => str%up |
---|
| 270 | k = k + 1 |
---|
| 271 | enddo |
---|
| 272 | |
---|
| 273 | END FUNCTION get_stratum |
---|
| 274 | |
---|
| 275 | !======================================================================= |
---|
[3286] | 276 | ! To initialize the layering |
---|
| 277 | SUBROUTINE ini_layering(this) |
---|
| 278 | |
---|
[3770] | 279 | use comsoil_h_PEM, only: soil_pem, nsoilmx_PEM, layer_PEM, index_breccia, index_bedrock |
---|
| 280 | |
---|
[3286] | 281 | implicit none |
---|
| 282 | |
---|
| 283 | !---- Arguments |
---|
| 284 | type(layering), intent(inout) :: this |
---|
| 285 | |
---|
[3770] | 286 | !---- Local variables |
---|
| 287 | integer :: i |
---|
| 288 | real :: h_soil, h_pore, h_tot |
---|
| 289 | |
---|
[3286] | 290 | !---- Code |
---|
[3770] | 291 | ! Creation of strata at the bottom of the layering to describe the sub-surface |
---|
| 292 | if (soil_pem) then |
---|
| 293 | do i = nsoilmx_PEM,index_bedrock,-1 |
---|
| 294 | h_soil = layer_PEM(i) - layer_PEM(i - 1) ! No porosity |
---|
| 295 | call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,0.,0.) |
---|
| 296 | enddo |
---|
| 297 | do i = index_bedrock - 1,index_breccia,-1 |
---|
| 298 | h_tot = layer_PEM(i) - layer_PEM(i - 1) |
---|
| 299 | h_pore = h_tot*breccia_porosity |
---|
| 300 | h_soil = h_tot - h_pore |
---|
| 301 | call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.) |
---|
| 302 | enddo |
---|
| 303 | do i = index_breccia - 1,2,-1 |
---|
| 304 | h_tot = layer_PEM(i) - layer_PEM(i - 1) |
---|
[3785] | 305 | h_pore = h_tot*regolith_porosity |
---|
[3770] | 306 | h_soil = h_tot - h_pore |
---|
| 307 | call add_stratum(this,-layer_PEM(i - 1),0.,0.,h_soil,h_pore,0.) |
---|
| 308 | enddo |
---|
[3785] | 309 | h_pore = layer_PEM(1)*regolith_porosity |
---|
[3770] | 310 | h_soil = layer_PEM(1) - h_pore |
---|
| 311 | call add_stratum(this,0.,0.,0.,h_soil,h_pore,0.) |
---|
| 312 | endif |
---|
[3286] | 313 | |
---|
| 314 | END SUBROUTINE ini_layering |
---|
| 315 | |
---|
| 316 | !======================================================================= |
---|
| 317 | ! To delete the entire layering |
---|
| 318 | SUBROUTINE del_layering(this) |
---|
| 319 | |
---|
| 320 | implicit none |
---|
| 321 | |
---|
| 322 | !---- Arguments |
---|
| 323 | type(layering), intent(inout) :: this |
---|
| 324 | |
---|
| 325 | !---- Local variables |
---|
| 326 | type(stratum), pointer :: current, next |
---|
| 327 | |
---|
| 328 | !---- Code |
---|
| 329 | current => this%bottom |
---|
| 330 | do while (associated(current)) |
---|
| 331 | next => current%up |
---|
| 332 | deallocate(current) |
---|
| 333 | nullify(current) |
---|
| 334 | current => next |
---|
| 335 | enddo |
---|
| 336 | this%nb_str = 0 |
---|
| 337 | |
---|
| 338 | END SUBROUTINE del_layering |
---|
| 339 | |
---|
| 340 | !======================================================================= |
---|
[3297] | 341 | ! To display a layering |
---|
| 342 | SUBROUTINE print_layering(this) |
---|
[3286] | 343 | |
---|
| 344 | implicit none |
---|
| 345 | |
---|
| 346 | !---- Arguments |
---|
[3297] | 347 | type(layering), intent(in) :: this |
---|
[3286] | 348 | |
---|
| 349 | !---- Local variables |
---|
[3297] | 350 | type(stratum), pointer :: current |
---|
| 351 | integer :: i |
---|
[3286] | 352 | |
---|
| 353 | !---- Code |
---|
[3297] | 354 | current => this%top |
---|
| 355 | |
---|
| 356 | i = this%nb_str |
---|
| 357 | do while (associated(current)) |
---|
[3673] | 358 | write(*,'(a,i4)') 'Stratum ', i |
---|
[3297] | 359 | call print_stratum(current) |
---|
| 360 | write(*,*) |
---|
| 361 | current => current%down |
---|
| 362 | i = i - 1 |
---|
[3286] | 363 | enddo |
---|
| 364 | |
---|
[3297] | 365 | END SUBROUTINE print_layering |
---|
[3286] | 366 | |
---|
| 367 | !======================================================================= |
---|
[3297] | 368 | ! To get the maximum number of "stratum" in the stratification (layerings) |
---|
[3778] | 369 | FUNCTION get_nb_str_max(layerings_map,ngrid,nslope) RESULT(nb_str_max) |
---|
[3286] | 370 | |
---|
| 371 | implicit none |
---|
| 372 | |
---|
| 373 | !---- Arguments |
---|
[3778] | 374 | type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map |
---|
[3297] | 375 | integer, intent(in) :: ngrid, nslope |
---|
| 376 | integer :: nb_str_max |
---|
[3286] | 377 | |
---|
| 378 | !---- Local variables |
---|
[3297] | 379 | integer :: ig, islope |
---|
| 380 | |
---|
| 381 | !---- Code |
---|
| 382 | nb_str_max = 0 |
---|
| 383 | do islope = 1,nslope |
---|
| 384 | do ig = 1,ngrid |
---|
[3778] | 385 | nb_str_max = max(layerings_map(ig,islope)%nb_str,nb_str_max) |
---|
[3297] | 386 | enddo |
---|
| 387 | enddo |
---|
| 388 | |
---|
| 389 | END FUNCTION get_nb_str_max |
---|
| 390 | |
---|
| 391 | !======================================================================= |
---|
[3778] | 392 | ! To convert the layerings map into an array able to be outputted |
---|
| 393 | SUBROUTINE stratif2array(layerings_map,ngrid,nslope,stratif_array) |
---|
[3297] | 394 | |
---|
| 395 | implicit none |
---|
| 396 | |
---|
| 397 | !---- Arguments |
---|
| 398 | integer, intent(in) :: ngrid, nslope |
---|
[3778] | 399 | type(layering), dimension(ngrid,nslope), intent(in) :: layerings_map |
---|
[3297] | 400 | real, dimension(:,:,:,:), allocatable, intent(inout) :: stratif_array |
---|
| 401 | |
---|
| 402 | !---- Local variables |
---|
[3286] | 403 | type(stratum), pointer :: current |
---|
[3297] | 404 | integer :: ig, islope, k |
---|
[3286] | 405 | |
---|
| 406 | !---- Code |
---|
[3297] | 407 | stratif_array = 0. |
---|
| 408 | do islope = 1,nslope |
---|
| 409 | do ig = 1,ngrid |
---|
[3778] | 410 | current => layerings_map(ig,islope)%bottom |
---|
[3297] | 411 | k = 1 |
---|
| 412 | do while (associated(current)) |
---|
[3770] | 413 | stratif_array(ig,islope,k,1) = current%top_elevation |
---|
| 414 | stratif_array(ig,islope,k,2) = current%h_co2ice |
---|
| 415 | stratif_array(ig,islope,k,3) = current%h_h2oice |
---|
| 416 | stratif_array(ig,islope,k,4) = current%h_dust |
---|
| 417 | stratif_array(ig,islope,k,5) = current%h_pore |
---|
| 418 | stratif_array(ig,islope,k,6) = current%poreice_volfrac |
---|
[3297] | 419 | current => current%up |
---|
| 420 | k = k + 1 |
---|
| 421 | enddo |
---|
| 422 | enddo |
---|
[3286] | 423 | enddo |
---|
| 424 | |
---|
[3297] | 425 | END SUBROUTINE stratif2array |
---|
[3286] | 426 | |
---|
| 427 | !======================================================================= |
---|
[3778] | 428 | ! To convert the stratification array into the layerings map |
---|
| 429 | SUBROUTINE array2stratif(stratif_array,ngrid,nslope,layerings_map) |
---|
[3286] | 430 | |
---|
| 431 | implicit none |
---|
| 432 | |
---|
| 433 | !---- Arguments |
---|
[3297] | 434 | integer, intent(in) :: ngrid, nslope |
---|
| 435 | real, dimension(:,:,:,:), allocatable, intent(in) :: stratif_array |
---|
[3778] | 436 | type(layering), dimension(ngrid,nslope), intent(inout) :: layerings_map |
---|
[3286] | 437 | |
---|
| 438 | !---- Local variables |
---|
[3430] | 439 | integer :: ig, islope, k |
---|
[3286] | 440 | |
---|
| 441 | !---- Code |
---|
[3297] | 442 | do islope = 1,nslope |
---|
| 443 | do ig = 1,ngrid |
---|
[3770] | 444 | do k = 1,size(stratif_array,3) |
---|
[3778] | 445 | call add_stratum(layerings_map(ig,islope),stratif_array(ig,islope,k,1),stratif_array(ig,islope,k,2),stratif_array(ig,islope,k,3),stratif_array(ig,islope,k,4),stratif_array(ig,islope,k,5),stratif_array(ig,islope,k,6)) |
---|
[3297] | 446 | enddo |
---|
| 447 | enddo |
---|
[3286] | 448 | enddo |
---|
| 449 | |
---|
[3297] | 450 | END SUBROUTINE array2stratif |
---|
[3286] | 451 | |
---|
| 452 | !======================================================================= |
---|
[3778] | 453 | ! To display the map of layerings |
---|
| 454 | SUBROUTINE print_layerings_map(this,ngrid,nslope) |
---|
[3286] | 455 | |
---|
| 456 | implicit none |
---|
| 457 | |
---|
| 458 | !---- Arguments |
---|
[3297] | 459 | integer, intent(in) :: ngrid, nslope |
---|
| 460 | type(layering), dimension(ngrid,nslope), intent(in) :: this |
---|
[3286] | 461 | |
---|
[3297] | 462 | !---- Local variables |
---|
| 463 | integer :: ig, islope |
---|
| 464 | |
---|
[3286] | 465 | !---- Code |
---|
[3297] | 466 | do ig = 1,ngrid |
---|
| 467 | write(*,'(a10,i4)') 'Grid cell ', ig |
---|
| 468 | do islope = 1,nslope |
---|
| 469 | write(*,'(a13,i1)') 'Slope number ', islope |
---|
| 470 | call print_layering(this(ig,islope)) |
---|
| 471 | write(*,*) '~~~~~~~~~~' |
---|
| 472 | enddo |
---|
| 473 | write(*,*) '--------------------' |
---|
| 474 | enddo |
---|
[3286] | 475 | |
---|
[3778] | 476 | END SUBROUTINE print_layerings_map |
---|
[3286] | 477 | |
---|
| 478 | !======================================================================= |
---|
[3770] | 479 | !======================================================================= |
---|
| 480 | ! To compute the thickness of a stratum |
---|
| 481 | FUNCTION thickness(this) RESULT(h_str) |
---|
| 482 | |
---|
| 483 | implicit none |
---|
| 484 | |
---|
| 485 | !---- Arguments |
---|
| 486 | type(stratum), intent(in) :: this |
---|
| 487 | real :: h_str |
---|
| 488 | |
---|
| 489 | !---- Code |
---|
| 490 | h_str = this%h_co2ice + this%h_h2oice + this%h_dust + this%h_pore |
---|
| 491 | |
---|
| 492 | END FUNCTION thickness |
---|
| 493 | |
---|
| 494 | !======================================================================= |
---|
| 495 | ! To know whether the stratum can be considered as a CO2 ice layer |
---|
| 496 | FUNCTION is_co2ice_str(str) RESULT(co2ice_str) |
---|
| 497 | |
---|
| 498 | implicit none |
---|
| 499 | |
---|
| 500 | !---- Arguments |
---|
| 501 | type(stratum), intent(in) :: str |
---|
| 502 | logical :: co2ice_str |
---|
| 503 | |
---|
| 504 | !---- Code |
---|
| 505 | co2ice_str = .false. |
---|
| 506 | if (str%h_co2ice > str%h_h2oice .and. str%h_co2ice > str%h_dust) co2ice_str = .true. |
---|
| 507 | |
---|
| 508 | END FUNCTION is_co2ice_str |
---|
| 509 | |
---|
| 510 | !======================================================================= |
---|
| 511 | ! To know whether the stratum can be considered as a H2O ice layer |
---|
| 512 | FUNCTION is_h2oice_str(str) RESULT(h2oice_str) |
---|
| 513 | |
---|
| 514 | implicit none |
---|
| 515 | |
---|
| 516 | !---- Arguments |
---|
| 517 | type(stratum), intent(in) :: str |
---|
| 518 | logical :: h2oice_str |
---|
| 519 | |
---|
| 520 | !---- Code |
---|
| 521 | h2oice_str = .false. |
---|
| 522 | if (str%h_h2oice > str%h_co2ice .and. str%h_h2oice > str%h_dust) h2oice_str = .true. |
---|
| 523 | |
---|
| 524 | END FUNCTION is_h2oice_str |
---|
| 525 | |
---|
| 526 | !======================================================================= |
---|
| 527 | ! To know whether the stratum can be considered as a dust layer |
---|
| 528 | FUNCTION is_dust_str(str) RESULT(dust_str) |
---|
| 529 | |
---|
| 530 | implicit none |
---|
| 531 | |
---|
| 532 | !---- Arguments |
---|
| 533 | type(stratum), intent(in) :: str |
---|
| 534 | logical :: dust_str |
---|
| 535 | |
---|
| 536 | !---- Code |
---|
| 537 | dust_str = .false. |
---|
| 538 | if (str%h_dust > str%h_co2ice .and. str%h_dust > str%h_h2oice) dust_str = .true. |
---|
| 539 | |
---|
| 540 | END FUNCTION is_dust_str |
---|
| 541 | |
---|
| 542 | !======================================================================= |
---|
| 543 | ! To know whether the stratum can be considered as a dust lag |
---|
| 544 | FUNCTION is_dust_lag(str) RESULT(dust_lag) |
---|
| 545 | |
---|
| 546 | implicit none |
---|
| 547 | |
---|
| 548 | !---- Arguments |
---|
| 549 | type(stratum), intent(in) :: str |
---|
| 550 | logical :: dust_lag |
---|
| 551 | |
---|
| 552 | !---- Code |
---|
| 553 | dust_lag = .false. |
---|
| 554 | if (str%h_co2ice < tol .and. str%h_h2oice < tol .and. str%poreice_volfrac < tol) dust_lag = .true. |
---|
| 555 | |
---|
| 556 | END FUNCTION is_dust_lag |
---|
| 557 | |
---|
| 558 | !======================================================================= |
---|
[3778] | 559 | ! To get data about possible subsurface ice in a layering |
---|
[3809] | 560 | SUBROUTINE subsurface_ice_layering(this,h_toplag,h2o_ice,co2_ice) |
---|
[3770] | 561 | |
---|
| 562 | implicit none |
---|
| 563 | |
---|
| 564 | !---- Arguments |
---|
| 565 | type(layering), intent(in) :: this |
---|
[3809] | 566 | real, intent(out) :: h2o_ice, co2_ice, h_toplag |
---|
[3770] | 567 | |
---|
| 568 | !---- Local variables |
---|
| 569 | type(stratum), pointer :: current |
---|
| 570 | |
---|
| 571 | !---- Code |
---|
| 572 | h_toplag = 0. |
---|
[3778] | 573 | h2o_ice = 0. |
---|
[3809] | 574 | co2_ice = 0. |
---|
[3770] | 575 | current => this%top |
---|
| 576 | ! Is the considered stratum a dust lag? |
---|
[3809] | 577 | if (.not. is_dust_lag(current)) return |
---|
[3821] | 578 | do while (is_dust_lag(current)) |
---|
[3770] | 579 | h_toplag = h_toplag + thickness(current) |
---|
[3821] | 580 | if (associated(current%down)) then |
---|
| 581 | current => current%down |
---|
| 582 | else |
---|
| 583 | exit |
---|
| 584 | endif |
---|
[3770] | 585 | enddo |
---|
[3821] | 586 | |
---|
| 587 | if (.not. associated(current%down)) then ! Bottom of the layering -> no ice below |
---|
| 588 | h_toplag = 0. |
---|
| 589 | return |
---|
| 590 | endif |
---|
| 591 | |
---|
[3778] | 592 | ! Is the stratum below the dust lag made of ice? |
---|
[3821] | 593 | if (current%h_h2oice > 0. .and. current%h_h2oice > current%h_co2ice) then ! Yes, there is more H2O ice than CO2 ice |
---|
[3809] | 594 | if (h_toplag < h_patchy_dust) then ! But the dust lag is too thin to considered ice as subsurface ice |
---|
| 595 | h_toplag = 0. |
---|
| 596 | h2o_ice = current%h_h2oice |
---|
| 597 | endif |
---|
[3821] | 598 | else if (current%h_co2ice > 0. .and. current%h_co2ice > current%h_h2oice) then ! No, there is more CO2 ice than H2O ice |
---|
[3809] | 599 | h_toplag = 0. ! Because it matters only for H2O ice |
---|
| 600 | if (h_toplag < h_patchy_ice) then ! But the dust lag is too thin to considered ice as subsurface ice |
---|
| 601 | co2_ice = current%h_co2ice |
---|
| 602 | endif |
---|
[3778] | 603 | endif |
---|
[3770] | 604 | |
---|
[3778] | 605 | END SUBROUTINE subsurface_ice_layering |
---|
[3770] | 606 | |
---|
| 607 | !======================================================================= |
---|
| 608 | !======================================================================= |
---|
[3785] | 609 | ! To lift dust in a stratum |
---|
| 610 | SUBROUTINE lift_dust(this,str,h2lift) |
---|
| 611 | |
---|
| 612 | implicit none |
---|
| 613 | |
---|
| 614 | !---- Arguments |
---|
| 615 | type(layering), intent(inout) :: this |
---|
| 616 | type(stratum), pointer, intent(inout) :: str |
---|
| 617 | real, intent(in) :: h2lift |
---|
| 618 | |
---|
| 619 | !---- Code |
---|
| 620 | ! Update of properties in the eroding dust lag |
---|
| 621 | str%top_elevation = str%top_elevation - h2lift*(1. + str%h_pore/str%h_dust) |
---|
| 622 | str%h_pore = str%h_pore - h2lift*str%h_pore/str%h_dust |
---|
| 623 | str%h_dust = str%h_dust - h2lift |
---|
| 624 | |
---|
| 625 | ! Remove the eroding dust lag if there is no dust anymore |
---|
| 626 | if (str%h_dust < tol) call del_stratum(this,str) |
---|
| 627 | |
---|
| 628 | END SUBROUTINE lift_dust |
---|
| 629 | |
---|
| 630 | !======================================================================= |
---|
| 631 | ! To sublimate ice in a stratum |
---|
| 632 | SUBROUTINE sublimate_ice(this,str,h2subl_co2ice,h2subl_h2oice,new_lag,zlag) |
---|
| 633 | |
---|
| 634 | implicit none |
---|
| 635 | |
---|
| 636 | !---- Arguments |
---|
| 637 | type(layering), intent(inout) :: this |
---|
| 638 | type(stratum), pointer, intent(inout) :: str |
---|
| 639 | logical, intent(inout) :: new_lag |
---|
| 640 | real, intent(inout) :: zlag |
---|
| 641 | real, intent(in) :: h2subl_co2ice, h2subl_h2oice |
---|
| 642 | |
---|
| 643 | !---- Local variables |
---|
[3789] | 644 | real :: h2subl, h_ice, h_pore, h_pore_new, h_tot |
---|
[3785] | 645 | real :: hlag_dust, hlag_h2oice |
---|
| 646 | type(stratum), pointer :: current |
---|
| 647 | |
---|
| 648 | !---- Code |
---|
| 649 | h_ice = str%h_co2ice + str%h_h2oice |
---|
[3789] | 650 | h_pore = str%h_pore |
---|
[3785] | 651 | h2subl = h2subl_co2ice + h2subl_h2oice |
---|
| 652 | |
---|
| 653 | ! How much dust and H2O ice does ice sublimation involve to create the lag? |
---|
| 654 | hlag_dust = str%h_dust*h2subl/h_ice |
---|
| 655 | hlag_h2oice = 0. |
---|
| 656 | if (h2subl_co2ice > 0. .and. h2subl_h2oice < tol) hlag_h2oice = str%h_h2oice*h2subl/h_ice |
---|
| 657 | |
---|
| 658 | ! Update of properties in the sublimating stratum |
---|
| 659 | str%top_elevation = str%top_elevation - h2subl*(1. + h_pore/h_ice) - hlag_dust - hlag_h2oice |
---|
| 660 | str%h_co2ice = str%h_co2ice - h2subl_co2ice |
---|
| 661 | str%h_h2oice = str%h_h2oice - h2subl_h2oice - hlag_h2oice |
---|
| 662 | str%h_dust = str%h_dust - hlag_dust |
---|
| 663 | str%h_pore = str%h_pore - h2subl*h_pore/h_ice |
---|
| 664 | |
---|
| 665 | ! Correct the value of top_elevation for the upper strata |
---|
| 666 | current => str%up |
---|
| 667 | do while (associated(current)) |
---|
| 668 | current%top_elevation = current%down%top_elevation + thickness(current) |
---|
| 669 | current => current%up |
---|
| 670 | enddo |
---|
| 671 | |
---|
| 672 | ! Remove the sublimating stratum if there is no ice anymore |
---|
| 673 | if (str%h_co2ice < tol .and. str%h_h2oice < tol) call del_stratum(this,str) |
---|
| 674 | |
---|
| 675 | ! Which porosity is considered? |
---|
| 676 | if (hlag_dust >= hlag_h2oice) then |
---|
[3789] | 677 | h_pore_new = hlag_dust*dry_lag_porosity/(1. - dry_lag_porosity) |
---|
[3785] | 678 | else |
---|
[3789] | 679 | h_pore_new = hlag_h2oice*wet_lag_porosity/(1. - wet_lag_porosity) |
---|
[3785] | 680 | endif |
---|
[3789] | 681 | h_tot = hlag_dust + hlag_h2oice + h_pore_new |
---|
[3785] | 682 | |
---|
| 683 | ! New stratum above the current stratum as a lag layer |
---|
| 684 | if (new_lag) then |
---|
[3789] | 685 | call insert_stratum(this,str,str%top_elevation + h_tot,0.,hlag_h2oice,hlag_dust,h_pore_new,0.) |
---|
[3785] | 686 | new_lag = .false. |
---|
| 687 | else |
---|
| 688 | str%up%top_elevation = str%up%top_elevation + h_tot |
---|
| 689 | str%up%h_h2oice = str%up%h_h2oice + hlag_h2oice |
---|
| 690 | str%up%h_dust = str%up%h_dust + hlag_dust |
---|
[3789] | 691 | str%up%h_pore = str%up%h_pore + h_pore_new |
---|
[3785] | 692 | endif |
---|
| 693 | |
---|
| 694 | zlag = zlag + h_tot |
---|
| 695 | |
---|
| 696 | END SUBROUTINE sublimate_ice |
---|
| 697 | |
---|
| 698 | !======================================================================= |
---|
[3286] | 699 | ! Layering algorithm |
---|
[3770] | 700 | SUBROUTINE make_layering(this,d_co2ice,d_h2oice,new_str,zshift_surf,new_lag,zlag,current) |
---|
[3286] | 701 | |
---|
| 702 | implicit none |
---|
| 703 | |
---|
| 704 | !---- Arguments |
---|
[3553] | 705 | ! Inputs |
---|
| 706 | |
---|
| 707 | ! Outputs |
---|
[3286] | 708 | type(layering), intent(inout) :: this |
---|
[3770] | 709 | type(stratum), pointer, intent(inout) :: current |
---|
| 710 | real, intent(inout) :: d_co2ice, d_h2oice |
---|
| 711 | logical, intent(inout) :: new_str, new_lag |
---|
[3553] | 712 | real, intent(out) :: zshift_surf, zlag |
---|
[3286] | 713 | |
---|
| 714 | !---- Local variables |
---|
[3785] | 715 | real :: dh_co2ice, dh_h2oice, dh_dust |
---|
| 716 | real :: h_co2ice_old, h_h2oice_old, h_dust_old |
---|
| 717 | real :: h2subl_co2ice, h2subl_h2oice, h2subl_co2ice_tot, h2subl_h2oice_tot |
---|
| 718 | real :: h2lift, h2lift_tot, h_pore, h_tot, h_toplag, R_subl |
---|
| 719 | logical :: unexpected |
---|
| 720 | type(stratum), pointer :: currentloc |
---|
[3286] | 721 | |
---|
| 722 | !---- Code |
---|
[3498] | 723 | dh_co2ice = d_co2ice/rho_co2ice |
---|
| 724 | dh_h2oice = d_h2oice/rho_h2oice |
---|
[3789] | 725 | ! To disable dust sedimenation when there is ice sublimation (because dust tendency is fixed) |
---|
[3821] | 726 | if (dh_co2ice >= 0. .and. dh_h2oice >= 0.) then |
---|
| 727 | if (impose_dust_ratio) then |
---|
| 728 | dh_dust = dust2ice_ratio*(dh_h2oice + dh_co2ice) |
---|
| 729 | else |
---|
| 730 | dh_dust = d_dust/rho_dust |
---|
| 731 | endif |
---|
[3789] | 732 | else |
---|
| 733 | dh_dust = 0. |
---|
| 734 | endif |
---|
[3770] | 735 | zshift_surf = this%top%top_elevation |
---|
[3553] | 736 | zlag = 0. |
---|
[3782] | 737 | unexpected = .false. |
---|
[3297] | 738 | |
---|
[3782] | 739 | !----------------------------------------------------------------------- |
---|
[3778] | 740 | ! NOTHING HAPPENS |
---|
| 741 | if (abs(dh_co2ice) < tol .and. abs(dh_h2oice) < tol .and. abs(dh_dust) < tol) then |
---|
| 742 | ! So we do nothing |
---|
| 743 | !----------------------------------------------------------------------- |
---|
[3782] | 744 | ! BUILDING A STRATUM (ice condensation and/or dust desimentation) |
---|
| 745 | else if (dh_co2ice >= 0. .and. dh_h2oice >= 0. .and. dh_dust >= 0.) then |
---|
[3770] | 746 | ! Which porosity is considered? |
---|
[3782] | 747 | if (dh_co2ice > dh_h2oice .and. dh_co2ice > dh_dust) then ! CO2 ice is dominant |
---|
[3770] | 748 | h_pore = dh_co2ice*co2ice_porosity/(1. - co2ice_porosity) |
---|
[3782] | 749 | else if (dh_h2oice > dh_co2ice .and. dh_h2oice > dh_dust) then ! H2O ice is dominant |
---|
[3770] | 750 | h_pore = dh_h2oice*h2oice_porosity/(1. - h2oice_porosity) |
---|
[3782] | 751 | else if (dh_dust > dh_co2ice .and. dh_dust > dh_h2oice) then ! Dust is dominant |
---|
[3785] | 752 | h_pore = dh_dust*regolith_porosity/(1. - regolith_porosity) |
---|
[3782] | 753 | else |
---|
| 754 | unexpected = .true. |
---|
[3286] | 755 | endif |
---|
[3770] | 756 | h_tot = dh_co2ice + dh_h2oice + dh_dust + h_pore |
---|
| 757 | ! New stratum at the top of the layering |
---|
| 758 | if (new_str) then |
---|
| 759 | call add_stratum(this,this%top%top_elevation + h_tot,dh_co2ice,dh_h2oice,dh_dust,h_pore,0.) |
---|
| 760 | new_str = .false. |
---|
| 761 | else |
---|
| 762 | this%top%top_elevation = this%top%top_elevation + h_tot |
---|
[3778] | 763 | this%top%h_co2ice = this%top%h_co2ice + dh_co2ice |
---|
| 764 | this%top%h_h2oice = this%top%h_h2oice + dh_h2oice |
---|
| 765 | this%top%h_dust = this%top%h_dust + dh_dust |
---|
| 766 | this%top%h_pore = this%top%h_pore + h_pore |
---|
[3770] | 767 | endif |
---|
| 768 | !----------------------------------------------------------------------- |
---|
[3785] | 769 | ! RETREATING A STRATUM (ice sublimation and/or dust lifting) |
---|
| 770 | else if (dh_co2ice <= 0. .and. dh_h2oice <= 0. .and. dh_dust <= 0.) then |
---|
[3782] | 771 | h2subl_co2ice_tot = abs(dh_co2ice) |
---|
| 772 | h2subl_h2oice_tot = abs(dh_h2oice) |
---|
[3785] | 773 | h2lift_tot = abs(dh_dust) |
---|
[3789] | 774 | do while ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0. .or. h2lift_tot > 0.) .and. associated(current)) |
---|
[3782] | 775 | h_co2ice_old = current%h_co2ice |
---|
| 776 | h_h2oice_old = current%h_h2oice |
---|
[3785] | 777 | |
---|
| 778 | ! How much is CO2 ice sublimation inhibited by the top dust lag? |
---|
| 779 | R_subl = 1. |
---|
| 780 | if (associated(current%up)) then ! If there is an upper stratum |
---|
| 781 | h_toplag = 0. |
---|
| 782 | currentloc => current%up |
---|
| 783 | do while (is_dust_lag(currentloc) .and. associated(currentloc%up)) |
---|
| 784 | h_toplag = h_toplag + thickness(currentloc) |
---|
| 785 | currentloc => currentloc%up |
---|
| 786 | enddo |
---|
[3789] | 787 | if (currentloc%h_h2oice >= h_patchy_ice) then |
---|
[3785] | 788 | R_subl = 0. |
---|
[3789] | 789 | else if (0. < currentloc%h_h2oice .and. currentloc%h_h2oice < h_patchy_ice) then |
---|
[3785] | 790 | h_toplag = h_toplag + thickness(currentloc) |
---|
[3809] | 791 | R_subl = fred_sublco2**(h_toplag/hmin_lag_co2) |
---|
[3785] | 792 | else |
---|
[3809] | 793 | R_subl = fred_sublco2**(h_toplag/hmin_lag_co2) |
---|
[3785] | 794 | endif |
---|
| 795 | endif |
---|
| 796 | h2subl_co2ice_tot = h2subl_co2ice_tot*R_subl |
---|
| 797 | |
---|
[3782] | 798 | ! Is there CO2 ice to sublimate? |
---|
| 799 | h2subl_co2ice = 0. |
---|
| 800 | if (h_co2ice_old > 0. .and. h2subl_co2ice_tot > 0.) then |
---|
| 801 | h2subl_co2ice = min(h2subl_co2ice_tot,h_co2ice_old) |
---|
| 802 | h2subl_co2ice_tot = h2subl_co2ice_tot - h2subl_co2ice |
---|
| 803 | endif |
---|
| 804 | ! Is there H2O ice to sublimate? |
---|
| 805 | h2subl_h2oice = 0. |
---|
| 806 | if (h_h2oice_old > 0. .and. h2subl_h2oice_tot > 0.) then |
---|
| 807 | h2subl_h2oice = min(h2subl_h2oice_tot,h_h2oice_old) |
---|
| 808 | h2subl_h2oice_tot = h2subl_h2oice_tot - h2subl_h2oice |
---|
| 809 | endif |
---|
[3785] | 810 | ! Ice sublimation |
---|
| 811 | if (h2subl_co2ice > 0. .or. h2subl_h2oice > 0.) call sublimate_ice(this,current,h2subl_co2ice,h2subl_h2oice,new_lag,zlag) |
---|
| 812 | |
---|
| 813 | ! Is there dust to lift? |
---|
| 814 | h2lift = 0. |
---|
| 815 | if (is_dust_lag(current) .and. h2lift_tot > 0.) then |
---|
| 816 | h2lift = min(h2lift_tot,current%h_dust) |
---|
| 817 | h2lift_tot = h2lift_tot - h2lift |
---|
| 818 | ! Dust lifting |
---|
| 819 | if (h2lift > 0.) call lift_dust(this,current,h2lift) |
---|
| 820 | else if (associated(current%up)) then |
---|
| 821 | if (is_dust_lag(current%up) .and. h2lift_tot > 0.) then |
---|
| 822 | h2lift = min(h2lift_tot,current%up%h_dust) |
---|
| 823 | h2lift_tot = h2lift_tot - h2lift |
---|
| 824 | ! Dust lifting |
---|
| 825 | if (h2lift > 0.) call lift_dust(this,current%up,h2lift) |
---|
| 826 | endif |
---|
| 827 | endif |
---|
| 828 | |
---|
| 829 | ! Still ice to be sublimated or dust to be lifted and no more ice in the current stratum so we move to the underlying stratum |
---|
| 830 | if ((h2subl_co2ice_tot > 0. .or. h2subl_h2oice_tot > 0. .or. h2lift_tot > 0.) .and. current%h_co2ice < tol .and. current%h_h2oice < tol .and. current%poreice_volfrac < tol) then |
---|
[3782] | 831 | current => current%down |
---|
| 832 | new_lag = .true. |
---|
| 833 | else |
---|
| 834 | exit |
---|
| 835 | endif |
---|
| 836 | enddo |
---|
| 837 | if (h2subl_co2ice_tot > 0.) dh_co2ice = 0. ! No CO2 ice available anymore so the tendency is set to 0 |
---|
| 838 | if (h2subl_h2oice_tot > 0.) dh_h2oice = 0. ! No H2O ice available anymore so the tendency is set to 0 |
---|
[3785] | 839 | if (h2lift_tot > 0.) dh_dust = 0. ! No dust available anymore so the tendency is set to 0 |
---|
[3782] | 840 | !----------------------------------------------------------------------- |
---|
| 841 | ! CO2 ICE SUBLIMATION + H2O ICE CONDENSATION |
---|
| 842 | else if (dh_co2ice < 0. .and. dh_h2oice > 0.) then |
---|
[3785] | 843 | !~ if (abs(dh_co2ice) > dh_h2oice .and. dh_co2ice < dh_dust .and. dh_dust <= 0.) then ! CO2 sublimation is dominant |
---|
[3286] | 844 | |
---|
[3785] | 845 | !~ else if (dh_h2oice > abs(dh_co2ice) .and. 0 <= dh_dust .and. dh_dust < dh_h2oice) then ! H2O condensation is dominant |
---|
| 846 | |
---|
[3782] | 847 | !~ else |
---|
| 848 | unexpected = .true. |
---|
| 849 | !~ endif |
---|
[3770] | 850 | !----------------------------------------------------------------------- |
---|
[3782] | 851 | ! NOT EXPECTED SITUATION |
---|
[3770] | 852 | else |
---|
[3782] | 853 | unexpected = .true. |
---|
| 854 | endif |
---|
| 855 | |
---|
| 856 | if (unexpected) then |
---|
[3770] | 857 | write(*,'(a)') 'Error: situation for the layering construction not managed!' |
---|
| 858 | write(*,'(a,es12.3)') ' > dh_co2ice [m/y] = ', dh_co2ice |
---|
| 859 | write(*,'(a,es12.3)') ' > dh_h2oice [m/y] = ', dh_h2oice |
---|
[3789] | 860 | write(*,'(a,es12.3)') ' > dh_dust [m/y] = ', dh_dust |
---|
[3782] | 861 | error stop |
---|
[3770] | 862 | endif |
---|
[3424] | 863 | |
---|
[3770] | 864 | zshift_surf = this%top%top_elevation - zshift_surf |
---|
[3424] | 865 | |
---|
[3286] | 866 | END SUBROUTINE make_layering |
---|
| 867 | |
---|
| 868 | END MODULE layering_mod |
---|