Changeset 3553 for trunk/LMDZ.COMMON
- Timestamp:
- Dec 16, 2024, 5:30:48 PM (7 days ago)
- Location:
- trunk/LMDZ.COMMON/libf/evolution
- Files:
-
- 9 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.COMMON/libf/evolution/changelog.txt
r3543 r3553 516 516 == 10/12/2024 == JBC 517 517 Using the correct subroutine to output 'ice_porefilling'. 518 519 == 16/12/2024 == JBC 520 Addition of the shifting of the soil temperature profile to follow the surface evolution due to ice condensation/sublimation + small cleanings/improvements. -
trunk/LMDZ.COMMON/libf/evolution/compute_soiltemp_mod.F90
r3532 r3553 216 216 END SUBROUTINE ini_tsoil_pem 217 217 218 !======================================================================= 219 220 SUBROUTINE shift_tsoil2surf(ngrid,nsoil,nslope,zshift_surf,zlag,tsurf,tsoil) 221 222 use comsoil_h_PEM, only: layer_PEM, mlayer_PEM, fluxgeo, thermdiff_PEM, mthermdiff_PEM 223 use math_mod, only: solve_steady_heat 224 225 implicit none 226 227 !----------------------------------------------------------------------- 228 ! Author: JBC 229 ! Purpose: Shifting the soil temperature profile to follow the surface evolution due to ice condensation/sublimation 230 !----------------------------------------------------------------------- 231 ! Inputs: 232 ! ------- 233 integer, intent(in) :: ngrid ! number of (horizontal) grid-points 234 integer, intent(in) :: nsoil ! number of soil layers 235 integer, intent(in) :: nslope ! number of sub-slopes 236 real, dimension(ngrid,nslope), intent(in) :: zshift_surf ! elevation shift for the surface [m] 237 real, dimension(ngrid,nslope), intent(in) :: zlag ! newly built lag thickness [m] 238 real, dimension(ngrid,nslope), intent(in) :: tsurf ! surface temperature [K] 239 ! Outputs: 240 ! -------- 241 real, dimension(ngrid,nsoil,nslope), intent(inout) :: tsoil ! soil (mid-layer) temperature [K] 242 ! Local: 243 ! ------ 244 integer :: ig, isoil, islope, ishift, ilag, iz 245 real :: a, z, zshift_surfloc 246 real, dimension(ngrid,nsoil,nslope) :: tsoil_old 247 248 write(*,*)"Shifting soil temperature to surface" 249 tsoil_old = tsoil 250 251 do ig = 1,ngrid 252 do islope = 1,nslope 253 zshift_surfloc = zshift_surf(ig,islope) 254 255 if (zshift_surfloc >= 0.) then ! In case of the surface is higher than initially 256 if (zshift_surfloc < mlayer_PEM(0)) then ! Surface change is too small to be taken into account 257 ! Nothing to do; we keep the soil temperature profile 258 else if (zshift_surfloc >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization 259 tsoil(ig,:,islope) = tsurf(ig,islope) 260 else 261 ishift = 0 262 do while (mlayer_PEM(ishift) <= zshift_surfloc) 263 ishift = ishift + 1 ! mlayer indices begin at 0 so this the good index for tsoil! 264 enddo 265 ! The "new soil" temperature is set to tsurf 266 tsoil(ig,:ishift,islope) = tsurf(ig,islope) 267 268 do isoil = ishift + 1,nsoil 269 ! Position in the old discretization of the depth 270 z = mlayer_PEM(isoil - 1) - zshift_surfloc 271 ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs 272 iz = 0 273 do while (mlayer_PEM(iz) <= z) 274 iz = iz + 1 275 enddo 276 ! Interpolation of the temperature profile from the old discretization 277 a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1)) 278 tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope) 279 enddo 280 endif 281 else ! In case of the surface is lower than initially 282 if (abs(zshift_surfloc) < mlayer_PEM(0)) then ! Surface change is too small to be taken into account 283 ! Nothing to do; we keep the soil temperature profile 284 else if (abs(zshift_surfloc) >= mlayer_PEM(nsoil - 1)) then ! Surface change is much larger than the discretization 285 call solve_steady_heat(nsoil,mlayer_PEM,layer_PEM,mthermdiff_PEM(ig,:),thermdiff_PEM(ig,:),tsurf(ig,islope),fluxgeo,tsoil(ig,:,islope)) 286 else 287 if (zlag(ig,islope) < mlayer_PEM(0)) then ! The lag is too thin to be taken into account 288 ilag = 0 289 else 290 ilag = 0 291 do while (mlayer_PEM(ilag) <= zlag(ig,islope)) 292 ilag = ilag + 1 ! mlayer indices begin at 0 so this the good index for tsoil! 293 enddo 294 ! Position of the lag bottom in the old discretization of the depth 295 z = zlag(ig,islope) - zshift_surfloc 296 ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs 297 iz = 0 298 do while (mlayer_PEM(iz) <= z) 299 iz = iz + 1 300 enddo 301 ! The "new lag" temperature is set to the ice temperature just below 302 a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1)) 303 tsoil(ig,:ilag,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope) 304 endif 305 306 ishift = nsoil - 1 307 z = mlayer_PEM(nsoil - 1) + zshift_surfloc 308 do while (mlayer_PEM(ishift) >= z) 309 ishift = ishift - 1 310 enddo 311 ishift = ishift + 1 ! Adding 1 is needed to match the good index for tsoil! 312 do isoil = ilag + 1,ishift 313 ! Position in the old discretization of the depth 314 z = mlayer_PEM(isoil - 1) - zshift_surfloc 315 ! Find the interval [mlayer_PEM(iz - 1),mlayer_PEM(iz)[ where the position z belongs 316 iz = 0 317 do while (mlayer_PEM(iz) <= z) 318 iz = iz + 1 319 enddo 320 ! Interpolation of the temperature profile from the old discretization 321 a = (tsoil_old(ig,iz + 1,islope) - tsoil_old(ig,iz,islope))/(mlayer_PEM(iz) - mlayer_PEM(iz - 1)) 322 tsoil(ig,isoil,islope) = a*(z - mlayer_PEM(iz - 1)) + tsoil_old(ig,iz,islope) 323 enddo 324 325 ! The "new deepest layers" temperature is set by solving the steady heat equation 326 call solve_steady_heat(nsoil - ishift + 1,mlayer_PEM(ishift - 1:),layer_PEM(ishift:),mthermdiff_PEM(ig,ishift - 1:),thermdiff_PEM(ig,ishift:),tsoil(ig,ishift,islope),fluxgeo,tsoil(ig,ishift:,islope)) 327 endif 328 endif 329 enddo 330 enddo 331 332 END SUBROUTINE shift_tsoil2surf 333 218 334 END MODULE compute_soiltemp_mod -
trunk/LMDZ.COMMON/libf/evolution/evol_ice_mod.F90
r3498 r3553 7 7 !======================================================================= 8 8 9 SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice )9 SUBROUTINE evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) 10 10 11 11 use time_evol_mod, only: dt 12 use glaciers_mod, only: rho_co2ice 12 13 13 14 implicit none … … 18 19 ! 19 20 !======================================================================= 20 ! 21 ! 22 ! 21 ! arguments: 22 ! ---------- 23 ! INPUT 23 24 integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes 24 25 25 ! OUTPUT 26 real, dimension(ngrid,nslope), intent(inout) :: co2_ice ! Previous and actual density of CO2 ice 27 real, dimension(ngrid,nslope), intent(inout) :: d_co2ice ! Evolution of perennial ice over one year 26 ! OUTPUT 27 real, dimension(ngrid,nslope), intent(inout) :: co2_ice ! Previous and actual density of CO2 ice 28 real, dimension(ngrid,nslope), intent(inout) :: d_co2ice ! Evolution of perennial ice over one year 29 real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m] 28 30 29 ! 30 ! 31 ! local: 32 ! ------ 31 33 real, dimension(ngrid,nslope) :: co2_ice_old ! Old density of CO2 ice 32 34 !======================================================================= … … 41 43 end where 42 44 45 zshift_surf = d_co2ice*dt/rho_co2ice 46 43 47 END SUBROUTINE evol_co2_ice 44 48 45 49 !======================================================================= 46 50 47 SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice, stopPEM)51 SUBROUTINE evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM) 48 52 49 use time_evol_mod, only: dt 50 use comslope_mod, only: subslope_dist, def_slope_mean 53 use time_evol_mod, only: dt 54 use comslope_mod, only: subslope_dist, def_slope_mean 55 use glaciers_mod, only: rho_h2oice 51 56 52 57 #ifndef CPP_STD … … 63 68 ! 64 69 !======================================================================= 65 ! 66 ! 67 ! 70 ! arguments: 71 ! ---------- 72 ! INPUT 68 73 integer, intent(in) :: ngrid, nslope ! # of grid points, # of subslopes 69 74 real, dimension(ngrid), intent(in) :: cell_area ! Area of each mesh grid (m^2) … … 71 76 real, dimension(ngrid), intent(in) :: delta_h2o_icetablesublim ! Mass of H2O that have condensed/sublimated at the ice table (kg/m^2) 72 77 73 ! OUTPUT 74 real, dimension(ngrid,nslope), intent(inout) :: h2o_ice ! Previous and actual density of h2o ice (kg/m^2) 75 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year) 76 integer, intent(inout) :: stopPEM ! Stopping criterion code 78 ! OUTPUT 79 real, dimension(ngrid,nslope), intent(inout) :: h2o_ice ! Previous and actual density of h2o ice (kg/m^2) 80 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! Evolution of perennial ice over one year (kg/m^2/year) 81 integer, intent(inout) :: stopPEM ! Stopping criterion code 82 real, dimension(ngrid,nslope), intent(out) :: zshift_surf ! Elevation shift for the surface [m] 77 83 78 ! 79 ! 84 ! local: 85 ! ------ 80 86 integer :: i, islope ! Loop variables 81 87 real :: pos_tend, neg_tend, real_coefficient, negative_part ! Variable to conserve h2o 82 real, dimension(ngrid,nslope) :: new_tend ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done88 real, dimension(ngrid,nslope) :: new_tend ! Tendencies computed in order to conserve h2o ice on the surface, only exchange between surface are done 83 89 !======================================================================= 84 90 if (ngrid /= 1) then ! Not in 1D … … 164 170 endif 165 171 172 zshift_surf = d_h2oice*dt/rho_h2oice 173 166 174 END SUBROUTINE evol_h2o_ice 167 175 -
trunk/LMDZ.COMMON/libf/evolution/glaciers_mod.F90
r3532 r3553 10 10 11 11 ! Thresholds for ice management 12 real, save :: inf_h2oice_threshold ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2] 13 real, save :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2] 14 real, save :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2] 12 real :: inf_h2oice_threshold ! To consider the amount of H2O ice as an infinite reservoir [kg.m-2] 13 real :: metam_h2oice_threshold ! To consider frost is becoming perennial H2O ice [kg.m-2] 14 real :: metam_co2ice_threshold ! To consider frost is becoming perennial CO2 ice [kg.m-2] 15 16 real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3] 17 real, parameter :: rho_h2oice = 920. ! Density of H2O ice [kg.m-3] 15 18 16 19 !======================================================================= -
trunk/LMDZ.COMMON/libf/evolution/layering_mod.F90
r3498 r3553 8 8 !======================================================================= 9 9 10 use glaciers_mod, only: rho_co2ice, rho_h2oice 10 11 11 12 implicit none … … 23 24 real, parameter :: co2ice_porosity = 0.3 ! Porosity of CO2 ice 24 25 real, parameter :: h2oice_porosity = 0.3 ! Porosity of H2O ice 25 real, parameter :: rho_co2ice = 1650. ! Density of CO2 ice [kg.m-3]26 real, parameter :: rho_h2oice = 920. ! Density of H2O ice [kg.m-3]27 26 real, parameter :: rho_dust = 2500. ! Density of dust [kg.m-3] 28 27 real, parameter :: rho_rock = 3200. ! Density of rock [kg.m-3] … … 452 451 !======================================================================= 453 452 ! Layering algorithm 454 SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,new_lag1,new_lag2,stopPEM,current1,current2) 455 456 implicit none 457 458 !---- Arguments 453 SUBROUTINE make_layering(this,d_co2ice,d_h2oice,d_dust,new_str,zshift_surf,new_lag1,new_lag2,zlag,stopPEM,current1,current2) 454 455 implicit none 456 457 !---- Arguments 458 ! Inputs 459 real, intent(in) :: d_co2ice, d_h2oice, d_dust 460 461 ! Outputs 459 462 type(layering), intent(inout) :: this 460 463 type(stratum), pointer, intent(inout) :: current1, current2 461 464 logical, intent(inout) :: new_str, new_lag1, new_lag2 462 465 integer, intent(inout) :: stopPEM 463 real, intent(in) :: d_co2ice, d_h2oice, d_dust466 real, intent(out) :: zshift_surf, zlag 464 467 465 468 !---- Local variables … … 474 477 dh_h2oice = d_h2oice/rho_h2oice 475 478 dh_dust = d_dust/rho_dust 479 zshift_surf = 0. 480 zlag = 0. 476 481 477 482 if (dh_dust < 0.) then ! Dust lifting only … … 504 509 stopPEM = 8 505 510 endif 511 zshift_surf = dh_dust + h2lift_tot 506 512 507 513 else … … 521 527 endif 522 528 endif 529 zshift_surf = thickness 523 530 524 531 !------ Condensation of CO2 ice + H2O ice … … 536 543 endif 537 544 endif 545 zshift_surf = thickness 538 546 539 547 !------ Sublimation of CO2 ice + Condensation of H2O ice … … 563 571 h2subl = h2subl_tot 564 572 h2subl_tot = 0. 565 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag1)573 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) 566 574 else if (h_co2ice_old < eps) then ! There is nothing to sublimate 567 575 ! Is the considered stratum a dust lag? … … 575 583 h2subl = h_co2ice_old 576 584 h2subl_tot = h2subl_tot - h2subl 577 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag1)585 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) 578 586 current1 => current1%down 579 587 new_lag1 = .true. … … 595 603 endif 596 604 endif 605 zshift_surf = zshift_surf + thickness 597 606 598 607 !------ Sublimation of CO2 ice + H2O ice … … 622 631 h2subl = h2subl_tot 623 632 h2subl_tot = 0. 624 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag1)633 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) 625 634 else if (h_co2ice_old < eps) then ! There is nothing to sublimate 626 635 ! Is the considered stratum a dust lag? … … 634 643 h2subl = h_co2ice_old 635 644 h2subl_tot = h2subl_tot - h2subl 636 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag1)645 call subl_co2ice_layering(this,current1,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag1,zlag) 637 646 current1 => current1%down 638 647 new_lag1 = .true. … … 666 675 h2subl = h2subl_tot 667 676 h2subl_tot = 0. 668 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag2)677 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag) 669 678 else if (h_h2oice_old < eps) then ! There is nothing to sublimate 670 679 ! Is the considered stratum a dust lag? … … 678 687 h2subl = h_h2oice_old 679 688 h2subl_tot = h2subl_tot - h2subl 680 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag2)689 call subl_h2oice_layering(this,current2,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag2,zlag) 681 690 current2 => current2%down 682 691 new_lag2 = .true. … … 698 707 !======================================================================= 699 708 ! To sublimate CO2 ice in the layering 700 SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag)709 SUBROUTINE subl_co2ice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag) 701 710 702 711 implicit none … … 706 715 type(stratum), pointer, intent(inout) :: str 707 716 logical, intent(inout) :: new_lag 717 real, intent(inout) :: zshift_surf, zlag 708 718 real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old 709 719 … … 719 729 720 730 str%thickness = str%thickness - h2subl/(1. - co2ice_porosity) - hsubl_dust 731 zshift_surf = zshift_surf - h2subl/(1. - co2ice_porosity) - hsubl_dust 721 732 if (str%thickness < tol) then 722 733 ! Remove the sublimating stratum if there is no ice anymore … … 748 759 endif 749 760 endif 761 zlag = zlag + h_lag 750 762 751 763 END SUBROUTINE subl_co2ice_layering … … 753 765 !======================================================================= 754 766 ! To sublimate H2O ice in the layering 755 SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old, new_lag)767 SUBROUTINE subl_h2oice_layering(this,str,h2subl,h_co2ice_old,h_h2oice_old,h_dust_old,zshift_surf,new_lag,zlag) 756 768 757 769 implicit none … … 761 773 type(stratum), pointer, intent(inout) :: str 762 774 logical, intent(inout) :: new_lag 775 real, intent(inout) :: zshift_surf, zlag 763 776 real, intent(in) :: h2subl, h_co2ice_old, h_h2oice_old, h_dust_old 764 777 … … 774 787 775 788 str%thickness = str%thickness - h2subl/(1. - h2oice_porosity) - hsubl_dust 789 zshift_surf = zshift_surf - h2subl/(1. - h2oice_porosity) - hsubl_dust 776 790 if (str%thickness < tol) then 777 791 ! Remove the sublimating stratum if there is no ice anymore … … 803 817 endif 804 818 endif 819 zlag = zlag + h_lag 805 820 806 821 END SUBROUTINE subl_h2oice_layering -
trunk/LMDZ.COMMON/libf/evolution/math_mod.F90
r3532 r3553 28 28 implicit none 29 29 30 integer, intent(in) :: nz ! number of layer 31 real, intent(in) :: z(nz) ! depth layer 32 real, intent(in) :: y(nz) ! function which needs to be derived 33 real, intent(in) :: y0,ybot ! boundary conditions 34 real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth 35 ! local 30 ! Inputs 31 !------- 32 integer, intent(in) :: nz ! number of layer 33 real, dimension(nz), intent(in) :: z ! depth layer 34 real, dimension(nz), intent(in) :: y ! function which needs to be derived 35 real, intent(in) :: y0, ybot ! boundary conditions 36 ! Outputs 37 !-------- 38 real, dimension(nz), intent(out) :: dzY ! derivative of y w.r.t depth 39 ! Local variables 40 !---------------- 36 41 integer :: j 37 42 real :: hm, hp, c1, c2, c3 … … 69 74 implicit none 70 75 71 integer, intent(in) :: nz 72 real, intent(in) :: z(nz), y(nz), y0, yNp1 73 real, intent(out) :: yp2(nz) 76 ! Inputs 77 !------- 78 integer, intent(in) :: nz 79 real, dimension(nz), intent(in) :: z, y 80 real, intent(in) :: y0, yNp1 81 ! Outputs 82 !-------- 83 real, dimension(nz), intent(out) :: yp2 84 ! Local variables 85 !---------------- 74 86 integer :: j 75 87 real :: hm, hp, c1, c2, c3 … … 102 114 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 103 115 104 implicit none 105 106 integer, intent(in) :: nz, j 107 real, intent(in) :: z(nz), y(nz) 116 implicit none 117 118 ! Inputs 119 !------- 120 integer, intent(in) :: nz, j 121 real, dimension(nz), intent(in) :: z, y 122 ! Outputs 123 !-------- 108 124 real, intent(out) :: dy_zj 125 ! Local viariables 126 !----------------- 109 127 real :: h1, h2, c1, c2, c3 110 128 … … 135 153 implicit none 136 154 137 integer, intent(in) :: nz, i1, i2 138 real, intent(in) :: y(nz), z(nz) 155 ! Inputs 156 !------- 157 integer, intent(in) :: nz, i1, i2 158 real, dimension(nz), intent(in) :: y, z 159 ! Outputs 160 !-------- 139 161 real, intent(out) :: integral 140 integer i 141 real dz(nz) 162 ! Local viariables 163 !----------------- 164 integer :: i 165 real, dimension(nz) :: dz 142 166 143 167 dz(1) = (z(2) - 0.)/2 … … 163 187 implicit none 164 188 189 ! Inputs 190 !------- 165 191 real, intent(in) :: y1, y2 ! difference between surface water density and at depth [kg/m^3] 166 192 real, intent(in) :: z1, z2 ! depth [m] 167 real, intent(out) :: zr ! depth at which we have zero 193 ! Outputs 194 !-------- 195 real, intent(out) :: zr ! depth at which we have zero 168 196 169 197 zr = (y1*z2 - y2*z1)/(y1 - y2) … … 171 199 END SUBROUTINE findroot 172 200 201 !======================================================================= 202 203 SUBROUTINE solve_tridiag(a,b,c,d,n,x,error) 204 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 205 !!! 206 !!! Purpose: Solve a tridiagonal system Ax = d using the Thomas' algorithm 207 !!! a: sub-diagonal 208 !!! b: main diagonal 209 !!! c: super-diagonal 210 !!! d: right-hand side 211 !!! x: solution 212 !!! 213 !!! Author: JBC 214 !!! 215 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 216 217 implicit none 218 219 ! Inputs 220 !------- 221 integer, intent(in) :: n 222 real, dimension(n), intent(in) :: b, d 223 real, dimension(n - 1), intent(in) :: a, c 224 ! Outputs 225 !-------- 226 real, dimension(n), intent(out) :: x 227 integer, intent(out) :: error 228 ! Local viariables 229 !----------------- 230 integer :: i 231 real :: m 232 real, dimension(n) :: cp, dp 233 234 ! Check stability: diagonally dominant condition 235 error = 0 236 if (abs(b(1)) < abs(c(1))) then 237 error = 1 238 return 239 endif 240 do i = 2,n - 1 241 if (abs(b(i)) < abs(a(i - 1)) + abs(c(i))) then 242 error = 1 243 return 244 endif 245 enddo 246 if (abs(b(n)) < abs(a(n - 1))) then 247 error = 1 248 return 249 endif 250 251 ! Initialization 252 cp(1) = c(1)/b(1) 253 dp(1) = d(1)/b(1) 254 255 ! Forward sweep 256 do i = 2,n - 1 257 m = b(i) - a(i - 1)*cp(i - 1) 258 cp(i) = c(i)/m 259 dp(i) = (d(i) - a(i - 1)*dp(i - 1))/m 260 enddo 261 m = b(n) - a(n - 1)*cp(n - 1) 262 dp(n) = (d(n) - a(n - 1)*dp(n - 1))/m 263 264 ! Backward substitution 265 x(n) = dp(n) 266 do i = n - 1,1,-1 267 x(i) = dp(i) - cp(i)*x(i + 1) 268 enddo 269 270 END SUBROUTINE solve_tridiag 271 272 !======================================================================= 273 274 SUBROUTINE solve_steady_heat(n,z,mz,kappa,mkappa,T_left,q_right,T) 275 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 276 !!! 277 !!! Purpose: Solve 1D steady-state heat equation with space-dependent diffusivity 278 !!! Left boudary condition : prescribed temperature T_left 279 !!! Right boudary condition: prescribed thermal flux q_right 280 !!! 281 !!! z : grid points 282 !!! mz : mid-grid points 283 !!! kappa : thermal diffusivity at grid points 284 !!! mkappa: thermal diffusivity at mid-grid points 285 !!! 286 !!! Author: JBC 287 !!! 288 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 289 290 use abort_pem_mod, only: abort_pem 291 292 implicit none 293 294 ! Inputs 295 !------- 296 integer, intent(in) :: n 297 real, dimension(n), intent(in) :: z, kappa 298 real, dimension(n - 1), intent(in) :: mz, mkappa 299 real, intent(in) :: T_left, q_right 300 ! Outputs 301 !-------- 302 real, dimension(n), intent(out) :: T 303 ! Local viariables 304 !----------------- 305 integer :: i, error 306 real, dimension(n) :: b, d 307 real, dimension(n - 1) :: a, c 308 309 ! Initialization 310 a = 0.; b = 0.; c = 0.; d = 0. 311 312 ! Left boundary condition (Dirichlet: prescribed temperature) 313 b(1) = 1. 314 d(1) = T_left 315 316 ! Internal points 317 do i = 2,n - 1 318 a(i - 1) = -mkappa(i - 1)/((mz(i) - mz(i - 1))*(z(i) - z(i - 1))) 319 c(i) = -mkappa(i)/((mz(i) - mz(i - 1))*(z(i + 1) - z(i))) 320 b(i) = -(a(i - 1) + c(i)) 321 enddo 322 323 ! Right boundary condition (Neumann: prescribed temperature) 324 a(n - 1) = kappa(n - 1)/(z(n) - z(n - 1)) 325 b(n) = -kappa(n)/(z(n) - z(n - 1)) 326 d(n) = q_right 327 328 ! Solve the tridiagonal linear system with the Thomas' algorithm 329 call solve_tridiag(a,b,c,d,n,T,error) 330 if (error /= 0) call abort_pem("solve_steady_heat","Unstable solving!",1) 331 332 END SUBROUTINE solve_steady_heat 333 173 334 end module math_mod -
trunk/LMDZ.COMMON/libf/evolution/pem.F90
r3543 r3553 44 44 use constants_marspem_mod, only: alpha_clap_co2, beta_clap_co2, alpha_clap_h2o, beta_clap_h2o, m_co2, m_noco2 45 45 use evol_ice_mod, only: evol_co2_ice, evol_h2o_ice 46 use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, &46 use comsoil_h_PEM, only: soil_pem, ini_comsoil_h_PEM, end_comsoil_h_PEM, nsoilmx_PEM, mlayer_PEM, & 47 47 TI_PEM, & ! soil thermal inertia 48 48 tsoil_PEM, layer_PEM, & ! Soil temp, number of subsurface layers, soil mid layer depths … … 50 50 use adsorption_mod, only: regolith_adsorption, adsorption_pem, & ! Bool to check if adsorption, main subroutine 51 51 ini_adsorption_h_PEM, end_adsorption_h_PEM, & ! Allocate arrays 52 co2_adsorb ded_phys, h2o_adsorbded_phys ! Mass of co2 and h2O adsorbded52 co2_adsorbed_phys, h2o_adsorbed_phys ! Mass of co2 and h2O adsorbed 53 53 use time_evol_mod, only: dt, evol_orbit_pem, Max_iter_pem, convert_years, year_bp_ini 54 54 use orbit_param_criterion_mod, only: orbit_param_criterion … … 65 65 use pemetat0_mod, only: pemetat0 66 66 use read_data_PCM_mod, only: read_data_PCM 67 use recomp_tend_co2_slope_mod, only: recomp_tend_co2 _slope68 use compute_soiltemp_mod, only: compute_tsoil_pem 67 use recomp_tend_co2_slope_mod, only: recomp_tend_co2 68 use compute_soiltemp_mod, only: compute_tsoil_pem, shift_tsoil2surf 69 69 use writediagpem_mod, only: writediagpem, writediagsoilpem 70 70 use co2condens_mod, only: CO2cond_ps … … 217 217 real, dimension(:,:,:), allocatable :: watersoil_density_PEM_avg ! Physic x Soil x Slopes, water soil density, yearly averaged [kg/m^3] 218 218 real, dimension(:,:), allocatable :: Tsurfavg_before_saved ! Surface temperature saved from previous time step [K] 219 real, dimension(:), allocatable :: delta_co2_adsorb ded ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2]220 real, dimension(:), allocatable :: delta_h2o_adsorb ded ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2]221 real :: totmassco2_adsorb ded ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg]222 real :: totmassh2o_adsorb ded ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg]219 real, dimension(:), allocatable :: delta_co2_adsorbed ! Physics: quantity of CO2 that is exchanged because of adsorption / desorption [kg/m^2] 220 real, dimension(:), allocatable :: delta_h2o_adsorbed ! Physics: quantity of H2O that is exchanged because of adsorption / desorption [kg/m^2] 221 real :: totmassco2_adsorbed ! Total mass of CO2 that is exchanged because of adsorption / desoprtion over the planets [kg] 222 real :: totmassh2o_adsorbed ! Total mass of H2O that is exchanged because of adsorption / desoprtion over the planets [kg] 223 223 logical :: bool_sublim ! logical to check if there is sublimation or not 224 224 logical, dimension(:,:), allocatable :: co2ice_disappeared ! logical to check if a co2 ice reservoir already disappeared at a previous timestep … … 228 228 real, dimension(:), allocatable :: porefill ! Pore filling (output) to compute the dynamic ice table 229 229 real :: ssi_depth ! Ice table depth (output) to compute the dynamic ice table 230 real, dimension(:,:), allocatable :: zshift_surf ! Elevation shift for the surface [m] 231 real, dimension(:,:), allocatable :: zlag ! Newly built lag thickness [m] 230 232 231 233 ! Some variables for the PEM run … … 541 543 allocate(vmr_co2_PCM(ngrid,timelen)) 542 544 allocate(ps_timeseries(ngrid,timelen)) 543 allocate(min_co2_ice(ngrid,nslope,2)) 544 allocate(min_h2o_ice(ngrid,nslope,2)) 545 allocate(min_co2_ice(ngrid,nslope,2),min_h2o_ice(ngrid,nslope,2)) 545 546 allocate(tsurf_avg_yr1(ngrid,nslope)) 546 547 allocate(tsurf_avg(ngrid,nslope)) 547 548 allocate(tsurf_PCM_timeseries(ngrid,nslope,timelen)) 548 549 allocate(tsoil_anom(ngrid,nsoilmx,nslope,timelen)) 549 allocate( q_co2_PEM_phys(ngrid,timelen))550 allocate(q_ h2o_PEM_phys(ngrid,timelen))550 allocate(zshift_surf(ngrid,nslope),zlag(ngrid,nslope)) 551 allocate(q_co2_PEM_phys(ngrid,timelen),q_h2o_PEM_phys(ngrid,timelen)) 551 552 allocate(watersurf_density_avg(ngrid,nslope)) 552 553 allocate(watersoil_density_timeseries(ngrid,nsoilmx,nslope,timelen)) … … 554 555 allocate(tsoil_phys_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 555 556 allocate(watersoil_density_PEM_timeseries(ngrid,nsoilmx_PEM,nslope,timelen)) 556 allocate(delta_co2_adsorb ded(ngrid))557 allocate(delta_co2_adsorbed(ngrid)) 557 558 allocate(co2ice_disappeared(ngrid,nslope)) 558 allocate(icetable_thickness_old(ngrid,nslope)) 559 allocate(ice_porefilling_old(ngrid,nsoilmx_PEM,nslope)) 560 allocate(delta_h2o_icetablesublim(ngrid)) 561 allocate(delta_h2o_adsorbded(ngrid)) 559 allocate(icetable_thickness_old(ngrid,nslope),ice_porefilling_old(ngrid,nsoilmx_PEM,nslope)) 560 allocate(delta_h2o_icetablesublim(ngrid),delta_h2o_adsorbed(ngrid)) 562 561 allocate(vmr_co2_PEM_phys(ngrid,timelen)) 563 562 … … 649 648 icetable_thickness,ice_porefilling,tsurf_avg_yr1,tsurf_avg,q_co2_PEM_phys,q_h2o_PEM_phys, & 650 649 ps_timeseries,tsoil_phys_PEM_timeseries,d_h2oice,d_co2ice,co2_ice,h2o_ice, & 651 global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorb ded_phys, &652 delta_co2_adsorb ded,h2o_adsorbded_phys,delta_h2o_adsorbded,stratif)650 global_avg_press_PCM,watersurf_density_avg,watersoil_density_PEM_avg,co2_adsorbed_phys, & 651 delta_co2_adsorbed,h2o_adsorbed_phys,delta_h2o_adsorbed,stratif) 653 652 654 653 ! We save the places where h2o ice is sublimating … … 680 679 681 680 if (adsorption_pem) then 682 totmassco2_adsorb ded = 0.683 totmassh2o_adsorb ded = 0.681 totmassco2_adsorbed = 0. 682 totmassh2o_adsorbed = 0. 684 683 do ig = 1,ngrid 685 684 do islope = 1,nslope 686 685 do l = 1,nsoilmx_PEM - 1 687 686 if (l==1) then 688 totmassco2_adsorb ded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &687 totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & 689 688 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 690 totmassh2o_adsorb ded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &689 totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & 691 690 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 692 691 else 693 totmassco2_adsorb ded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &692 totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & 694 693 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 695 totmassh2o_adsorb ded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* &694 totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l-1))* & 696 695 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 697 696 endif … … 699 698 enddo 700 699 enddo 701 write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorb ded702 write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorb ded700 write(*,*) "Tot mass of CO2 in the regolith =", totmassco2_adsorbed 701 write(*,*) "Tot mass of H2O in the regolith =", totmassh2o_adsorbed 703 702 endif ! adsorption 704 703 deallocate(tsurf_avg_yr1) … … 750 749 if (adsorption_pem) then 751 750 do i = 1,ngrid 752 global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorb ded(i)/Total_surface751 global_avg_press_new = global_avg_press_new - g*cell_area(i)*delta_co2_adsorbed(i)/Total_surface 753 752 enddo 754 753 endif … … 819 818 ! II_b Evolution of ice 820 819 !------------------------ 821 call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorb ded,delta_h2o_icetablesublim,h2o_ice,d_h2oice,stopPEM)822 call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice )820 call evol_h2o_ice(ngrid,nslope,cell_area,delta_h2o_adsorbed,delta_h2o_icetablesublim,h2o_ice,d_h2oice,zshift_surf,stopPEM) 821 call evol_co2_ice(ngrid,nslope,co2_ice,d_co2ice,zshift_surf) 823 822 if (layering_algo) then 824 823 do islope = 1,nslope 825 824 do ig = 1,ngrid 826 call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope), new_lag1(ig,islope),new_lag2(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p)825 call make_layering(stratif(ig,islope),d_co2ice(ig,islope),d_h2oice(ig,islope),d_dust,new_str(ig,islope),zshift_surf(ig,islope),new_lag1(ig,islope),new_lag2(ig,islope),zlag(ig,islope),stopPEM,current1(ig,islope)%p,current2(ig,islope)%p) 827 826 enddo 828 827 enddo 828 else 829 zlag = 0. 829 830 endif 830 831 … … 887 888 tsurf_avg(ig,islope) = ave/timelen 888 889 endif 890 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! SUBLIMATION??? !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 889 891 enddo 890 892 enddo … … 901 903 902 904 if (soil_pem) then 903 904 ! II_d.2 Update soil temperature 905 ! II_d.2 Shifting soil temperature to surface 906 call shift_tsoil2surf(ngrid,nsoilmx_PEM,nslope,zshift_surf,zlag,tsurf_avg,tsoil_PEM) 907 908 ! II_d.3 Update soil temperature 905 909 write(*,*)"Updating soil temperature" 906 910 allocate(Tsoil_locslope(ngrid,nsoilmx_PEM)) … … 926 930 deallocate(Tsoil_locslope) 927 931 928 ! II_d. 3Update the ice table932 ! II_d.4 Update the ice table 929 933 if (icetable_equilibrium) then 930 934 write(*,*) "Compute ice table (equilibrium method)" … … 947 951 endif 948 952 949 ! II_d. 4Update the soil thermal properties953 ! II_d.5 Update the soil thermal properties 950 954 call update_soil_thermalproperties(ngrid,nslope,nsoilmx_PEM,d_h2oice,h2o_ice,global_avg_press_new,icetable_depth,icetable_thickness,ice_porefilling,icetable_equilibrium,icetable_dynamic,TI_PEM) 951 955 952 ! II_d. 5Update the mass of the regolith adsorbed956 ! II_d.6 Update the mass of the regolith adsorbed 953 957 if (adsorption_pem) then 954 958 call regolith_adsorption(ngrid,nslope,nsoilmx_PEM,timelen,d_h2oice,d_co2ice, & 955 959 h2o_ice,co2_ice,tsoil_PEM,TI_PEM,ps_timeseries,q_co2_PEM_phys,q_h2o_PEM_phys, & 956 h2o_adsorb ded_phys,delta_h2o_adsorbded,co2_adsorbded_phys,delta_co2_adsorbded)957 958 totmassco2_adsorb ded = 0.959 totmassh2o_adsorb ded = 0.960 h2o_adsorbed_phys,delta_h2o_adsorbed,co2_adsorbed_phys,delta_co2_adsorbed) 961 962 totmassco2_adsorbed = 0. 963 totmassh2o_adsorbed = 0. 960 964 do ig = 1,ngrid 961 965 do islope = 1,nslope 962 966 do l = 1,nsoilmx_PEM 963 967 if (l == 1) then 964 totmassco2_adsorb ded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &968 totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & 965 969 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 966 totmassh2o_adsorb ded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l))* &970 totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l))* & 967 971 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 968 972 else 969 totmassco2_adsorb ded = totmassco2_adsorbded + co2_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &973 totmassco2_adsorbed = totmassco2_adsorbed + co2_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & 970 974 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 971 totmassh2o_adsorb ded = totmassh2o_adsorbded + h2o_adsorbded_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* &975 totmassh2o_adsorbed = totmassh2o_adsorbed + h2o_adsorbed_phys(ig,l,islope)*(layer_PEM(l) - layer_PEM(l - 1))* & 972 976 subslope_dist(ig,islope)/cos(pi*def_slope_mean(islope)/180.)*cell_area(ig) 973 977 endif … … 975 979 enddo 976 980 enddo 977 write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorb ded978 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorb ded981 write(*,*) "Tot mass of CO2 in the regolith=", totmassco2_adsorbed 982 write(*,*) "Tot mass of H2O in the regolith=", totmassh2o_adsorbed 979 983 endif 980 984 endif !soil_pem … … 1005 1009 if (icetable_dynamic) call writediagsoilpem(ngrid,'ice_porefilling'//str2,'ice pore filling','-',3,ice_porefilling(:,:,islope)) 1006 1010 if (adsorption_pem) then 1007 call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorb ded_phys(:,:,islope))1008 call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorb ded_phys(:,:,islope))1011 call writediagsoilpem(ngrid,'co2_ads_slope'//str2,'co2_ads','K',3,co2_adsorbed_phys(:,:,islope)) 1012 call writediagsoilpem(ngrid,'h2o_ads_slope'//str2,'h2o_ads','K',3,h2o_adsorbed_phys(:,:,islope)) 1009 1013 endif 1010 1014 endif … … 1015 1019 ! II_f Update the tendencies 1016 1020 !------------------------ 1017 call recomp_tend_co2 _slope(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new)1021 call recomp_tend_co2(ngrid,nslope,timelen,d_co2ice,d_co2ice_ini,co2_ice,emis,vmr_co2_PCM,vmr_co2_PEM_phys,ps_timeseries,global_avg_press_PCM,global_avg_press_new) 1018 1022 1019 1023 !------------------------ … … 1021 1025 ! II_g Checking the stopping criterion 1022 1026 !------------------------ 1023 1024 1027 write(*,*) "Checking the stopping criteria..." 1025 1028 call stopping_crit_h2o_ice(cell_area,h2oice_ini_surf,ini_h2oice_sublim,h2o_ice,stopPEM,ngrid) … … 1210 1213 call pemdem0("restartpem.nc",longitude,latitude,cell_area,ngrid,nslope,def_slope,subslope_dist) 1211 1214 call pemdem1("restartpem.nc",i_myear,nsoilmx_PEM,ngrid,nslope,tsoil_PEM,TI_PEM,icetable_depth,icetable_thickness,ice_porefilling, & 1212 co2_adsorb ded_phys,h2o_adsorbded_phys,h2o_ice,stratif)1215 co2_adsorbed_phys,h2o_adsorbed_phys,h2o_ice,stratif) 1213 1216 write(*,*) "restartpem.nc has been written" 1214 1217 … … 1231 1234 deallocate(watersurf_density_avg,watersoil_density_timeseries,Tsurfavg_before_saved) 1232 1235 deallocate(tsoil_phys_PEM_timeseries,watersoil_density_PEM_timeseries,watersoil_density_PEM_avg) 1233 deallocate(delta_co2_adsorb ded,delta_h2o_adsorbded,vmr_co2_PEM_phys,delta_h2o_icetablesublim)1234 deallocate(icetable_thickness_old,ice_porefilling_old )1236 deallocate(delta_co2_adsorbed,delta_h2o_adsorbed,vmr_co2_PEM_phys,delta_h2o_icetablesublim) 1237 deallocate(icetable_thickness_old,ice_porefilling_old,zshift_surf,zlag) 1235 1238 deallocate(is_co2ice_ini,co2ice_disappeared,ini_co2ice_sublim,ini_h2oice_sublim,stratif) 1236 1239 !----------------------------- END OUTPUT ------------------------------ -
trunk/LMDZ.COMMON/libf/evolution/recomp_orb_param_mod.F90
r3498 r3553 65 65 write(*,*) 'recomp_orb_param, Old ecc. =', e_elips 66 66 write(*,*) 'recomp_orb_param, Old Lsp =', Lsp 67 write(*,*) 'recomp_orb_param, New year =', Year 67 68 68 69 found_year = .false. … … 116 117 #endif 117 118 118 write(*,*) 'recomp_orb_param, New year =', Year119 119 write(*,*) 'recomp_orb_param, New obl. =', obliquit 120 120 write(*,*) 'recomp_orb_param, New ecc. =', e_elips -
trunk/LMDZ.COMMON/libf/evolution/recomp_tend_co2_slope_mod.F90
r3532 r3553 7 7 !======================================================================= 8 8 9 SUBROUTINE recomp_tend_co2 _slope(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_phys_ini,co2ice_slope,emissivity_slope, &10 vmr_co2_PCM,vmr_co2_pem,ps_PCM_2,global_avg_press_PCM,global_avg_press_new)9 SUBROUTINE recomp_tend_co2(ngrid,nslope,timelen,d_co2ice_phys,d_co2ice_ini,co2ice,emissivity, & 10 vmr_co2_PCM,vmr_co2_PEM,ps_PCM_2,global_avg_press_PCM,global_avg_press_new) 11 11 12 12 use constants_marspem_mod, only : alpha_clap_co2, beta_clap_co2, sigmaB, Lco2, sols_per_my, sec_per_sol … … 24 24 ! INPUT 25 25 integer, intent(in) :: timelen, ngrid, nslope 26 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM 27 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_ pem! physical point field: Volume mixing ratio of co2 in the first layer28 real, dimension(ngrid,timelen), intent(in) :: ps_PCM_2 29 real, intent(in) :: global_avg_press_PCM 30 real, intent(in) :: global_avg_press_new 31 real, dimension(ngrid,nslope), intent(in) :: d_co2ice_ phys_ini! physical point field: Evolution of perennial ice over one year32 real, dimension(ngrid,nslope), intent(in) :: co2ice _slope! CO2 ice per mesh and sub-grid slope (kg/m^2)33 real, dimension(ngrid,nslope), intent(in) :: emissivity _slope! Emissivity per mesh and sub-grid slope(1)26 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PCM ! physical point field: Volume mixing ratio of co2 in the first layer 27 real, dimension(ngrid,timelen), intent(in) :: vmr_co2_PEM ! physical point field: Volume mixing ratio of co2 in the first layer 28 real, dimension(ngrid,timelen), intent(in) :: ps_PCM_2 ! physical point field: Surface pressure in the PCM 29 real, intent(in) :: global_avg_press_PCM ! global averaged pressure at previous timestep 30 real, intent(in) :: global_avg_press_new ! global averaged pressure at current timestep 31 real, dimension(ngrid,nslope), intent(in) :: d_co2ice_ini ! physical point field: Evolution of perennial ice over one year 32 real, dimension(ngrid,nslope), intent(in) :: co2ice ! CO2 ice per mesh and sub-grid slope (kg/m^2) 33 real, dimension(ngrid,nslope), intent(in) :: emissivity ! Emissivity per mesh and sub-grid slope(1) 34 34 ! OUTPUT 35 real, dimension(ngrid,nslope), intent(inout) :: 35 real, dimension(ngrid,nslope), intent(inout) :: d_co2ice_phys ! physical point field: Evolution of perennial ice over one year 36 36 37 37 ! local: … … 45 45 do i = 1,ngrid 46 46 do islope = 1,nslope 47 coef = sols_per_my*sec_per_sol*emissivity _slope(i,islope)*sigmaB/Lco247 coef = sols_per_my*sec_per_sol*emissivity(i,islope)*sigmaB/Lco2 48 48 ave = 0. 49 if (co2ice _slope(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then50 do t =1,timelen49 if (co2ice(i,islope) > 1.e-4 .and. abs(d_co2ice_phys(i,islope)) > 1.e-5) then 50 do t = 1,timelen 51 51 ave = ave + (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PCM(i,t)*ps_PCM_2(i,t)/100.)))**4 & 52 - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_ pem(i,t)*ps_PCM_2(i,t)*(global_avg_press_new/global_avg_press_PCM)/100.)))**452 - (beta_clap_co2/(alpha_clap_co2-log(vmr_co2_PEM(i,t)*ps_PCM_2(i,t)*(global_avg_press_new/global_avg_press_PCM)/100.)))**4 53 53 enddo 54 if (ave < 1 e-4) ave = 0.55 d_co2ice_phys(i,islope) = d_co2ice_ phys_ini(i,islope) - coef*ave/timelen54 if (ave < 1.e-4) ave = 0. 55 d_co2ice_phys(i,islope) = d_co2ice_ini(i,islope) - coef*ave/timelen 56 56 endif 57 57 enddo 58 58 enddo 59 59 60 END SUBROUTINE recomp_tend_co2_slope 60 END SUBROUTINE recomp_tend_co2 61 !======================================================================= 62 63 SUBROUTINE recomp_tend_h2o(ngrid,nslope,timelen,d_h2oice,PCM_temp,PEM_temp) 64 65 implicit none 66 67 !======================================================================= 68 ! 69 ! Routine that compute the evolution of the tendencie for h2o ice 70 ! 71 !======================================================================= 72 73 ! arguments: 74 ! ---------- 75 ! INPUT 76 integer, intent(in) :: timelen, ngrid, nslope 77 real, dimension(:), intent(in) :: PCM_temp, PEM_temp 78 79 ! OUTPUT 80 real, dimension(ngrid,nslope), intent(inout) :: d_h2oice ! physical point field: Evolution of perennial ice over one year 81 82 ! local: 83 ! ------ 84 real :: coef, ave, Rz_old, Rz_new, R_dec, soil_psv_old, soil_psv_new, hum_dec, h2oice_depth_new, h2oice_depth_old 85 86 write(*,*) "Update of the H2O tendency due to lag layer" 87 88 ! Flux correction due to lag layer 89 !~ Rz_old = h2oice_depth_old*0.0325/4.e-4 ! resistance from PCM 90 !~ Rz_new = h2oice_depth_new*0.0325/4.e-4 ! new resistance based on new depth 91 !~ R_dec = (1./Rz_old)/(1./Rz_new) ! decrease because of resistance 92 !~ soil_psv_old = psv(max(PCM_temp(h2oice_depth_old))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the GCM run temperature at the old ice location 93 !~ soil_psv_new = psv(max(PEM_temp(h2oice_depth_new))) ! the maxmimum annual mean saturation vapor pressure at the temperature of the PEM run temperature at the new ice location 94 !~ hum_dec = soil_psv_old/soil_psv_new ! decrease because of lower water vapor pressure at the new depth 95 !~ d_h2oice = d_h2oice*R_dec*hum_dec ! decrease of flux 96 97 END SUBROUTINE recomp_tend_h2o 61 98 62 99 END MODULE recomp_tend_co2_slope_mod
Note: See TracChangeset
for help on using the changeset viewer.